zchkgb.f

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

*                    the matrix to test that INFO is returned correctly.
*
                     IZERO = 0
                     IF( ZEROT ) THEN
                        IF( IMAT.EQ.2 ) THEN
                           IZERO = 1
                        ELSE IF( IMAT.EQ.3 ) THEN
                           IZERO = MIN( M, N )
                        ELSE
                           IZERO = MIN( M, N ) / 2 + 1
                        END IF
                        IOFF = ( IZERO-1 )*LDA
                        IF( IMAT.LT.4 ) THEN
*
*                          Store the column to be zeroed out in B.
*
                           I1 = MAX( 1, KU+2-IZERO )
                           I2 = MIN( KL+KU+1, KU+1+( M-IZERO ) )
                           CALL ZCOPY( I2-I1+1, A( IOFF+I1 ), 1, B, 1 )
*
                           DO 30 I = I1, I2
                              A( IOFF+I ) = ZERO
   30                      CONTINUE
                        ELSE
                           DO 50 J = IZERO, N
                              DO 40 I = MAX( 1, KU+2-J ),
     $                                MIN( KL+KU+1, KU+1+( M-J ) )
                                 A( IOFF+I ) = ZERO
   40                         CONTINUE
                              IOFF = IOFF + LDA
   50                      CONTINUE
                        END IF
                     END IF
*
*                    These lines, if used in place of the calls in the
*                    loop over INB, cause the code to bomb on a Sun
*                    SPARCstation.
*
*                     ANORMO = ZLANGB( 'O', N, KL, KU, A, LDA, RWORK )
*                     ANORMI = ZLANGB( 'I', N, KL, KU, A, LDA, RWORK )
*
*                    Do for each blocksize in NBVAL
*
                     DO 110 INB = 1, NNB
                        NB = NBVAL( INB )
                        CALL XLAENV( 1, NB )
*
*                       Compute the LU factorization of the band matrix.
*
                        IF( M.GT.0 .AND. N.GT.0 )
     $                     CALL ZLACPY( 'Full', KL+KU+1, N, A, LDA,
     $                                  AFAC( KL+1 ), LDAFAC )
                        SRNAMT = 'ZGBTRF'
                        CALL ZGBTRF( M, N, KL, KU, AFAC, LDAFAC, IWORK,
     $                               INFO )
*
*                       Check error code from ZGBTRF.
*
                        IF( INFO.NE.IZERO )
     $                     CALL ALAERH( PATH, 'ZGBTRF', INFO, IZERO,
     $                                  ' ', M, N, KL, KU, NB, IMAT,
     $                                  NFAIL, NERRS, NOUT )
                        TRFCON = .FALSE.
*
*+    TEST 1
*                       Reconstruct matrix from factors and compute
*                       residual.
*
                        CALL ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC,
     $                               IWORK, WORK, RESULT( 1 ) )
*
*                       Print information about the tests so far that
*                       did not pass the threshold.
*
                        IF( RESULT( 1 ).GE.THRESH ) THEN
                           IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                        CALL ALAHD( NOUT, PATH )
                           WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB,
     $                        IMAT, 1, RESULT( 1 )
                           NFAIL = NFAIL + 1
                        END IF
                        NRUN = NRUN + 1
*
*                       Skip the remaining tests if this is not the
*                       first block size or if M .ne. N.
*
                        IF( INB.GT.1 .OR. M.NE.N )
     $                     GO TO 110
*
                        ANORMO = ZLANGB( 'O', N, KL, KU, A, LDA, RWORK )
                        ANORMI = ZLANGB( 'I', N, KL, KU, A, LDA, RWORK )
*
                        IF( INFO.EQ.0 ) THEN
*
*                          Form the inverse of A so we can get a good
*                          estimate of CNDNUM = norm(A) * norm(inv(A)).
*
                           LDB = MAX( 1, N )
                           CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ),
     $                                  DCMPLX( ONE ), WORK, LDB )
                           SRNAMT = 'ZGBTRS'
                           CALL ZGBTRS( 'No transpose', N, KL, KU, N,
     $                                  AFAC, LDAFAC, IWORK, WORK, LDB,
     $                                  INFO )
