LAPACK: TESTING/LIN/cdrvrfp.f Source File

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