LAPACK: TESTING/LIN/ddrvpp.f Source File

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