LAPACK: TESTING/LIN/zdrvgt.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZDRVGT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, AF,
00002      $                   B, X, XACT, WORK, RWORK, IWORK, NOUT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            TSTERR
00010       INTEGER            NN, NOUT, NRHS
00011       DOUBLE PRECISION   THRESH
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            DOTYPE( * )
00015       INTEGER            IWORK( * ), NVAL( * )
00016       DOUBLE PRECISION   RWORK( * )
00017       COMPLEX*16         A( * ), AF( * ), B( * ), WORK( * ), X( * ),
00018      $                   XACT( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  ZDRVGT tests ZGTSV and -SVX.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00030 *          The matrix types to be used for testing.  Matrices of type j
00031 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
00032 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
00033 *
00034 *  NN      (input) INTEGER
00035 *          The number of values of N contained in the vector NVAL.
00036 *
00037 *  NVAL    (input) INTEGER array, dimension (NN)
00038 *          The values of the matrix dimension N.
00039 *
00040 *  THRESH  (input) DOUBLE PRECISION
00041 *          The threshold value for the test ratios.  A result is
00042 *          included in the output file if RESULT >= THRESH.  To have
00043 *          every test ratio printed, use THRESH = 0.
00044 *
00045 *  TSTERR  (input) LOGICAL
00046 *          Flag that indicates whether error exits are to be tested.
00047 *
00048 *  A       (workspace) COMPLEX*16 array, dimension (NMAX*4)
00049 *
00050 *  AF      (workspace) COMPLEX*16 array, dimension (NMAX*4)
00051 *
00052 *  B       (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
00053 *
00054 *  X       (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
00055 *
00056 *  XACT    (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
00057 *
00058 *  WORK    (workspace) COMPLEX*16 array, dimension
00059 *                      (NMAX*max(3,NRHS))
00060 *
00061 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
00062 *
00063 *  IWORK   (workspace) INTEGER array, dimension (2*NMAX)
00064 *
00065 *  NOUT    (input) INTEGER
00066 *          The unit number for output.
00067 *
00068 *  =====================================================================
00069 *
00070 *     .. Parameters ..
00071       DOUBLE PRECISION   ONE, ZERO
00072       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00073       INTEGER            NTYPES
00074       PARAMETER          ( NTYPES = 12 )
00075       INTEGER            NTESTS
00076       PARAMETER          ( NTESTS = 6 )
00077 *     ..
00078 *     .. Local Scalars ..
00079       LOGICAL            TRFCON, ZEROT
00080       CHARACTER          DIST, FACT, TRANS, TYPE
00081       CHARACTER*3        PATH
00082       INTEGER            I, IFACT, IMAT, IN, INFO, ITRAN, IX, IZERO, J,
00083      $                   K, K1, KL, KOFF, KU, LDA, M, MODE, N, NERRS,
00084      $                   NFAIL, NIMAT, NRUN, NT
00085       DOUBLE PRECISION   AINVNM, ANORM, ANORMI, ANORMO, COND, RCOND,
00086      $                   RCONDC, RCONDI, RCONDO
00087 *     ..
00088 *     .. Local Arrays ..
00089       CHARACTER          TRANSS( 3 )
00090       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00091       DOUBLE PRECISION   RESULT( NTESTS ), Z( 3 )
00092 *     ..
00093 *     .. External Functions ..
00094       DOUBLE PRECISION   DGET06, DZASUM, ZLANGT
00095       EXTERNAL           DGET06, DZASUM, ZLANGT
00096 *     ..
00097 *     .. External Subroutines ..
00098       EXTERNAL           ALADHD, ALAERH, ALASVM, ZCOPY, ZDSCAL, ZERRVX,
00099      $                   ZGET04, ZGTSV, ZGTSVX, ZGTT01, ZGTT02, ZGTT05,
00100      $                   ZGTTRF, ZGTTRS, ZLACPY, ZLAGTM, ZLARNV, ZLASET,
00101      $                   ZLATB4, ZLATMS
00102 *     ..
00103 *     .. Intrinsic Functions ..
00104       INTRINSIC          DCMPLX, MAX
00105 *     ..
00106 *     .. Scalars in Common ..
00107       LOGICAL            LERR, OK
00108       CHARACTER*32       SRNAMT
00109       INTEGER            INFOT, NUNIT
00110 *     ..
00111 *     .. Common blocks ..
00112       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00113       COMMON             / SRNAMC / SRNAMT
00114 *     ..
00115 *     .. Data statements ..
00116       DATA               ISEEDY / 0, 0, 0, 1 / , TRANSS / 'N', 'T',
00117      $                   'C' /
00118 *     ..
00119 *     .. Executable Statements ..
00120 *
00121       PATH( 1: 1 ) = 'Zomplex precision'
00122       PATH( 2: 3 ) = 'GT'
00123       NRUN = 0
00124       NFAIL = 0
00125       NERRS = 0
00126       DO 10 I = 1, 4
00127          ISEED( I ) = ISEEDY( I )
00128    10 CONTINUE
00129 *
00130 *     Test the error exits
00131 *
00132       IF( TSTERR )
00133      $   CALL ZERRVX( PATH, NOUT )
00134       INFOT = 0
00135 *
00136       DO 140 IN = 1, NN
00137 *
00138 *        Do for each value of N in NVAL.
00139 *
00140          N = NVAL( IN )
00141          M = MAX( N-1, 0 )
00142          LDA = MAX( 1, N )
00143          NIMAT = NTYPES
00144          IF( N.LE.0 )
00145      $      NIMAT = 1
00146 *
00147          DO 130 IMAT = 1, NIMAT
00148 *
00149 *           Do the tests only if DOTYPE( IMAT ) is true.
00150 *
00151             IF( .NOT.DOTYPE( IMAT ) )
00152      $         GO TO 130
00153 *
00154 *           Set up parameters with ZLATB4.
00155 *
00156             CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00157      $                   COND, DIST )
00158 *
00159             ZEROT = IMAT.GE.8 .AND. IMAT.LE.10
00160             IF( IMAT.LE.6 ) THEN
00161 *
00162 *              Types 1-6:  generate matrices of known condition number.
00163 *
00164                KOFF = MAX( 2-KU, 3-MAX( 1, N ) )
00165                SRNAMT = 'ZLATMS'
00166                CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND,
00167      $                      ANORM, KL, KU, 'Z', AF( KOFF ), 3, WORK,
00168      $                      INFO )
00169 *
00170 *              Check the error code from ZLATMS.
00171 *
00172                IF( INFO.NE.0 ) THEN
00173                   CALL ALAERH( PATH, 'ZLATMS', INFO, 0, ' ', N, N, KL,
00174      $                         KU, -1, IMAT, NFAIL, NERRS, NOUT )
00175                   GO TO 130
00176                END IF
00177                IZERO = 0
00178 *
00179                IF( N.GT.1 ) THEN
00180                   CALL ZCOPY( N-1, AF( 4 ), 3, A, 1 )
00181                   CALL ZCOPY( N-1, AF( 3 ), 3, A( N+M+1 ), 1 )
00182                END IF
00183                CALL ZCOPY( N, AF( 2 ), 3, A( M+1 ), 1 )
00184             ELSE
00185 *
00186 *              Types 7-12:  generate tridiagonal matrices with
00187 *              unknown condition numbers.
00188 *
00189                IF( .NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN
00190 *
00191 *                 Generate a matrix with elements from [-1,1].
00192 *
00193                   CALL ZLARNV( 2, ISEED, N+2*M, A )
00194                   IF( ANORM.NE.ONE )
00195      $               CALL ZDSCAL( N+2*M, ANORM, A, 1 )
00196                ELSE IF( IZERO.GT.0 ) THEN
00197 *
00198 *                 Reuse the last matrix by copying back the zeroed out
00199 *                 elements.
00200 *
00201                   IF( IZERO.EQ.1 ) THEN
00202                      A( N ) = Z( 2 )
00203                      IF( N.GT.1 )
00204      $                  A( 1 ) = Z( 3 )
00205                   ELSE IF( IZERO.EQ.N ) THEN
00206                      A( 3*N-2 ) = Z( 1 )
00207                      A( 2*N-1 ) = Z( 2 )
00208                   ELSE
00209                      A( 2*N-2+IZERO ) = Z( 1 )
00210                      A( N-1+IZERO ) = Z( 2 )
00211                      A( IZERO ) = Z( 3 )
00212                   END IF
00213                END IF
00214 *
00215 *              If IMAT > 7, set one column of the matrix to 0.
00216 *
00217                IF( .NOT.ZEROT ) THEN
00218                   IZERO = 0
00219                ELSE IF( IMAT.EQ.8 ) THEN
00220                   IZERO = 1
00221                   Z( 2 ) = A( N )
00222                   A( N ) = ZERO
00223                   IF( N.GT.1 ) THEN
00224                      Z( 3 ) = A( 1 )
00225                      A( 1 ) = ZERO
00226                   END IF
00227                ELSE IF( IMAT.EQ.9 ) THEN
00228                   IZERO = N
00229                   Z( 1 ) = A( 3*N-2 )
00230                   Z( 2 ) = A( 2*N-1 )
00231                   A( 3*N-2 ) = ZERO
00232                   A( 2*N-1 ) = ZERO
00233                ELSE
00234                   IZERO = ( N+1 ) / 2
00235                   DO 20 I = IZERO, N - 1
00236                      A( 2*N-2+I ) = ZERO
00237                      A( N-1+I ) = ZERO
00238                      A( I ) = ZERO
00239    20             CONTINUE
00240                   A( 3*N-2 ) = ZERO
00241                   A( 2*N-1 ) = ZERO
00242                END IF
00243             END IF
00244 *
00245             DO 120 IFACT = 1, 2
00246                IF( IFACT.EQ.1 ) THEN
00247                   FACT = 'F'
00248                ELSE
00249                   FACT = 'N'
00250                END IF
00251 *
00252 *              Compute the condition number for comparison with
00253 *              the value returned by ZGTSVX.
00254 *
00255                IF( ZEROT ) THEN
00256                   IF( IFACT.EQ.1 )
00257      $               GO TO 120
00258                   RCONDO = ZERO
00259                   RCONDI = ZERO
00260 *
00261                ELSE IF( IFACT.EQ.1 ) THEN
00262                   CALL ZCOPY( N+2*M, A, 1, AF, 1 )
00263 *
00264 *                 Compute the 1-norm and infinity-norm of A.
00265 *
00266                   ANORMO = ZLANGT( '1', N, A, A( M+1 ), A( N+M+1 ) )
00267                   ANORMI = ZLANGT( 'I', N, A, A( M+1 ), A( N+M+1 ) )
00268 *
00269 *                 Factor the matrix A.
00270 *
00271                   CALL ZGTTRF( N, AF, AF( M+1 ), AF( N+M+1 ),
00272      $                         AF( N+2*M+1 ), IWORK, INFO )
00273 *
00274 *                 Use ZGTTRS to solve for one column at a time of
00275 *                 inv(A), computing the maximum column sum as we go.
00276 *
00277                   AINVNM = ZERO
00278                   DO 40 I = 1, N
00279                      DO 30 J = 1, N
00280                         X( J ) = ZERO
00281    30                CONTINUE
00282                      X( I ) = ONE
00283                      CALL ZGTTRS( 'No transpose', N, 1, AF, AF( M+1 ),
00284      $                            AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
00285      $                            LDA, INFO )
00286                      AINVNM = MAX( AINVNM, DZASUM( N, X, 1 ) )
00287    40             CONTINUE
00288 *
00289 *                 Compute the 1-norm condition number of A.
00290 *
00291                   IF( ANORMO.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00292                      RCONDO = ONE
00293                   ELSE
00294                      RCONDO = ( ONE / ANORMO ) / AINVNM
00295                   END IF
00296 *
00297 *                 Use ZGTTRS to solve for one column at a time of
00298 *                 inv(A'), computing the maximum column sum as we go.
00299 *
00300                   AINVNM = ZERO
00301                   DO 60 I = 1, N
00302                      DO 50 J = 1, N
00303                         X( J ) = ZERO
00304    50                CONTINUE
00305                      X( I ) = ONE
00306                      CALL ZGTTRS( 'Conjugate transpose', N, 1, AF,
00307      $                            AF( M+1 ), AF( N+M+1 ), AF( N+2*M+1 ),
00308      $                            IWORK, X, LDA, INFO )
00309                      AINVNM = MAX( AINVNM, DZASUM( N, X, 1 ) )
00310    60             CONTINUE
00311 *
00312 *                 Compute the infinity-norm condition number of A.
00313 *
00314                   IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00315                      RCONDI = ONE
00316                   ELSE
00317                      RCONDI = ( ONE / ANORMI ) / AINVNM
00318                   END IF
00319                END IF
00320 *
00321                DO 110 ITRAN = 1, 3
00322                   TRANS = TRANSS( ITRAN )
00323                   IF( ITRAN.EQ.1 ) THEN
00324                      RCONDC = RCONDO
00325                   ELSE
00326                      RCONDC = RCONDI
00327                   END IF
00328 *
00329 *                 Generate NRHS random solution vectors.
00330 *
00331                   IX = 1
00332                   DO 70 J = 1, NRHS
00333                      CALL ZLARNV( 2, ISEED, N, XACT( IX ) )
00334                      IX = IX + LDA
00335    70             CONTINUE
00336 *
00337 *                 Set the right hand side.
00338 *
00339                   CALL ZLAGTM( TRANS, N, NRHS, ONE, A, A( M+1 ),
00340      $                         A( N+M+1 ), XACT, LDA, ZERO, B, LDA )
00341 *
00342                   IF( IFACT.EQ.2 .AND. ITRAN.EQ.1 ) THEN
00343 *
00344 *                    --- Test ZGTSV  ---
00345 *
00346 *                    Solve the system using Gaussian elimination with
00347 *                    partial pivoting.
00348 *
00349                      CALL ZCOPY( N+2*M, A, 1, AF, 1 )
00350                      CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00351 *
00352                      SRNAMT = 'ZGTSV '
00353                      CALL ZGTSV( N, NRHS, AF, AF( M+1 ), AF( N+M+1 ), X,
00354      $                           LDA, INFO )
00355 *
00356 *                    Check error code from ZGTSV .
00357 *
00358                      IF( INFO.NE.IZERO )
00359      $                  CALL ALAERH( PATH, 'ZGTSV ', INFO, IZERO, ' ',
00360      $                               N, N, 1, 1, NRHS, IMAT, NFAIL,
00361      $                               NERRS, NOUT )
00362                      NT = 1
00363                      IF( IZERO.EQ.0 ) THEN
00364 *
00365 *                       Check residual of computed solution.
00366 *
00367                         CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK,
00368      $                               LDA )
00369                         CALL ZGTT02( TRANS, N, NRHS, A, A( M+1 ),
00370      $                               A( N+M+1 ), X, LDA, WORK, LDA,
00371      $                               RWORK, RESULT( 2 ) )
00372 *
00373 *                       Check solution from generated exact solution.
00374 *
00375                         CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00376      $                               RESULT( 3 ) )
00377                         NT = 3
00378                      END IF
00379 *
00380 *                    Print information about the tests that did not pass
00381 *                    the threshold.
00382 *
00383                      DO 80 K = 2, NT
00384                         IF( RESULT( K ).GE.THRESH ) THEN
00385                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00386      $                        CALL ALADHD( NOUT, PATH )
00387                            WRITE( NOUT, FMT = 9999 )'ZGTSV ', N, IMAT,
00388      $                        K, RESULT( K )
00389                            NFAIL = NFAIL + 1
00390                         END IF
00391    80                CONTINUE
00392                      NRUN = NRUN + NT - 1
00393                   END IF
00394 *
00395 *                 --- Test ZGTSVX ---
00396 *
00397                   IF( IFACT.GT.1 ) THEN
00398 *
00399 *                    Initialize AF to zero.
00400 *
00401                      DO 90 I = 1, 3*N - 2
00402                         AF( I ) = ZERO
00403    90                CONTINUE
00404                   END IF
00405                   CALL ZLASET( 'Full', N, NRHS, DCMPLX( ZERO ),
00406      $                         DCMPLX( ZERO ), X, LDA )
00407 *
00408 *                 Solve the system and compute the condition number and
00409 *                 error bounds using ZGTSVX.
00410 *
00411                   SRNAMT = 'ZGTSVX'
00412                   CALL ZGTSVX( FACT, TRANS, N, NRHS, A, A( M+1 ),
00413      $                         A( N+M+1 ), AF, AF( M+1 ), AF( N+M+1 ),
00414      $                         AF( N+2*M+1 ), IWORK, B, LDA, X, LDA,
00415      $                         RCOND, RWORK, RWORK( NRHS+1 ), WORK,
00416      $                         RWORK( 2*NRHS+1 ), INFO )
00417 *
00418 *                 Check the error code from ZGTSVX.
00419 *
00420                   IF( INFO.NE.IZERO )
00421      $               CALL ALAERH( PATH, 'ZGTSVX', INFO, IZERO,
00422      $                            FACT // TRANS, N, N, 1, 1, NRHS, IMAT,
00423      $                            NFAIL, NERRS, NOUT )
00424 *
00425                   IF( IFACT.GE.2 ) THEN
00426 *
00427 *                    Reconstruct matrix from factors and compute
00428 *                    residual.
00429 *
00430                      CALL ZGTT01( N, A, A( M+1 ), A( N+M+1 ), AF,
00431      $                            AF( M+1 ), AF( N+M+1 ), AF( N+2*M+1 ),
00432      $                            IWORK, WORK, LDA, RWORK, RESULT( 1 ) )
00433                      K1 = 1
00434                   ELSE
00435                      K1 = 2
00436                   END IF
00437 *
00438                   IF( INFO.EQ.0 ) THEN
00439                      TRFCON = .FALSE.
00440 *
00441 *                    Check residual of computed solution.
00442 *
00443                      CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00444                      CALL ZGTT02( TRANS, N, NRHS, A, A( M+1 ),
00445      $                            A( N+M+1 ), X, LDA, WORK, LDA, RWORK,
00446      $                            RESULT( 2 ) )
00447 *
00448 *                    Check solution from generated exact solution.
00449 *
00450                      CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00451      $                            RESULT( 3 ) )
00452 *
00453 *                    Check the error bounds from iterative refinement.
00454 *
00455                      CALL ZGTT05( TRANS, N, NRHS, A, A( M+1 ),
00456      $                            A( N+M+1 ), B, LDA, X, LDA, XACT, LDA,
00457      $                            RWORK, RWORK( NRHS+1 ), RESULT( 4 ) )
00458                      NT = 5
00459                   END IF
00460 *
00461 *                 Print information about the tests that did not pass
00462 *                 the threshold.
00463 *
00464                   DO 100 K = K1, NT
00465                      IF( RESULT( K ).GE.THRESH ) THEN
00466                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00467      $                     CALL ALADHD( NOUT, PATH )
00468                         WRITE( NOUT, FMT = 9998 )'ZGTSVX', FACT, TRANS,
00469      $                     N, IMAT, K, RESULT( K )
00470                         NFAIL = NFAIL + 1
00471                      END IF
00472   100             CONTINUE
00473 *
00474 *                 Check the reciprocal of the condition number.
00475 *
00476                   RESULT( 6 ) = DGET06( RCOND, RCONDC )
00477                   IF( RESULT( 6 ).GE.THRESH ) THEN
00478                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00479      $                  CALL ALADHD( NOUT, PATH )
00480                      WRITE( NOUT, FMT = 9998 )'ZGTSVX', FACT, TRANS, N,
00481      $                  IMAT, K, RESULT( K )
00482                      NFAIL = NFAIL + 1
00483                   END IF
00484                   NRUN = NRUN + NT - K1 + 2
00485 *
00486   110          CONTINUE
00487   120       CONTINUE
00488   130    CONTINUE
00489   140 CONTINUE
00490 *
00491 *     Print a summary of the results.
00492 *
00493       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
00494 *
00495  9999 FORMAT( 1X, A, ', N =', I5, ', type ', I2, ', test ', I2,
00496      $      ', ratio = ', G12.5 )
00497  9998 FORMAT( 1X, A, ', FACT=''', A1, ''', TRANS=''', A1, ''', N =',
00498      $      I5, ', type ', I2, ', test ', I2, ', ratio = ', G12.5 )
00499       RETURN
00500 *
00501 *     End of ZDRVGT
00502 *
00503       END