LAPACK: TESTING/LIN/cdrvhp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CDRVHP( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
00002      $                   A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
00003      $                   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, NOUT, NRHS
00012       REAL               THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            DOTYPE( * )
00016       INTEGER            IWORK( * ), NVAL( * )
00017       REAL               RWORK( * )
00018       COMPLEX            A( * ), AFAC( * ), AINV( * ), B( * ),
00019      $                   WORK( * ), X( * ), XACT( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  CDRVHP tests the driver routines CHPSV and -SVX.
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 *  NRHS    (input) INTEGER
00042 *          The number of right hand side vectors to be generated for
00043 *          each linear system.
00044 *
00045 *  THRESH  (input) REAL
00046 *          The threshold value for the test ratios.  A result is
00047 *          included in the output file if RESULT >= THRESH.  To have
00048 *          every test ratio printed, use THRESH = 0.
00049 *
00050 *  TSTERR  (input) LOGICAL
00051 *          Flag that indicates whether error exits are to be tested.
00052 *
00053 *  NMAX    (input) INTEGER
00054 *          The maximum value permitted for N, used in dimensioning the
00055 *          work arrays.
00056 *
00057 *  A       (workspace) COMPLEX array, dimension
00058 *                      (NMAX*(NMAX+1)/2)
00059 *
00060 *  AFAC    (workspace) COMPLEX array, dimension
00061 *                      (NMAX*(NMAX+1)/2)
00062 *
00063 *  AINV    (workspace) COMPLEX array, dimension
00064 *                      (NMAX*(NMAX+1)/2)
00065 *
00066 *  B       (workspace) COMPLEX array, dimension (NMAX*NRHS)
00067 *
00068 *  X       (workspace) COMPLEX array, dimension (NMAX*NRHS)
00069 *
00070 *  XACT    (workspace) COMPLEX array, dimension (NMAX*NRHS)
00071 *
00072 *  WORK    (workspace) COMPLEX array, dimension
00073 *                      (NMAX*max(2,NRHS))
00074 *
00075 *  RWORK   (workspace) REAL array, dimension (NMAX+2*NRHS)
00076 *
00077 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
00078 *
00079 *  NOUT    (input) INTEGER
00080 *          The unit number for output.
00081 *
00082 *  =====================================================================
00083 *
00084 *     .. Parameters ..
00085       REAL               ONE, ZERO
00086       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00087       INTEGER            NTYPES, NTESTS
00088       PARAMETER          ( NTYPES = 10, NTESTS = 6 )
00089       INTEGER            NFACT
00090       PARAMETER          ( NFACT = 2 )
00091 *     ..
00092 *     .. Local Scalars ..
00093       LOGICAL            ZEROT
00094       CHARACTER          DIST, FACT, PACKIT, TYPE, UPLO, XTYPE
00095       CHARACTER*3        PATH
00096       INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
00097      $                   IZERO, J, K, K1, KL, KU, LDA, MODE, N, NB,
00098      $                   NBMIN, NERRS, NFAIL, NIMAT, NPP, NRUN, NT
00099       REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
00100 *     ..
00101 *     .. Local Arrays ..
00102       CHARACTER          FACTS( NFACT )
00103       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00104       REAL               RESULT( NTESTS )
00105 *     ..
00106 *     .. External Functions ..
00107       REAL               CLANHP, SGET06
00108       EXTERNAL           CLANHP, SGET06
00109 *     ..
00110 *     .. External Subroutines ..
00111       EXTERNAL           ALADHD, ALAERH, ALASVM, CCOPY, CERRVX, CGET04,
00112      $                   CHPSV, CHPSVX, CHPT01, CHPTRF, CHPTRI, CLACPY,
00113      $                   CLAIPD, CLARHS, CLASET, CLATB4, CLATMS, CPPT02,
00114      $                   CPPT05, XLAENV
00115 *     ..
00116 *     .. Scalars in Common ..
00117       LOGICAL            LERR, OK
00118       CHARACTER*32       SRNAMT
00119       INTEGER            INFOT, NUNIT
00120 *     ..
00121 *     .. Common blocks ..
00122       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
00123       COMMON             / SRNAMC / SRNAMT
00124 *     ..
00125 *     .. Intrinsic Functions ..
00126       INTRINSIC          CMPLX, MAX, MIN
00127 *     ..
00128 *     .. Data statements ..
00129       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00130       DATA               FACTS / 'F', 'N' /
00131 *     ..
00132 *     .. Executable Statements ..
00133 *
00134 *     Initialize constants and the random number seed.
00135 *
00136       PATH( 1: 1 ) = 'C'
00137       PATH( 2: 3 ) = 'HP'
00138       NRUN = 0
00139       NFAIL = 0
00140       NERRS = 0
00141       DO 10 I = 1, 4
00142          ISEED( I ) = ISEEDY( I )
00143    10 CONTINUE
00144 *
00145 *     Test the error exits
00146 *
00147       IF( TSTERR )
00148      $   CALL CERRVX( PATH, NOUT )
00149       INFOT = 0
00150 *
00151 *     Set the block size and minimum block size for testing.
00152 *
00153       NB = 1
00154       NBMIN = 2
00155       CALL XLAENV( 1, NB )
00156       CALL XLAENV( 2, NBMIN )
00157 *
00158 *     Do for each value of N in NVAL
00159 *
00160       DO 180 IN = 1, NN
00161          N = NVAL( IN )
00162          LDA = MAX( N, 1 )
00163          NPP = N*( N+1 ) / 2
00164          XTYPE = 'N'
00165          NIMAT = NTYPES
00166          IF( N.LE.0 )
00167      $      NIMAT = 1
00168 *
00169          DO 170 IMAT = 1, NIMAT
00170 *
00171 *           Do the tests only if DOTYPE( IMAT ) is true.
00172 *
00173             IF( .NOT.DOTYPE( IMAT ) )
00174      $         GO TO 170
00175 *
00176 *           Skip types 3, 4, 5, or 6 if the matrix size is too small.
00177 *
00178             ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
00179             IF( ZEROT .AND. N.LT.IMAT-2 )
00180      $         GO TO 170
00181 *
00182 *           Do first for UPLO = 'U', then for UPLO = 'L'
00183 *
00184             DO 160 IUPLO = 1, 2
00185                IF( IUPLO.EQ.1 ) THEN
00186                   UPLO = 'U'
00187                   PACKIT = 'C'
00188                ELSE
00189                   UPLO = 'L'
00190                   PACKIT = 'R'
00191                END IF
00192 *
00193 *              Set up parameters with CLATB4 and generate a test matrix
00194 *              with CLATMS.
00195 *
00196                CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00197      $                      CNDNUM, DIST )
00198 *
00199                SRNAMT = 'CLATMS'
00200                CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00201      $                      CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK,
00202      $                      INFO )
00203 *
00204 *              Check error code from CLATMS.
00205 *
00206                IF( INFO.NE.0 ) THEN
00207                   CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1,
00208      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00209                   GO TO 160
00210                END IF
00211 *
00212 *              For types 3-6, zero one or more rows and columns of the
00213 *              matrix to test that INFO is returned correctly.
00214 *
00215                IF( ZEROT ) THEN
00216                   IF( IMAT.EQ.3 ) THEN
00217                      IZERO = 1
00218                   ELSE IF( IMAT.EQ.4 ) THEN
00219                      IZERO = N
00220                   ELSE
00221                      IZERO = N / 2 + 1
00222                   END IF
00223 *
00224                   IF( IMAT.LT.6 ) THEN
00225 *
00226 *                    Set row and column IZERO to zero.
00227 *
00228                      IF( IUPLO.EQ.1 ) THEN
00229                         IOFF = ( IZERO-1 )*IZERO / 2
00230                         DO 20 I = 1, IZERO - 1
00231                            A( IOFF+I ) = ZERO
00232    20                   CONTINUE
00233                         IOFF = IOFF + IZERO
00234                         DO 30 I = IZERO, N
00235                            A( IOFF ) = ZERO
00236                            IOFF = IOFF + I
00237    30                   CONTINUE
00238                      ELSE
00239                         IOFF = IZERO
00240                         DO 40 I = 1, IZERO - 1
00241                            A( IOFF ) = ZERO
00242                            IOFF = IOFF + N - I
00243    40                   CONTINUE
00244                         IOFF = IOFF - IZERO
00245                         DO 50 I = IZERO, N
00246                            A( IOFF+I ) = ZERO
00247    50                   CONTINUE
00248                      END IF
00249                   ELSE
00250                      IOFF = 0
00251                      IF( IUPLO.EQ.1 ) THEN
00252 *
00253 *                       Set the first IZERO rows and columns to zero.
00254 *
00255                         DO 70 J = 1, N
00256                            I2 = MIN( J, IZERO )
00257                            DO 60 I = 1, I2
00258                               A( IOFF+I ) = ZERO
00259    60                      CONTINUE
00260                            IOFF = IOFF + J
00261    70                   CONTINUE
00262                      ELSE
00263 *
00264 *                       Set the last IZERO rows and columns to zero.
00265 *
00266                         DO 90 J = 1, N
00267                            I1 = MAX( J, IZERO )
00268                            DO 80 I = I1, N
00269                               A( IOFF+I ) = ZERO
00270    80                      CONTINUE
00271                            IOFF = IOFF + N - J
00272    90                   CONTINUE
00273                      END IF
00274                   END IF
00275                ELSE
00276                   IZERO = 0
00277                END IF
00278 *
00279 *              Set the imaginary part of the diagonals.
00280 *
00281                IF( IUPLO.EQ.1 ) THEN
00282                   CALL CLAIPD( N, A, 2, 1 )
00283                ELSE
00284                   CALL CLAIPD( N, A, N, -1 )
00285                END IF
00286 *
00287                DO 150 IFACT = 1, NFACT
00288 *
00289 *                 Do first for FACT = 'F', then for other values.
00290 *
00291                   FACT = FACTS( IFACT )
00292 *
00293 *                 Compute the condition number for comparison with
00294 *                 the value returned by CHPSVX.
00295 *
00296                   IF( ZEROT ) THEN
00297                      IF( IFACT.EQ.1 )
00298      $                  GO TO 150
00299                      RCONDC = ZERO
00300 *
00301                   ELSE IF( IFACT.EQ.1 ) THEN
00302 *
00303 *                    Compute the 1-norm of A.
00304 *
00305                      ANORM = CLANHP( '1', UPLO, N, A, RWORK )
00306 *
00307 *                    Factor the matrix A.
00308 *
00309                      CALL CCOPY( NPP, A, 1, AFAC, 1 )
00310                      CALL CHPTRF( UPLO, N, AFAC, IWORK, INFO )
00311 *
00312 *                    Compute inv(A) and take its norm.
00313 *
00314                      CALL CCOPY( NPP, AFAC, 1, AINV, 1 )
00315                      CALL CHPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
00316                      AINVNM = CLANHP( '1', UPLO, N, AINV, RWORK )
00317 *
00318 *                    Compute the 1-norm condition number of A.
00319 *
00320                      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00321                         RCONDC = ONE
00322                      ELSE
00323                         RCONDC = ( ONE / ANORM ) / AINVNM
00324                      END IF
00325                   END IF
00326 *
00327 *                 Form an exact solution and set the right hand side.
00328 *
00329                   SRNAMT = 'CLARHS'
00330                   CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
00331      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
00332      $                         INFO )
00333                   XTYPE = 'C'
00334 *
00335 *                 --- Test CHPSV  ---
00336 *
00337                   IF( IFACT.EQ.2 ) THEN
00338                      CALL CCOPY( NPP, A, 1, AFAC, 1 )
00339                      CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00340 *
00341 *                    Factor the matrix and solve the system using CHPSV.
00342 *
00343                      SRNAMT = 'CHPSV '
00344                      CALL CHPSV( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
00345      $                           INFO )
00346 *
00347 *                    Adjust the expected value of INFO to account for
00348 *                    pivoting.
00349 *
00350                      K = IZERO
00351                      IF( K.GT.0 ) THEN
00352   100                   CONTINUE
00353                         IF( IWORK( K ).LT.0 ) THEN
00354                            IF( IWORK( K ).NE.-K ) THEN
00355                               K = -IWORK( K )
00356                               GO TO 100
00357                            END IF
00358                         ELSE IF( IWORK( K ).NE.K ) THEN
00359                            K = IWORK( K )
00360                            GO TO 100
00361                         END IF
00362                      END IF
00363 *
00364 *                    Check error code from CHPSV .
00365 *
00366                      IF( INFO.NE.K ) THEN
00367                         CALL ALAERH( PATH, 'CHPSV ', INFO, K, UPLO, N,
00368      $                               N, -1, -1, NRHS, IMAT, NFAIL,
00369      $                               NERRS, NOUT )
00370                         GO TO 120
00371                      ELSE IF( INFO.NE.0 ) THEN
00372                         GO TO 120
00373                      END IF
00374 *
00375 *                    Reconstruct matrix from factors and compute
00376 *                    residual.
00377 *
00378                      CALL CHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA,
00379      $                            RWORK, RESULT( 1 ) )
00380 *
00381 *                    Compute residual of the computed solution.
00382 *
00383                      CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00384                      CALL CPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
00385      $                            RWORK, RESULT( 2 ) )
00386 *
00387 *                    Check solution from generated exact solution.
00388 *
00389                      CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00390      $                            RESULT( 3 ) )
00391                      NT = 3
00392 *
00393 *                    Print information about the tests that did not pass
00394 *                    the threshold.
00395 *
00396                      DO 110 K = 1, NT
00397                         IF( RESULT( K ).GE.THRESH ) THEN
00398                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00399      $                        CALL ALADHD( NOUT, PATH )
00400                            WRITE( NOUT, FMT = 9999 )'CHPSV ', UPLO, N,
00401      $                        IMAT, K, RESULT( K )
00402                            NFAIL = NFAIL + 1
00403                         END IF
00404   110                CONTINUE
00405                      NRUN = NRUN + NT
00406   120                CONTINUE
00407                   END IF
00408 *
00409 *                 --- Test CHPSVX ---
00410 *
00411                   IF( IFACT.EQ.2 .AND. NPP.GT.0 )
00412      $               CALL CLASET( 'Full', NPP, 1, CMPLX( ZERO ),
00413      $                            CMPLX( ZERO ), AFAC, NPP )
00414                   CALL CLASET( 'Full', N, NRHS, CMPLX( ZERO ),
00415      $                         CMPLX( ZERO ), X, LDA )
00416 *
00417 *                 Solve the system and compute the condition number and
00418 *                 error bounds using CHPSVX.
00419 *
00420                   SRNAMT = 'CHPSVX'
00421                   CALL CHPSVX( FACT, UPLO, N, NRHS, A, AFAC, IWORK, B,
00422      $                         LDA, X, LDA, RCOND, RWORK,
00423      $                         RWORK( NRHS+1 ), WORK, RWORK( 2*NRHS+1 ),
00424      $                         INFO )
00425 *
00426 *                 Adjust the expected value of INFO to account for
00427 *                 pivoting.
00428 *
00429                   K = IZERO
00430                   IF( K.GT.0 ) THEN
00431   130                CONTINUE
00432                      IF( IWORK( K ).LT.0 ) THEN
00433                         IF( IWORK( K ).NE.-K ) THEN
00434                            K = -IWORK( K )
00435                            GO TO 130
00436                         END IF
00437                      ELSE IF( IWORK( K ).NE.K ) THEN
00438                         K = IWORK( K )
00439                         GO TO 130
00440                      END IF
00441                   END IF
00442 *
00443 *                 Check the error code from CHPSVX.
00444 *
00445                   IF( INFO.NE.K ) THEN
00446                      CALL ALAERH( PATH, 'CHPSVX', INFO, K, FACT // UPLO,
00447      $                            N, N, -1, -1, NRHS, IMAT, NFAIL,
00448      $                            NERRS, NOUT )
00449                      GO TO 150
00450                   END IF
00451 *
00452                   IF( INFO.EQ.0 ) THEN
00453                      IF( IFACT.GE.2 ) THEN
00454 *
00455 *                       Reconstruct matrix from factors and compute
00456 *                       residual.
00457 *
00458                         CALL CHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA,
00459      $                               RWORK( 2*NRHS+1 ), RESULT( 1 ) )
00460                         K1 = 1
00461                      ELSE
00462                         K1 = 2
00463                      END IF
00464 *
00465 *                    Compute residual of the computed solution.
00466 *
00467                      CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
00468                      CALL CPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
00469      $                            RWORK( 2*NRHS+1 ), RESULT( 2 ) )
00470 *
00471 *                    Check solution from generated exact solution.
00472 *
00473                      CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00474      $                            RESULT( 3 ) )
00475 *
00476 *                    Check the error bounds from iterative refinement.
00477 *
00478                      CALL CPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA,
00479      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
00480      $                            RESULT( 4 ) )
00481                   ELSE
00482                      K1 = 6
00483                   END IF
00484 *
00485 *                 Compare RCOND from CHPSVX with the computed value
00486 *                 in RCONDC.
00487 *
00488                   RESULT( 6 ) = SGET06( RCOND, RCONDC )
00489 *
00490 *                 Print information about the tests that did not pass
00491 *                 the threshold.
00492 *
00493                   DO 140 K = K1, 6
00494                      IF( RESULT( K ).GE.THRESH ) THEN
00495                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00496      $                     CALL ALADHD( NOUT, PATH )
00497                         WRITE( NOUT, FMT = 9998 )'CHPSVX', FACT, UPLO,
00498      $                     N, IMAT, K, RESULT( K )
00499                         NFAIL = NFAIL + 1
00500                      END IF
00501   140             CONTINUE
00502                   NRUN = NRUN + 7 - K1
00503 *
00504   150          CONTINUE
00505 *
00506   160       CONTINUE
00507   170    CONTINUE
00508   180 CONTINUE
00509 *
00510 *     Print a summary of the results.
00511 *
00512       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
00513 *
00514  9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
00515      $      ', test ', I2, ', ratio =', G12.5 )
00516  9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N =', I5,
00517      $      ', type ', I2, ', test ', I2, ', ratio =', G12.5 )
00518       RETURN
00519 *
00520 *     End of CDRVHP
00521 *
00522       END