LAPACK: TESTING/LIN/cdrvpp.f Source File

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