LAPACK: TESTING/LIN/cdrvpo.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CDRVPO( 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 *  CDRVPO tests the driver routines CPOSV 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)
00058 *
00059 *  AFAC    (workspace) COMPLEX array, dimension (NMAX*NMAX)
00060 *
00061 *  ASAV    (workspace) COMPLEX array, dimension (NMAX*NMAX)
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, 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, NB, NBMIN,
00097      $                   NERRS, NFACT, NFAIL, NIMAT, NRUN, NT
00098       REAL               AINVNM, AMAX, ANORM, CNDNUM, RCOND, RCONDC,
00099      $                   ROLDC, SCOND
00100 *     ..
00101 *     .. Local Arrays ..
00102       CHARACTER          EQUEDS( 2 ), FACTS( 3 ), UPLOS( 2 )
00103       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00104       REAL               RESULT( NTESTS )
00105 *     ..
00106 *     .. External Functions ..
00107       LOGICAL            LSAME
00108       REAL               CLANHE, SGET06
00109       EXTERNAL           LSAME, CLANHE, SGET06
00110 *     ..
00111 *     .. External Subroutines ..
00112       EXTERNAL           ALADHD, ALAERH, ALASVM, CERRVX, CGET04, CLACPY,
00113      $                   CLAIPD, CLAQHE, CLARHS, CLASET, CLATB4, CLATMS,
00114      $                   CPOEQU, CPOSV, CPOSVX, CPOT01, CPOT02, CPOT05,
00115      $                   CPOTRF, CPOTRI, XLAENV
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' /
00132       DATA               FACTS / 'F', 'N', 'E' /
00133       DATA               EQUEDS / 'N', 'Y' /
00134 *     ..
00135 *     .. Executable Statements ..
00136 *
00137 *     Initialize constants and the random number seed.
00138 *
00139       PATH( 1: 1 ) = 'Complex precision'
00140       PATH( 2: 3 ) = 'PO'
00141       NRUN = 0
00142       NFAIL = 0
00143       NERRS = 0
00144       DO 10 I = 1, 4
00145          ISEED( I ) = ISEEDY( I )
00146    10 CONTINUE
00147 *
00148 *     Test the error exits
00149 *
00150       IF( TSTERR )
00151      $   CALL CERRVX( PATH, NOUT )
00152       INFOT = 0
00153 *
00154 *     Set the block size and minimum block size for testing.
00155 *
00156       NB = 1
00157       NBMIN = 2
00158       CALL XLAENV( 1, NB )
00159       CALL XLAENV( 2, NBMIN )
00160 *
00161 *     Do for each value of N in NVAL
00162 *
00163       DO 130 IN = 1, NN
00164          N = NVAL( IN )
00165          LDA = MAX( N, 1 )
00166          XTYPE = 'N'
00167          NIMAT = NTYPES
00168          IF( N.LE.0 )
00169      $      NIMAT = 1
00170 *
00171          DO 120 IMAT = 1, NIMAT
00172 *
00173 *           Do the tests only if DOTYPE( IMAT ) is true.
00174 *
00175             IF( .NOT.DOTYPE( IMAT ) )
00176      $         GO TO 120
00177 *
00178 *           Skip types 3, 4, or 5 if the matrix size is too small.
00179 *
00180             ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
00181             IF( ZEROT .AND. N.LT.IMAT-2 )
00182      $         GO TO 120
00183 *
00184 *           Do first for UPLO = 'U', then for UPLO = 'L'
00185 *
00186             DO 110 IUPLO = 1, 2
00187                UPLO = UPLOS( IUPLO )
00188 *
00189 *              Set up parameters with CLATB4 and generate a test matrix
00190 *              with CLATMS.
00191 *
00192                CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
00193      $                      CNDNUM, DIST )
00194 *
00195                SRNAMT = 'CLATMS'
00196                CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
00197      $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
00198      $                      INFO )
00199 *
00200 *              Check error code from CLATMS.
00201 *
00202                IF( INFO.NE.0 ) THEN
00203                   CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1,
00204      $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
00205                   GO TO 110
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                   IOFF = ( IZERO-1 )*LDA
00220 *
00221 *                 Set row and column IZERO of A to 0.
00222 *
00223                   IF( IUPLO.EQ.1 ) THEN
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 + LDA
00231    30                CONTINUE
00232                   ELSE
00233                      IOFF = IZERO
00234                      DO 40 I = 1, IZERO - 1
00235                         A( IOFF ) = ZERO
00236                         IOFF = IOFF + LDA
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 *              Set the imaginary part of the diagonals.
00248 *
00249                CALL CLAIPD( N, A, LDA+1, 0 )
00250 *
00251 *              Save a copy of the matrix A in ASAV.
00252 *
00253                CALL CLACPY( UPLO, N, N, A, LDA, ASAV, LDA )
00254 *
00255                DO 100 IEQUED = 1, 2
00256                   EQUED = EQUEDS( IEQUED )
00257                   IF( IEQUED.EQ.1 ) THEN
00258                      NFACT = 3
00259                   ELSE
00260                      NFACT = 1
00261                   END IF
00262 *
00263                   DO 90 IFACT = 1, NFACT
00264                      FACT = FACTS( IFACT )
00265                      PREFAC = LSAME( FACT, 'F' )
00266                      NOFACT = LSAME( FACT, 'N' )
00267                      EQUIL = LSAME( FACT, 'E' )
00268 *
00269                      IF( ZEROT ) THEN
00270                         IF( PREFAC )
00271      $                     GO TO 90
00272                         RCONDC = ZERO
00273 *
00274                      ELSE IF( .NOT.LSAME( FACT, 'N' ) ) THEN
00275 *
00276 *                       Compute the condition number for comparison with
00277 *                       the value returned by CPOSVX (FACT = 'N' reuses
00278 *                       the condition number from the previous iteration
00279 *                       with FACT = 'F').
00280 *
00281                         CALL CLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA )
00282                         IF( EQUIL .OR. IEQUED.GT.1 ) THEN
00283 *
00284 *                          Compute row and column scale factors to
00285 *                          equilibrate the matrix A.
00286 *
00287                            CALL CPOEQU( N, AFAC, LDA, S, SCOND, AMAX,
00288      $                                  INFO )
00289                            IF( INFO.EQ.0 .AND. N.GT.0 ) THEN
00290                               IF( IEQUED.GT.1 )
00291      $                           SCOND = ZERO
00292 *
00293 *                             Equilibrate the matrix.
00294 *
00295                               CALL CLAQHE( UPLO, N, AFAC, LDA, S, SCOND,
00296      $                                     AMAX, EQUED )
00297                            END IF
00298                         END IF
00299 *
00300 *                       Save the condition number of the
00301 *                       non-equilibrated system for use in CGET04.
00302 *
00303                         IF( EQUIL )
00304      $                     ROLDC = RCONDC
00305 *
00306 *                       Compute the 1-norm of A.
00307 *
00308                         ANORM = CLANHE( '1', UPLO, N, AFAC, LDA, RWORK )
00309 *
00310 *                       Factor the matrix A.
00311 *
00312                         CALL CPOTRF( UPLO, N, AFAC, LDA, INFO )
00313 *
00314 *                       Form the inverse of A.
00315 *
00316                         CALL CLACPY( UPLO, N, N, AFAC, LDA, A, LDA )
00317                         CALL CPOTRI( UPLO, N, A, LDA, INFO )
00318 *
00319 *                       Compute the 1-norm condition number of A.
00320 *
00321                         AINVNM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
00322                         IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
00323                            RCONDC = ONE
00324                         ELSE
00325                            RCONDC = ( ONE / ANORM ) / AINVNM
00326                         END IF
00327                      END IF
00328 *
00329 *                    Restore the matrix A.
00330 *
00331                      CALL CLACPY( UPLO, N, N, ASAV, LDA, A, LDA )
00332 *
00333 *                    Form an exact solution and set the right hand side.
00334 *
00335                      SRNAMT = 'CLARHS'
00336                      CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
00337      $                            NRHS, A, LDA, XACT, LDA, B, LDA,
00338      $                            ISEED, INFO )
00339                      XTYPE = 'C'
00340                      CALL CLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA )
00341 *
00342                      IF( NOFACT ) THEN
00343 *
00344 *                       --- Test CPOSV  ---
00345 *
00346 *                       Compute the L*L' or U'*U factorization of the
00347 *                       matrix and solve the system.
00348 *
00349                         CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
00350                         CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
00351 *
00352                         SRNAMT = 'CPOSV '
00353                         CALL CPOSV( UPLO, N, NRHS, AFAC, LDA, X, LDA,
00354      $                              INFO )
00355 *
00356 *                       Check error code from CPOSV .
00357 *
00358                         IF( INFO.NE.IZERO ) THEN
00359                            CALL ALAERH( PATH, 'CPOSV ', INFO, IZERO,
00360      $                                  UPLO, N, N, -1, -1, NRHS, IMAT,
00361      $                                  NFAIL, NERRS, NOUT )
00362                            GO TO 70
00363                         ELSE IF( INFO.NE.0 ) THEN
00364                            GO TO 70
00365                         END IF
00366 *
00367 *                       Reconstruct matrix from factors and compute
00368 *                       residual.
00369 *
00370                         CALL CPOT01( UPLO, N, A, LDA, AFAC, LDA, RWORK,
00371      $                               RESULT( 1 ) )
00372 *
00373 *                       Compute residual of the computed solution.
00374 *
00375                         CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK,
00376      $                               LDA )
00377                         CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA,
00378      $                               WORK, LDA, RWORK, RESULT( 2 ) )
00379 *
00380 *                       Check solution from generated exact solution.
00381 *
00382                         CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00383      $                               RESULT( 3 ) )
00384                         NT = 3
00385 *
00386 *                       Print information about the tests that did not
00387 *                       pass the threshold.
00388 *
00389                         DO 60 K = 1, NT
00390                            IF( RESULT( K ).GE.THRESH ) THEN
00391                               IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00392      $                           CALL ALADHD( NOUT, PATH )
00393                               WRITE( NOUT, FMT = 9999 )'CPOSV ', UPLO,
00394      $                           N, IMAT, K, RESULT( K )
00395                               NFAIL = NFAIL + 1
00396                            END IF
00397    60                   CONTINUE
00398                         NRUN = NRUN + NT
00399    70                   CONTINUE
00400                      END IF
00401 *
00402 *                    --- Test CPOSVX ---
00403 *
00404                      IF( .NOT.PREFAC )
00405      $                  CALL CLASET( UPLO, N, N, CMPLX( ZERO ),
00406      $                               CMPLX( ZERO ), AFAC, LDA )
00407                      CALL CLASET( 'Full', N, NRHS, CMPLX( ZERO ),
00408      $                            CMPLX( ZERO ), X, LDA )
00409                      IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
00410 *
00411 *                       Equilibrate the matrix if FACT='F' and
00412 *                       EQUED='Y'.
00413 *
00414                         CALL CLAQHE( UPLO, N, A, LDA, S, SCOND, AMAX,
00415      $                               EQUED )
00416                      END IF
00417 *
00418 *                    Solve the system and compute the condition number
00419 *                    and error bounds using CPOSVX.
00420 *
00421                      SRNAMT = 'CPOSVX'
00422                      CALL CPOSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC,
00423      $                            LDA, EQUED, S, B, LDA, X, LDA, RCOND,
00424      $                            RWORK, RWORK( NRHS+1 ), WORK,
00425      $                            RWORK( 2*NRHS+1 ), INFO )
00426 *
00427 *                    Check the error code from CPOSVX.
00428 *
00429                      IF( INFO.NE.IZERO ) THEN
00430                         CALL ALAERH( PATH, 'CPOSVX', INFO, IZERO,
00431      $                               FACT // UPLO, N, N, -1, -1, NRHS,
00432      $                               IMAT, NFAIL, NERRS, NOUT )
00433                         GO TO 90
00434                      END IF
00435 *
00436                      IF( INFO.EQ.0 ) THEN
00437                         IF( .NOT.PREFAC ) THEN
00438 *
00439 *                          Reconstruct matrix from factors and compute
00440 *                          residual.
00441 *
00442                            CALL CPOT01( UPLO, N, A, LDA, AFAC, LDA,
00443      $                                  RWORK( 2*NRHS+1 ), RESULT( 1 ) )
00444                            K1 = 1
00445                         ELSE
00446                            K1 = 2
00447                         END IF
00448 *
00449 *                       Compute residual of the computed solution.
00450 *
00451                         CALL CLACPY( 'Full', N, NRHS, BSAV, LDA, WORK,
00452      $                               LDA )
00453                         CALL CPOT02( UPLO, N, NRHS, ASAV, LDA, X, LDA,
00454      $                               WORK, LDA, RWORK( 2*NRHS+1 ),
00455      $                               RESULT( 2 ) )
00456 *
00457 *                       Check solution from generated exact solution.
00458 *
00459                         IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED,
00460      $                      'N' ) ) ) THEN
00461                            CALL CGET04( N, NRHS, X, LDA, XACT, LDA,
00462      $                                  RCONDC, RESULT( 3 ) )
00463                         ELSE
00464                            CALL CGET04( N, NRHS, X, LDA, XACT, LDA,
00465      $                                  ROLDC, RESULT( 3 ) )
00466                         END IF
00467 *
00468 *                       Check the error bounds from iterative
00469 *                       refinement.
00470 *
00471                         CALL CPOT05( UPLO, N, NRHS, ASAV, LDA, B, LDA,
00472      $                               X, LDA, XACT, LDA, RWORK,
00473      $                               RWORK( NRHS+1 ), RESULT( 4 ) )
00474                      ELSE
00475                         K1 = 6
00476                      END IF
00477 *
00478 *                    Compare RCOND from CPOSVX with the computed value
00479 *                    in RCONDC.
00480 *
00481                      RESULT( 6 ) = SGET06( RCOND, RCONDC )
00482 *
00483 *                    Print information about the tests that did not pass
00484 *                    the threshold.
00485 *
00486                      DO 80 K = K1, 6
00487                         IF( RESULT( K ).GE.THRESH ) THEN
00488                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00489      $                        CALL ALADHD( NOUT, PATH )
00490                            IF( PREFAC ) THEN
00491                               WRITE( NOUT, FMT = 9997 )'CPOSVX', FACT,
00492      $                           UPLO, N, EQUED, IMAT, K, RESULT( K )
00493                            ELSE
00494                               WRITE( NOUT, FMT = 9998 )'CPOSVX', FACT,
00495      $                           UPLO, N, IMAT, K, RESULT( K )
00496                            END IF
00497                            NFAIL = NFAIL + 1
00498                         END IF
00499    80                CONTINUE
00500                      NRUN = NRUN + 7 - K1
00501    90             CONTINUE
00502   100          CONTINUE
00503   110       CONTINUE
00504   120    CONTINUE
00505   130 CONTINUE
00506 *
00507 *     Print a summary of the results.
00508 *
00509       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
00510 *
00511  9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I1,
00512      $      ', test(', I1, ')=', G12.5 )
00513  9998 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5,
00514      $      ', type ', I1, ', test(', I1, ')=', G12.5 )
00515  9997 FORMAT( 1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N=', I5,
00516      $      ', EQUED=''', A1, ''', type ', I1, ', test(', I1, ') =',
00517      $      G12.5 )
00518       RETURN
00519 *
00520 *     End of CDRVPO
00521 *
00522       END