LAPACK: TESTING/LIN/zdrvrfp.f Source File

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