cchkgt.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 468 行 · 第 1/2 页

                     A( 1 ) = ZERO
                  END IF
               ELSE IF( IMAT.EQ.9 ) THEN
                  IZERO = N
                  Z( 1 ) = A( 3*N-2 )
                  Z( 2 ) = A( 2*N-1 )
                  A( 3*N-2 ) = ZERO
                  A( 2*N-1 ) = ZERO
               ELSE
                  IZERO = ( N+1 ) / 2
                  DO 20 I = IZERO, N - 1
                     A( 2*N-2+I ) = ZERO
                     A( N-1+I ) = ZERO
                     A( I ) = ZERO
   20             CONTINUE
                  A( 3*N-2 ) = ZERO
                  A( 2*N-1 ) = ZERO
               END IF
            END IF
*
*+    TEST 1
*           Factor A as L*U and compute the ratio
*              norm(L*U - A) / (n * norm(A) * EPS )
*
            CALL CCOPY( N+2*M, A, 1, AF, 1 )
            SRNAMT = 'CGTTRF'
            CALL CGTTRF( N, AF, AF( M+1 ), AF( N+M+1 ), AF( N+2*M+1 ),
     $                   IWORK, INFO )
*
*           Check error code from CGTTRF.
*
            IF( INFO.NE.IZERO )
     $         CALL ALAERH( PATH, 'CGTTRF', INFO, IZERO, ' ', N, N, 1,
     $                      1, -1, IMAT, NFAIL, NERRS, NOUT )
            TRFCON = INFO.NE.0
*
            CALL CGTT01( N, A, A( M+1 ), A( N+M+1 ), AF, AF( M+1 ),
     $                   AF( N+M+1 ), AF( N+2*M+1 ), IWORK, WORK, LDA,
     $                   RWORK, RESULT( 1 ) )
*
*           Print the test ratio if it is .GE. THRESH.
*
            IF( RESULT( 1 ).GE.THRESH ) THEN
               IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $            CALL ALAHD( NOUT, PATH )
               WRITE( NOUT, FMT = 9999 )N, IMAT, 1, RESULT( 1 )
               NFAIL = NFAIL + 1
            END IF
            NRUN = NRUN + 1
*
            DO 50 ITRAN = 1, 2
               TRANS = TRANSS( ITRAN )
               IF( ITRAN.EQ.1 ) THEN
                  NORM = 'O'
               ELSE
                  NORM = 'I'
               END IF
               ANORM = CLANGT( NORM, N, A, A( M+1 ), A( N+M+1 ) )
*
               IF( .NOT.TRFCON ) THEN
*
*                 Use CGTTRS to solve for one column at a time of
*                 inv(A), computing the maximum column sum as we go.
*
                  AINVNM = ZERO
                  DO 40 I = 1, N
                     DO 30 J = 1, N
                        X( J ) = ZERO
   30                CONTINUE
                     X( I ) = ONE
                     CALL CGTTRS( TRANS, N, 1, AF, AF( M+1 ),
     $                            AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
     $                            LDA, INFO )
                     AINVNM = MAX( AINVNM, SCASUM( N, X, 1 ) )
   40             CONTINUE