*
*                          Compute the 1-norm condition number of A.
*
                           AINVNM = ZLANGE( 'O', N, N, WORK, LDB,
     $                              RWORK )
                           IF( ANORMO.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
                              RCONDO = ONE
                           ELSE
                              RCONDO = ( ONE / ANORMO ) / AINVNM
                           END IF
*
*                          Compute the infinity-norm condition number of
*                          A.
*
                           AINVNM = ZLANGE( 'I', N, N, WORK, LDB,
     $                              RWORK )
                           IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
                              RCONDI = ONE
                           ELSE
                              RCONDI = ( ONE / ANORMI ) / AINVNM
                           END IF
                        ELSE
*
*                          Do only the condition estimate if INFO.NE.0.
*
                           TRFCON = .TRUE.
                           RCONDO = ZERO
                           RCONDI = ZERO
                        END IF
*
*                       Skip the solve tests if the matrix is singular.
*
                        IF( TRFCON )
     $                     GO TO 90
*
                        DO 80 IRHS = 1, NNS
                           NRHS = NSVAL( IRHS )
                           XTYPE = 'N'
*
                           DO 70 ITRAN = 1, NTRAN
                              TRANS = TRANSS( ITRAN )
                              IF( ITRAN.EQ.1 ) THEN
                                 RCONDC = RCONDO
                                 NORM = 'O'
                              ELSE
                                 RCONDC = RCONDI
                                 NORM = 'I'
                              END IF
*
*+    TEST 2:
*                             Solve and compute residual for A * X = B.
*
                              SRNAMT = 'ZLARHS'
                              CALL ZLARHS( PATH, XTYPE, ' ', TRANS, N,
     $                                     N, KL, KU, NRHS, A, LDA,
     $                                     XACT, LDB, B, LDB, ISEED,
     $                                     INFO )
                              XTYPE = 'C'
                              CALL ZLACPY( 'Full', N, NRHS, B, LDB, X,
     $                                     LDB )
*
                              SRNAMT = 'ZGBTRS'
                              CALL ZGBTRS( TRANS, N, KL, KU, NRHS, AFAC,
     $                                     LDAFAC, IWORK, X, LDB, INFO )
*
*                             Check error code from ZGBTRS.
*
                              IF( INFO.NE.0 )
     $                           CALL ALAERH( PATH, 'ZGBTRS', INFO, 0,
     $                                        TRANS, N, N, KL, KU, -1,
     $                                        IMAT, NFAIL, NERRS, NOUT )
*
                              CALL ZLACPY( 'Full', N, NRHS, B, LDB,
     $                                     WORK, LDB )
                              CALL ZGBT02( TRANS, M, N, KL, KU, NRHS, A,
     $                                     LDA, X, LDB, WORK, LDB,
     $                                     RESULT( 2 ) )
*
*+    TEST 3:
*                             Check solution from generated exact
*                             solution.
*
                              CALL ZGET04( N, NRHS, X, LDB, XACT, LDB,
     $                                     RCONDC, RESULT( 3 ) )
*
*+    TESTS 4, 5, 6:
*                             Use iterative refinement to improve the
*                             solution.
*
                              SRNAMT = 'ZGBRFS'
                              CALL ZGBRFS( TRANS, N, KL, KU, NRHS, A,
     $                                     LDA, AFAC, LDAFAC, IWORK, B,
     $                                     LDB, X, LDB, RWORK,
     $                                     RWORK( NRHS+1 ), WORK,
     $                                     RWORK( 2*NRHS+1 ), INFO )
*
*                             Check error code from ZGBRFS.
*
                              IF( INFO.NE.0 )
     $                           CALL ALAERH( PATH, 'ZGBRFS', INFO, 0,
     $                                        TRANS, N, N, KL, KU, NRHS,
     $                                        IMAT, NFAIL, NERRS, NOUT )
*
                              CALL ZGET04( N, NRHS, X, LDB, XACT, LDB,
     $                                     RCONDC, RESULT( 4 ) )
                              CALL ZGBT05( TRANS, N, KL, KU, NRHS, A,
     $                                     LDA, B, LDB, X, LDB, XACT,
     $                                     LDB, RWORK, RWORK( NRHS+1 ),
     $                                     RESULT( 5 ) )
*
*                             Print information about the tests that did
*                             not pass the threshold.
*
                              DO 60 K = 2, 6
                                 IF( RESULT( K ).GE.THRESH ) THEN
                                    IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                                 CALL ALAHD( NOUT, PATH )
                                    WRITE( NOUT, FMT = 9996 )TRANS, N,
     $                                 KL, KU, NRHS, IMAT, K,
     $                                 RESULT( K )
                                    NFAIL = NFAIL + 1
                                 END IF
   60                         CONTINUE
                              NRUN = NRUN + 5
   70                      CONTINUE
   80                   CONTINUE
*
*+    TEST 7:
*                          Get an estimate of RCOND = 1/CNDNUM.
*
   90                   CONTINUE
                        DO 100 ITRAN = 1, 2
                           IF( ITRAN.EQ.1 ) THEN
                              ANORM = ANORMO
                              RCONDC = RCONDO
                              NORM = 'O'
                           ELSE
                              ANORM = ANORMI
                              RCONDC = RCONDI
                              NORM = 'I'
                           END IF
                           SRNAMT = 'ZGBCON'
                           CALL ZGBCON( NORM, N, KL, KU, AFAC, LDAFAC,
     $                                  IWORK, ANORM, RCOND, WORK,
     $                                  RWORK, INFO )
*
*                             Check error code from ZGBCON.
*
                           IF( INFO.NE.0 )
     $                        CALL ALAERH( PATH, 'ZGBCON', INFO, 0,
     $                                     NORM, N, N, KL, KU, -1, IMAT,
     $                                     NFAIL, NERRS, NOUT )
*
                           RESULT( 7 ) = DGET06( RCOND, RCONDC )
*
*                          Print information about the tests that did
*                          not pass the threshold.
*
                           IF( RESULT( 7 ).GE.THRESH ) THEN
                              IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                           CALL ALAHD( NOUT, PATH )
                              WRITE( NOUT, FMT = 9995 )NORM, N, KL, KU,
     $                           IMAT, 7, RESULT( 7 )
                              NFAIL = NFAIL + 1
                           END IF
                           NRUN = NRUN + 1
  100                   CONTINUE
  110                CONTINUE
  120             CONTINUE
  130          CONTINUE
  140       CONTINUE
  150    CONTINUE
  160 CONTINUE
*
*     Print a summary of the results.
*
      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
 9999 FORMAT( ' *** In ZCHKGB, LA=', I5, ' is too small for M=', I5,
     $      ', N=', I5, ', KL=', I4, ', KU=', I4,
     $      / ' ==> Increase LA to at least ', I5 )
 9998 FORMAT( ' *** In ZCHKGB, LAFAC=', I5, ' is too small for M=', I5,
     $      ', N=', I5, ', KL=', I4, ', KU=', I4,
     $      / ' ==> Increase LAFAC to at least ', I5 )
 9997 FORMAT( ' M =', I5, ', N =', I5, ', KL=', I5, ', KU=', I5,
     $      ', NB =', I4, ', type ', I1, ', test(', I1, ')=', G12.5 )
 9996 FORMAT( ' TRANS=''', A1, ''', N=', I5, ', KL=', I5, ', KU=', I5,
     $      ', NRHS=', I3, ', type ', I1, ', test(', I1, ')=', G12.5 )
 9995 FORMAT( ' NORM =''', A1, ''', N=', I5, ', KL=', I5, ', KU=', I5,
     $      ',', 10X, ' type ', I1, ', test(', I1, ')=', G12.5 )
*
      RETURN
*
*     End of ZCHKGB
*
      END