LAPACK: TESTING/LIN/sdrvrfp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL,
00002      +              THRESH, A, ASAV, AFAC, AINV, B,
00003      +              BSAV, XACT, X, ARF, ARFINV,
00004      +              S_WORK_SLATMS, S_WORK_SPOT01, S_TEMP_SPOT02,
00005      +              S_TEMP_SPOT03, S_WORK_SLANSY,
00006      +              S_WORK_SPOT02, S_WORK_SPOT03 )
00007 *
00008 *  -- LAPACK test routine (version 3.2.0) --
00009 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00010 *     November 2008
00011 *
00012 *     .. Scalar Arguments ..
00013       INTEGER            NN, NNS, NNT, NOUT
00014       REAL               THRESH
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
00018       REAL               A( * )
00019       REAL               AINV( * )
00020       REAL               ASAV( * )
00021       REAL               B( * )
00022       REAL               BSAV( * )
00023       REAL               AFAC( * )
00024       REAL               ARF( * )
00025       REAL               ARFINV( * )
00026       REAL               XACT( * )
00027       REAL               X( * )
00028       REAL               S_WORK_SLATMS( * )
00029       REAL               S_WORK_SPOT01( * )
00030       REAL               S_TEMP_SPOT02( * )
00031       REAL               S_TEMP_SPOT03( * )
00032       REAL               S_WORK_SLANSY( * )
00033       REAL               S_WORK_SPOT02( * )
00034       REAL               S_WORK_SPOT03( * )
00035 *     ..
00036 *
00037 *  Purpose
00038 *  =======
00039 *
00040 *  SDRVRFP tests the LAPACK RFP routines:
00041 *      SPFTRF, SPFTRS, and SPFTRI.
00042 *
00043 *  This testing routine follow the same tests as DDRVPO (test for the full
00044 *  format Symmetric Positive Definite solver).
00045 *
00046 *  The tests are performed in Full Format, convertion back and forth from
00047 *  full format to RFP format are performed using the routines STRTTF and
00048 *  STFTTR.
00049 *
00050 *  First, a specific matrix A of size N is created. There is nine types of 
00051 *  different matrixes possible.
00052 *   1. Diagonal                        6. Random, CNDNUM = sqrt(0.1/EPS)
00053 *   2. Random, CNDNUM = 2              7. Random, CNDNUM = 0.1/EPS
00054 *  *3. First row and column zero       8. Scaled near underflow
00055 *  *4. Last row and column zero        9. Scaled near overflow
00056 *  *5. Middle row and column zero
00057 *  (* - tests error exits from SPFTRF, no test ratios are computed)
00058 *  A solution XACT of size N-by-NRHS is created and the associated right
00059 *  hand side B as well. Then SPFTRF is called to compute L (or U), the
00060 *  Cholesky factor of A. Then L (or U) is used to solve the linear system
00061 *  of equations AX = B. This gives X. Then L (or U) is used to compute the
00062 *  inverse of A, AINV. The following four tests are then performed:
00063 *  (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00064 *      norm( U'*U - A ) / ( N * norm(A) * EPS ),
00065 *  (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
00066 *  (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
00067 *  (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
00068 *  where EPS is the machine precision, RCOND the condition number of A, and
00069 *  norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4).
00070 *  Errors occur when INFO parameter is not as expected. Failures occur when
00071 *  a test ratios is greater than THRES.
00072 *
00073 *  Arguments
00074 *  =========
00075 *
00076 *  NOUT          (input) INTEGER
00077 *                The unit number for output.
00078 *
00079 *  NN            (input) INTEGER
00080 *                The number of values of N contained in the vector NVAL.
00081 *
00082 *  NVAL          (input) INTEGER array, dimension (NN)
00083 *                The values of the matrix dimension N.
00084 *
00085 *  NNS           (input) INTEGER
00086 *                The number of values of NRHS contained in the vector NSVAL.
00087 *
00088 *  NSVAL         (input) INTEGER array, dimension (NNS)
00089 *                The values of the number of right-hand sides NRHS.
00090 *
00091 *  NNT           (input) INTEGER
00092 *                The number of values of MATRIX TYPE contained in the vector NTVAL.
00093 *
00094 *  NTVAL         (input) INTEGER array, dimension (NNT)
00095 *                The values of matrix type (between 0 and 9 for PO/PP/PF matrices).
00096 *
00097 *  THRESH        (input) REAL
00098 *                The threshold value for the test ratios.  A result is
00099 *                included in the output file if RESULT >= THRESH.  To have
00100 *                every test ratio printed, use THRESH = 0.
00101 *
00102 *  A             (workspace) REAL array, dimension (NMAX*NMAX)
00103 *
00104 *  ASAV          (workspace) REAL array, dimension (NMAX*NMAX)
00105 *
00106 *  AFAC          (workspace) REAL array, dimension (NMAX*NMAX)
00107 *
00108 *  AINV          (workspace) REAL array, dimension (NMAX*NMAX)
00109 *
00110 *  B             (workspace) REAL array, dimension (NMAX*MAXRHS)
00111 *
00112 *  BSAV          (workspace) REAL array, dimension (NMAX*MAXRHS)
00113 *
00114 *  XACT          (workspace) REAL array, dimension (NMAX*MAXRHS)
00115 *
00116 *  X             (workspace) REAL array, dimension (NMAX*MAXRHS)
00117 *
00118 *  ARF           (workspace) REAL array, dimension ((NMAX*(NMAX+1))/2)
00119 *
00120 *  ARFINV        (workspace) REAL array, dimension ((NMAX*(NMAX+1))/2)
00121 *
00122 *  S_WORK_SLATMS (workspace) REAL array, dimension ( 3*NMAX )
00123 *
00124 *  S_WORK_SPOT01 (workspace) REAL array, dimension ( NMAX )
00125 *
00126 *  S_TEMP_SPOT02 (workspace) REAL array, dimension ( NMAX*MAXRHS )
00127 *
00128 *  S_TEMP_SPOT03 (workspace) REAL array, dimension ( NMAX*NMAX )
00129 *
00130 *  S_WORK_SLATMS (workspace) REAL array, dimension ( NMAX )
00131 *
00132 *  S_WORK_SLANSY (workspace) REAL array, dimension ( NMAX )
00133 *
00134 *  S_WORK_SPOT02 (workspace) REAL array, dimension ( NMAX )
00135 *
00136 *  S_WORK_SPOT03 (workspace) REAL array, dimension ( NMAX )
00137 *
00138 *  =====================================================================
00139 *
00140 *     .. Parameters ..
00141       REAL               ONE, ZERO
00142       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00143       INTEGER            NTESTS
00144       PARAMETER          ( NTESTS = 4 )
00145 *     ..
00146 *     .. Local Scalars ..
00147       LOGICAL            ZEROT
00148       INTEGER            I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL,
00149      +                   NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN,
00150      +                   IIT, IIS
00151       CHARACTER          DIST, CTYPE, UPLO, CFORM
00152       INTEGER            KL, KU, MODE
00153       REAL               ANORM, AINVNM, CNDNUM, RCONDC
00154 *     ..
00155 *     .. Local Arrays ..
00156       CHARACTER          UPLOS( 2 ), FORMS( 2 )
00157       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00158       REAL               RESULT( NTESTS )
00159 *     ..
00160 *     .. External Functions ..
00161       REAL               SLANSY
00162       EXTERNAL           SLANSY
00163 *     ..
00164 *     .. External Subroutines ..
00165       EXTERNAL           ALADHD, ALAERH, ALASVM, SGET04, STFTTR, SLACPY,
00166      +                   SLARHS, SLATB4, SLATMS, SPFTRI, SPFTRF, SPFTRS,
00167      +                   SPOT01, SPOT02, SPOT03, SPOTRI, SPOTRF, STRTTF
00168 *     ..
00169 *     .. Scalars in Common ..
00170       CHARACTER*32       SRNAMT
00171 *     ..
00172 *     .. Common blocks ..
00173       COMMON             / SRNAMC / SRNAMT
00174 *     ..
00175 *     .. Data statements ..
00176       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00177       DATA               UPLOS / 'U', 'L' /
00178       DATA               FORMS / 'N', 'T' /
00179 *     ..
00180 *     .. Executable Statements ..
00181 *
00182 *     Initialize constants and the random number seed.
00183 *
00184       NRUN = 0
00185       NFAIL = 0
00186       NERRS = 0
00187       DO 10 I = 1, 4
00188          ISEED( I ) = ISEEDY( I )
00189    10 CONTINUE
00190 *
00191       DO 130 IIN = 1, NN
00192 *
00193          N = NVAL( IIN )
00194          LDA = MAX( N, 1 )
00195          LDB = MAX( N, 1 )
00196 *
00197          DO 980 IIS = 1, NNS
00198 *
00199             NRHS = NSVAL( IIS )
00200 *
00201             DO 120 IIT = 1, NNT
00202 *
00203                IMAT = NTVAL( IIT )
00204 *
00205 *              If N.EQ.0, only consider the first type
00206 *
00207                IF( N.EQ.0 .AND. IIT.GT.1 ) GO TO 120
00208 *
00209 *              Skip types 3, 4, or 5 if the matrix size is too small.
00210 *
00211                IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120
00212                IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120
00213 *
00214 *              Do first for UPLO = 'U', then for UPLO = 'L'
00215 *
00216                DO 110 IUPLO = 1, 2
00217                   UPLO = UPLOS( IUPLO )
00218 *
00219 *                 Do first for CFORM = 'N', then for CFORM = 'C'
00220 *
00221                   DO 100 IFORM = 1, 2
00222                      CFORM = FORMS( IFORM )
00223 *
00224 *                    Set up parameters with SLATB4 and generate a test
00225 *                    matrix with SLATMS.
00226 *
00227                      CALL SLATB4( 'SPO', IMAT, N, N, CTYPE, KL, KU,
00228      +                            ANORM, MODE, CNDNUM, DIST )
00229 *
00230                      SRNAMT = 'SLATMS'
00231                      CALL SLATMS( N, N, DIST, ISEED, CTYPE,
00232      +                            S_WORK_SLATMS,
00233      +                            MODE, CNDNUM, ANORM, KL, KU, UPLO, A,
00234      +                            LDA, S_WORK_SLATMS, INFO )
00235 *
00236 *                    Check error code from SLATMS.
00237 *
00238                      IF( INFO.NE.0 ) THEN
00239                         CALL ALAERH( 'SPF', 'SLATMS', INFO, 0, UPLO, N,
00240      +                               N, -1, -1, -1, IIT, NFAIL, NERRS,
00241      +                               NOUT )
00242                         GO TO 100
00243                      END IF
00244 *
00245 *                    For types 3-5, zero one row and column of the matrix to
00246 *                    test that INFO is returned correctly.
00247 *
00248                      ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
00249                      IF( ZEROT ) THEN
00250                         IF( IIT.EQ.3 ) THEN
00251                            IZERO = 1
00252                         ELSE IF( IIT.EQ.4 ) THEN
00253                            IZERO = N
00254                         ELSE
00255                            IZERO = N / 2 + 1
00256                         END IF
00257                         IOFF = ( IZERO-1 )*LDA
00258 *
00259 *                       Set row and column IZERO of A to 0.
00260 *
00261                         IF( IUPLO.EQ.1 ) THEN
00262                            DO 20 I = 1, IZERO - 1
00263                               A( IOFF+I ) = ZERO
00264    20                      CONTINUE
00265                            IOFF = IOFF + IZERO
00266                            DO 30 I = IZERO, N
00267                               A( IOFF ) = ZERO
00268                               IOFF = IOFF + LDA
00269    30                      CONTINUE
00270                         ELSE
00271                            IOFF = IZERO
00272                            DO 40 I = 1, IZERO - 1
00273                               A( IOFF ) = ZERO
00274                               IOFF = IOFF + LDA
00275    40                      CONTINUE
00276                            IOFF = IOFF - IZERO
00277                            DO 50 I = IZERO, N
00278                               A( IOFF+I ) = ZERO
00279    50                      CONTINUE
00280                         END IF
00281                      ELSE
00282                         IZERO = 0
00283                      END IF
00284 *
00285 *                    Save a copy of the matrix A in ASAV.
00286 *
00287                      CALL SLACPY( UPLO, N, N, A, LDA, ASAV, LDA )
00288 *
00289 *                    Compute the condition number of A (RCONDC).
00290 *
00291                      IF( ZEROT ) THEN
00292                         RCONDC = ZERO
00293                      ELSE
00294 *
00295 *                       Compute the 1-norm of A.
00296 *
00297                         ANORM = SLANSY( '1', UPLO, N, A, LDA,
00298      +                         S_WORK_SLANSY )
00299 *
00300 *                       Factor the matrix A.
00301 *
00302                         CALL SPOTRF( UPLO, N, A, LDA, INFO )
00303 *
00304 *                       Form the inverse of A.
00305 *
00306                         CALL SPOTRI( UPLO, N, A, LDA, INFO )
00307 *
00308 *                       Compute the 1-norm condition number of A.
00309 *
00310                         AINVNM = SLANSY( '1', UPLO, N, A, LDA,
00311      +                           S_WORK_SLANSY )
00312                         RCONDC = ( ONE / ANORM ) / AINVNM
00313 *
00314 *                       Restore the matrix A.
00315 *
00316                         CALL SLACPY( UPLO, N, N, ASAV, LDA, A, LDA )
00317 *
00318                      END IF
00319 *
00320 *                    Form an exact solution and set the right hand side.
00321 *
00322                      SRNAMT = 'SLARHS'
00323                      CALL SLARHS( 'SPO', 'N', UPLO, ' ', N, N, KL, KU,
00324      +                            NRHS, A, LDA, XACT, LDA, B, LDA,
00325      +                            ISEED, INFO )
00326                      CALL SLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA )
00327 *
00328 *                    Compute the L*L' or U'*U factorization of the
00329 *                    matrix and solve the system.
00330 *
00331                      CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
00332                      CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDB )
00333 *
00334                      SRNAMT = 'STRTTF'
00335                      CALL STRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO )
00336                      SRNAMT = 'SPFTRF'
00337                      CALL SPFTRF( CFORM, UPLO, N, ARF, INFO )
00338 *
00339 *                    Check error code from SPFTRF.
00340 *
00341                      IF( INFO.NE.IZERO ) THEN
00342 *
00343 *                       LANGOU: there is a small hick here: IZERO should
00344 *                       always be INFO however if INFO is ZERO, ALAERH does not
00345 *                       complain.
00346 *
00347                          CALL ALAERH( 'SPF', 'SPFSV ', INFO, IZERO,
00348      +                                UPLO, N, N, -1, -1, NRHS, IIT,
00349      +                                NFAIL, NERRS, NOUT )
00350                          GO TO 100
00351                       END IF
00352 *
00353 *                    Skip the tests if INFO is not 0.
00354 *
00355                      IF( INFO.NE.0 ) THEN
00356                         GO TO 100
00357                      END IF
00358 *
00359                      SRNAMT = 'SPFTRS'
00360                      CALL SPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB,
00361      +                            INFO )
00362 *
00363                      SRNAMT = 'STFTTR'
00364                      CALL STFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO )
00365 *
00366 *                    Reconstruct matrix from factors and compute
00367 *                    residual.
00368 *
00369                      CALL SLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA )
00370                      CALL SPOT01( UPLO, N, A, LDA, AFAC, LDA,
00371      +                             S_WORK_SPOT01, RESULT( 1 ) )
00372                      CALL SLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA )
00373 *
00374 *                    Form the inverse and compute the residual.
00375 *
00376                      IF(MOD(N,2).EQ.0)THEN 
00377                         CALL SLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV,
00378      +                               N+1 )
00379                      ELSE
00380                         CALL SLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV,
00381      +                               N )
00382                      END IF
00383 *
00384                      SRNAMT = 'SPFTRI'
00385                      CALL SPFTRI( CFORM, UPLO, N, ARFINV , INFO )
00386 *
00387                      SRNAMT = 'STFTTR'
00388                      CALL STFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA,
00389      +                            INFO )
00390 *
00391 *                    Check error code from SPFTRI.
00392 *
00393                      IF( INFO.NE.0 )
00394      +                  CALL ALAERH( 'SPO', 'SPFTRI', INFO, 0, UPLO, N,
00395      +                               N, -1, -1, -1, IMAT, NFAIL, NERRS,
00396      +                               NOUT )
00397 *
00398                      CALL SPOT03( UPLO, N, A, LDA, AINV, LDA,
00399      +                            S_TEMP_SPOT03, LDA, S_WORK_SPOT03,
00400      +                            RCONDC, RESULT( 2 ) )
00401 *
00402 *                    Compute residual of the computed solution.
00403 *
00404                      CALL SLACPY( 'Full', N, NRHS, B, LDA,
00405      +                            S_TEMP_SPOT02, LDA )
00406                      CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA,
00407      +                            S_TEMP_SPOT02, LDA, S_WORK_SPOT02,
00408      +                            RESULT( 3 ) )
00409 *
00410 *                    Check solution from generated exact solution.
00411  
00412                      CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
00413      +                         RESULT( 4 ) )
00414                      NT = 4
00415 *
00416 *                    Print information about the tests that did not
00417 *                    pass the threshold.
00418 *
00419                      DO 60 K = 1, NT
00420                         IF( RESULT( K ).GE.THRESH ) THEN
00421                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
00422      +                        CALL ALADHD( NOUT, 'SPF' )
00423                            WRITE( NOUT, FMT = 9999 )'SPFSV ', UPLO,
00424      +                            N, IIT, K, RESULT( K )
00425                            NFAIL = NFAIL + 1
00426                         END IF
00427    60                CONTINUE
00428                      NRUN = NRUN + NT
00429   100             CONTINUE
00430   110          CONTINUE
00431   120       CONTINUE
00432   980    CONTINUE
00433   130 CONTINUE
00434 *
00435 *     Print a summary of the results.
00436 *
00437       CALL ALASVM( 'SPF', NOUT, NFAIL, NRUN, NERRS )
00438 *
00439  9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1,
00440      +      ', test(', I1, ')=', G12.5 )
00441 *
00442       RETURN
00443 *
00444 *     End of SDRVRFP
00445 *
00446       END