*
*                 Compute RCONDC = 1 / (norm(A) * norm(inv(A))
*
                  IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
                     RCONDC = ONE
                  ELSE
                     RCONDC = ( ONE / ANORM ) / AINVNM
                  END IF
                  IF( ITRAN.EQ.1 ) THEN
                     RCONDO = RCONDC
                  ELSE
                     RCONDI = RCONDC
                  END IF
               ELSE
                  RCONDC = ZERO
               END IF
*
*+    TEST 7
*              Estimate the reciprocal of the condition number of the
*              matrix.
*
               SRNAMT = 'CGTCON'
               CALL CGTCON( NORM, N, AF, AF( M+1 ), AF( N+M+1 ),
     $                      AF( N+2*M+1 ), IWORK, ANORM, RCOND, WORK,
     $                      INFO )
*
*              Check error code from CGTCON.
*
               IF( INFO.NE.0 )
     $            CALL ALAERH( PATH, 'CGTCON', INFO, 0, NORM, N, N, -1,
     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
*
               RESULT( 7 ) = SGET06( RCOND, RCONDC )
*
*              Print the test ratio if it is .GE. THRESH.
*
               IF( RESULT( 7 ).GE.THRESH ) THEN
                  IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $               CALL ALAHD( NOUT, PATH )
                  WRITE( NOUT, FMT = 9997 )NORM, N, IMAT, 7,
     $               RESULT( 7 )
                  NFAIL = NFAIL + 1
               END IF
               NRUN = NRUN + 1
   50       CONTINUE
*
*           Skip the remaining tests if the matrix is singular.
*
            IF( TRFCON )
     $         GO TO 100
*
            DO 90 IRHS = 1, NNS
               NRHS = NSVAL( IRHS )
*
*              Generate NRHS random solution vectors.
*
               IX = 1
               DO 60 J = 1, NRHS
                  CALL CLARNV( 2, ISEED, N, XACT( IX ) )
                  IX = IX + LDA
   60          CONTINUE
*
               DO 80 ITRAN = 1, 3
                  TRANS = TRANSS( ITRAN )
                  IF( ITRAN.EQ.1 ) THEN
                     RCONDC = RCONDO
                  ELSE
                     RCONDC = RCONDI
                  END IF
*
*                 Set the right hand side.
*
                  CALL CLAGTM( TRANS, N, NRHS, ONE, A,
     $                         A( M+1 ), A( N+M+1 ), XACT, LDA,
     $                         ZERO, B, LDA )
*
*+    TEST 2
*              Solve op(A) * X = B and compute the residual.
*
                  CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
                  SRNAMT = 'CGTTRS'
                  CALL CGTTRS( TRANS, N, NRHS, AF, AF( M+1 ),
     $                         AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
     $                         LDA, INFO )
*
*              Check error code from CGTTRS.
*
                  IF( INFO.NE.0 )
     $               CALL ALAERH( PATH, 'CGTTRS', INFO, 0, TRANS, N, N,
     $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
     $                            NOUT )
*
                  CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
                  CALL CGTT02( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
     $                         X, LDA, WORK, LDA, RWORK, RESULT( 2 ) )
*
*+    TEST 3
*              Check solution from generated exact solution.
*
                  CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
     $                         RESULT( 3 ) )
*
*+    TESTS 4, 5, and 6
*              Use iterative refinement to improve the solution.
*
                  SRNAMT = 'CGTRFS'
                  CALL CGTRFS( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
     $                         AF, AF( M+1 ), AF( N+M+1 ),
     $                         AF( N+2*M+1 ), IWORK, B, LDA, X, LDA,
     $                         RWORK, RWORK( NRHS+1 ), WORK,
     $                         RWORK( 2*NRHS+1 ), INFO )
*
*              Check error code from CGTRFS.
*
                  IF( INFO.NE.0 )
     $               CALL ALAERH( PATH, 'CGTRFS', INFO, 0, TRANS, N, N,
     $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
     $                            NOUT )
*
                  CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
     $                         RESULT( 4 ) )
                  CALL CGTT05( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
     $                         B, LDA, X, LDA, XACT, LDA, RWORK,
     $                         RWORK( NRHS+1 ), RESULT( 5 ) )
*
*              Print information about the tests that did not pass the
*              threshold.
*
                  DO 70 K = 2, 6
                     IF( RESULT( K ).GE.THRESH ) THEN
                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                     CALL ALAHD( NOUT, PATH )
                        WRITE( NOUT, FMT = 9998 )TRANS, N, NRHS, IMAT,
     $                     K, RESULT( K )
                        NFAIL = NFAIL + 1
                     END IF
   70             CONTINUE
                  NRUN = NRUN + 5
   80          CONTINUE
   90       CONTINUE
  100    CONTINUE
  110 CONTINUE
*
*     Print a summary of the results.
*
      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
 9999 FORMAT( 12X, 'N =', I5, ',', 10X, ' type ', I2, ', test(', I2,
     $      ') = ', G12.5 )
 9998 FORMAT( ' TRANS=''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
     $      I2, ', test(', I2, ') = ', G12.5 )
 9997 FORMAT( ' NORM =''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
     $      ', test(', I2, ') = ', G12.5 )
      RETURN
*
*     End of CCHKGT
*
      END