LAPACK: TESTING/LIN/zchkhp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZCHKHP( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
00002      $                   NMAX, A, AFAC, AINV, B, X, XACT, WORK, RWORK,
00003      $                   IWORK, NOUT )
00004 *
00005 *  -- LAPACK test routine (version 3.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            TSTERR
00011       INTEGER            NMAX, NN, NNS, NOUT
00012       DOUBLE PRECISION   THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
00017       DOUBLE PRECISION   RWORK( * )
00018       COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
00019      $                   WORK( * ), X( * ), XACT( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  ZCHKHP tests ZHPTRF, -TRI, -TRS, -RFS, and -CON
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00031 *          The matrix types to be used for testing.  Matrices of type j
00032 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
00033 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
00034 *
00035 *  NN      (input) INTEGER
00036 *          The number of values of N contained in the vector NVAL.
00037 *
00038 *  NVAL    (input) INTEGER array, dimension (NN)
00039 *          The values of the matrix dimension N.
00040 *
00041 *  NNS     (input) INTEGER
00042 *          The number of values of NRHS contained in the vector NSVAL.
00043 *
00044 *  NSVAL   (input) INTEGER array, dimension (NNS)
00045 *          The values of the number of right hand sides NRHS.
00046 *
00047 *  THRESH  (input) DOUBLE PRECISION
00048 *          The threshold value for the test ratios.  A result is
00049 *          included in the output file if RESULT >= THRESH.  To have
00050 *          every test ratio printed, use THRESH = 0.
00051 *
00052 *  TSTERR  (input) LOGICAL
00053 *          Flag that indicates whether error exits are to be tested.
00054 *
00055 *  NMAX    (input) INTEGER
00056 *          The maximum value permitted for N, used in dimensioning the
00057 *          work arrays.
00058 *
00059 *  A       (workspace) COMPLEX*16 array, dimension
00060 *                      (NMAX*(NMAX+1)/2)
00061 *
00062 *  AFAC    (workspace) COMPLEX*16 array, dimension
00063 *                      (NMAX*(NMAX+1)/2)
00064 *
00065 *  AINV    (workspace) COMPLEX*16 array, dimension
00066 *                      (NMAX*(NMAX+1)/2)
00067 *
00068 *  B       (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX)
00069 *          where NSMAX is the largest entry in NSVAL.
00070 *
00071 *  X       (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX)
00072 *
00073 *  XACT    (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX)
00074 *
00075 *  WORK    (workspace) COMPLEX*16 array, dimension
00076 *                      (NMAX*max(2,NSMAX))
00077 *
00078 *  RWORK   (workspace) DOUBLE PRECISION array,
00079 *                                 dimension (NMAX+2*NSMAX)
00080 *
00081 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
00082 *
00083 *  NOUT    (input) INTEGER
00084 *          The unit number for output.
00085 *
00086 *  =====================================================================
00087 *
00088 *     .. Parameters ..
00089       DOUBLE PRECISION   ZERO
00090       PARAMETER          ( ZERO = 0.0D+0 )
00091       INTEGER            NTYPES
00092       PARAMETER          ( NTYPES = 10 )
00093       INTEGER            NTESTS
00094       PARAMETER          ( NTESTS = 8 )
00095 *     ..
00096 *     .. Local Scalars ..
00097       LOGICAL            TRFCON, ZEROT
00098       CHARACTER          DIST, PACKIT, TYPE, UPLO, XTYPE
00099       CHARACTER*3        PATH
00100       INTEGER            I, I1, I2, IMAT, IN, INFO, IOFF, IRHS, IUPLO,
00101      $                   IZERO, J, K, KL, KU, LDA, MODE, N, NERRS,
00102      $                   NFAIL, NIMAT, NPP, NRHS, NRUN, NT
00103       DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
00104 *     ..
00105 *     .. Local Arrays ..
00106       CHARACTER          UPLOS( 2 )
00107       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00108       DOUBLE PRECISION   RESULT( NTESTS )
00109 *     ..
00110 *     .. External Functions ..
00111       LOGICAL            LSAME
00112       DOUBLE PRECISION   DGET06, ZLANHP
00113       EXTERNAL           LSAME, DGET06, ZLANHP
00114 *     ..
00115 *     .. External Subroutines ..
00116       EXTERNAL           ALAERH, ALAHD, ALASUM, ZCOPY, ZERRSY, ZGET04,
00117      $                   ZHPCON, ZHPRFS, ZHPT01, ZHPTRF, ZHPTRI, ZHPTRS,
00118      $                   ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS, ZPPT02,
00119      $                   ZPPT03, ZPPT05
00120 *     ..
00121 *     .. Intrinsic Functions ..
00122       INTRINSIC          MAX, MIN
00123 *     ..
00124 *     .. Scalars in Common ..
00125       LOGICAL            LERR, OK
00126       CHARACTER*32       SRNAMT
00127       INTEGER            INFOT, NUNIT
00128 *     ..
00129 *     .. Common blocks ..
00130       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00131       COMMON             / SRNAMC / SRNAMT
00132 *     ..
00133 *     .. Data statements ..
00134       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00135       DATA               UPLOS / 'U', 'L' /
00136 *     ..
00137 *     .. Executable Statements ..
00138 *
00139 *     Initialize constants and the random number seed.
00140 *
00141       PATH( 1: 1 ) = 'Zomplex precision'
00142       PATH( 2: 3 ) = 'HP'
00143       NRUN = 0
00144       NFAIL = 0
00145       NERRS = 0
00146       DO 10 I = 1, 4
00147          ISEED( I ) = ISEEDY( I )
00148    10 CONTINUE
00149 *
00150 *     Test the error exits
00151 *
00152       IF( TSTERR )
00153      $   CALL ZERRSY( PATH, NOUT )
00154       INFOT = 0
00155 *
00156 *     Do for each value of N in NVAL
00157 *
00158       DO 170 IN = 1, NN
00159          N = NVAL( IN )
00160          LDA = MAX( N, 1 )
00161          XTYPE = 'N'
00162          NIMAT = NTYPES
00163          IF( N.LE.0 )
00164      $      NIMAT = 1
00165 *
00166          IZERO = 0
00167          DO 160 IMAT = 1, NIMAT
00168 *
00169 *           Do the tests only if DOTYPE( IMAT ) is true.
00170 *
00171             IF( .NOT.DOTYPE( IMAT ) )
00172      $         GO TO 160
00173 *
00174 *           Skip types 3, 4, 5, or 6 if the matrix size is too small.
00175 *
00176             ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
00177             IF( ZEROT .AND. N.LT.IMAT-2 )
00178      $         GO TO 160
00179 *
00180 *           Do first for UPLO = 'U', then for UPLO = 'L'
00181 *
00182             DO 150 IUPLO = 1, 2
00183                UPLO = UPLOS( IUPLO )
00184                IF( LSAME( UPLO, 'U' ) ) THEN
00185                   PACKIT = 'C'
00186                ELSE
00187                   PACKIT = 'R'
00188                END IF
00189 *
00190 *              Set up parameters with ZLATB4 and generate a test matrix
00191 *              with ZLATMS.
00192 *
00193                CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00194      $                      CNDNUM, DIST )
00195 *
00196                SRNAMT = 'ZLATMS'
00197                CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00198      $                      CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK,
00199      $                      INFO )
00200 *
00201 *              Check error code from ZLATMS.
00202 *
00203                IF( INFO.NE.0 ) THEN
00204                   CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1,
00205      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00206                   GO TO 150
00207                END IF
00208 *
00209 *              For types 3-6, zero one or more rows and columns of
00210 *              the matrix to test that INFO is returned correctly.
00211 *
00212                IF( ZEROT ) THEN
00213                   IF( IMAT.EQ.3 ) THEN
00214                      IZERO = 1
00215                   ELSE IF( IMAT.EQ.4 ) THEN
00216                      IZERO = N
00217                   ELSE
00218                      IZERO = N / 2 + 1
00219                   END IF
00220 *
00221                   IF( IMAT.LT.6 ) THEN
00222 *
00223 *                    Set row and column IZERO to zero.
00224 *
00225                      IF( IUPLO.EQ.1 ) THEN
00226                         IOFF = ( IZERO-1 )*IZERO / 2
00227                         DO 20 I = 1, IZERO - 1
00228                            A( IOFF+I ) = ZERO
00229    20                   CONTINUE
00230                         IOFF = IOFF + IZERO
00231                         DO 30 I = IZERO, N
00232                            A( IOFF ) = ZERO
00233                            IOFF = IOFF + I
00234    30                   CONTINUE
00235                      ELSE
00236                         IOFF = IZERO
00237                         DO 40 I = 1, IZERO - 1
00238                            A( IOFF ) = ZERO
00239                            IOFF = IOFF + N - I
00240    40                   CONTINUE
00241                         IOFF = IOFF - IZERO
00242                         DO 50 I = IZERO, N
00243                            A( IOFF+I ) = ZERO
00244    50                   CONTINUE
00245                      END IF
00246                   ELSE
00247                      IOFF = 0
00248                      IF( IUPLO.EQ.1 ) THEN
00249 *
00250 *                       Set the first IZERO rows and columns to zero.
00251 *
00252                         DO 70 J = 1, N
00253                            I2 = MIN( J, IZERO )
00254                            DO 60 I = 1, I2
00255                               A( IOFF+I ) = ZERO
00256    60                      CONTINUE
00257                            IOFF = IOFF + J
00258    70                   CONTINUE
00259                      ELSE
00260 *
00261 *                       Set the last IZERO rows and columns to zero.
00262 *
00263                         DO 90 J = 1, N
00264                            I1 = MAX( J, IZERO )
00265                            DO 80 I = I1, N
00266                               A( IOFF+I ) = ZERO
00267    80                      CONTINUE
00268                            IOFF = IOFF + N - J
00269    90                   CONTINUE
00270                      END IF
00271                   END IF
00272                ELSE
00273                   IZERO = 0
00274                END IF
00275 *
00276 *              Set the imaginary part of the diagonals.
00277 *
00278                IF( IUPLO.EQ.1 ) THEN
00279                   CALL ZLAIPD( N, A, 2, 1 )
00280                ELSE
00281                   CALL ZLAIPD( N, A, N, -1 )
00282                END IF
00283 *
00284 *              Compute the L*D*L' or U*D*U' factorization of the matrix.
00285 *
00286                NPP = N*( N+1 ) / 2
00287                CALL ZCOPY( NPP, A, 1, AFAC, 1 )
00288                SRNAMT = 'ZHPTRF'
00289                CALL ZHPTRF( UPLO, N, AFAC, IWORK, INFO )
00290 *
00291 *              Adjust the expected value of INFO to account for
00292 *              pivoting.
00293 *
00294                K = IZERO
00295                IF( K.GT.0 ) THEN
00296   100             CONTINUE
00297                   IF( IWORK( K ).LT.0 ) THEN
00298                      IF( IWORK( K ).NE.-K ) THEN
00299                         K = -IWORK( K )
00300                         GO TO 100
00301                      END IF
00302                   ELSE IF( IWORK( K ).NE.K ) THEN
00303                      K = IWORK( K )
00304                      GO TO 100
00305                   END IF
00306                END IF
00307 *
00308 *              Check error code from ZHPTRF.
00309 *
00310                IF( INFO.NE.K )
00311      $            CALL ALAERH( PATH, 'ZHPTRF', INFO, K, UPLO, N, N, -1,
00312      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00313                IF( INFO.NE.0 ) THEN
00314                   TRFCON = .TRUE.
00315                ELSE
00316                   TRFCON = .FALSE.
00317                END IF
00318 *
00319 *+    TEST 1
00320 *              Reconstruct matrix from factors and compute residual.
00321 *
00322                CALL ZHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK,
00323      $                      RESULT( 1 ) )
00324                NT = 1
00325 *
00326 *+    TEST 2
00327 *              Form the inverse and compute the residual.
00328 *
00329                IF( .NOT.TRFCON ) THEN
00330                   CALL ZCOPY( NPP, AFAC, 1, AINV, 1 )
00331                   SRNAMT = 'ZHPTRI'
00332                   CALL ZHPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
00333 *
00334 *              Check error code from ZHPTRI.
00335 *
00336                   IF( INFO.NE.0 )
00337      $               CALL ALAERH( PATH, 'ZHPTRI', INFO, 0, UPLO, N, N,
00338      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
00339 *
00340                   CALL ZPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK,
00341      $                         RCONDC, RESULT( 2 ) )
00342                   NT = 2
00343                END IF
00344 *
00345 *              Print information about the tests that did not pass
00346 *              the threshold.
00347 *
00348                DO 110 K = 1, NT
00349                   IF( RESULT( K ).GE.THRESH ) THEN
00350                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00351      $                  CALL ALAHD( NOUT, PATH )
00352                      WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, K,
00353      $                  RESULT( K )
00354                      NFAIL = NFAIL + 1
00355                   END IF
00356   110          CONTINUE
00357                NRUN = NRUN + NT
00358 *
00359 *              Do only the condition estimate if INFO is not 0.
00360 *
00361                IF( TRFCON ) THEN
00362                   RCONDC = ZERO
00363                   GO TO 140
00364                END IF
00365 *
00366                DO 130 IRHS = 1, NNS
00367                   NRHS = NSVAL( IRHS )
00368 *
00369 *+    TEST 3
00370 *              Solve and compute residual for  A * X = B.
00371 *
00372                   SRNAMT = 'ZLARHS'
00373                   CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
00374      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
00375      $                         INFO )
00376                   XTYPE = 'C'
00377                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00378 *
00379                   SRNAMT = 'ZHPTRS'
00380                   CALL ZHPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
00381      $                         INFO )
00382 *
00383 *              Check error code from ZHPTRS.
00384 *
00385                   IF( INFO.NE.0 )
00386      $               CALL ALAERH( PATH, 'ZHPTRS', INFO, 0, UPLO, N, N,
00387      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00388      $                            NOUT )
00389 *
00390                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00391                   CALL ZPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
00392      $                         RWORK, RESULT( 3 ) )
00393 *
00394 *+    TEST 4
00395 *              Check solution from generated exact solution.
00396 *
00397                   CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00398      $                         RESULT( 4 ) )
00399 *
00400 *+    TESTS 5, 6, and 7
00401 *              Use iterative refinement to improve the solution.
00402 *
00403                   SRNAMT = 'ZHPRFS'
00404                   CALL ZHPRFS( UPLO, N, NRHS, A, AFAC, IWORK, B, LDA, X,
00405      $                         LDA, RWORK, RWORK( NRHS+1 ), WORK,
00406      $                         RWORK( 2*NRHS+1 ), INFO )
00407 *
00408 *              Check error code from ZHPRFS.
00409 *
00410                   IF( INFO.NE.0 )
00411      $               CALL ALAERH( PATH, 'ZHPRFS', INFO, 0, UPLO, N, N,
00412      $                            -1, -1, NRHS, IMAT, NFAIL, NERRS,
00413      $                            NOUT )
00414 *
00415                   CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00416      $                         RESULT( 5 ) )
00417                   CALL ZPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA, XACT,
00418      $                         LDA, RWORK, RWORK( NRHS+1 ),
00419      $                         RESULT( 6 ) )
00420 *
00421 *                 Print information about the tests that did not pass
00422 *                 the threshold.
00423 *
00424                   DO 120 K = 3, 7
00425                      IF( RESULT( K ).GE.THRESH ) THEN
00426                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00427      $                     CALL ALAHD( NOUT, PATH )
00428                         WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT,
00429      $                     K, RESULT( K )
00430                         NFAIL = NFAIL + 1
00431                      END IF
00432   120             CONTINUE
00433                   NRUN = NRUN + 5
00434   130          CONTINUE
00435 *
00436 *+    TEST 8
00437 *              Get an estimate of RCOND = 1/CNDNUM.
00438 *
00439   140          CONTINUE
00440                ANORM = ZLANHP( '1', UPLO, N, A, RWORK )
00441                SRNAMT = 'ZHPCON'
00442                CALL ZHPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK,
00443      $                      INFO )
00444 *
00445 *              Check error code from ZHPCON.
00446 *
00447                IF( INFO.NE.0 )
00448      $            CALL ALAERH( PATH, 'ZHPCON', INFO, 0, UPLO, N, N, -1,
00449      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00450 *
00451                RESULT( 8 ) = DGET06( RCOND, RCONDC )
00452 *
00453 *              Print the test ratio if it is .GE. THRESH.
00454 *
00455                IF( RESULT( 8 ).GE.THRESH ) THEN
00456                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00457      $               CALL ALAHD( NOUT, PATH )
00458                   WRITE( NOUT, FMT = 9999 )UPLO, N, IMAT, 8,
00459      $               RESULT( 8 )
00460                   NFAIL = NFAIL + 1
00461                END IF
00462                NRUN = NRUN + 1
00463   150       CONTINUE
00464   160    CONTINUE
00465   170 CONTINUE
00466 *
00467 *     Print a summary of the results.
00468 *
00469       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
00470 *
00471  9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', type ', I2, ', test ',
00472      $      I2, ', ratio =', G12.5 )
00473  9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
00474      $      I2, ', test(', I2, ') =', G12.5 )
00475       RETURN
00476 *
00477 *     End of ZCHKHP
00478 *
00479       END