LAPACK: TESTING/LIN/ddrvrf4.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DDRVRF4( NOUT, NN, NVAL, THRESH, C1, C2, LDC, CRF, A,
00002      +                    LDA, D_WORK_DLANGE )
00003 *
00004 *  -- LAPACK test routine (version 3.2.0) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2008
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            LDA, LDC, NN, NOUT
00010       DOUBLE PRECISION   THRESH
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            NVAL( NN )
00014       DOUBLE PRECISION   A( LDA, * ), C1( LDC, * ), C2( LDC, *),
00015      +                   CRF( * ), D_WORK_DLANGE( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DDRVRF4 tests the LAPACK RFP routines:
00022 *      DSFRK
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  NOUT          (input) INTEGER
00028 *                The unit number for output.
00029 *
00030 *  NN            (input) INTEGER
00031 *                The number of values of N contained in the vector NVAL.
00032 *
00033 *  NVAL          (input) INTEGER array, dimension (NN)
00034 *                The values of the matrix dimension N.
00035 *
00036 *  THRESH        (input) DOUBLE PRECISION
00037 *                The threshold value for the test ratios.  A result is
00038 *                included in the output file if RESULT >= THRESH.  To
00039 *                have every test ratio printed, use THRESH = 0.
00040 *
00041 *  C1            (workspace) DOUBLE PRECISION array,
00042 *                dimension (LDC,NMAX)
00043 *
00044 *  C2            (workspace) DOUBLE PRECISION array,
00045 *                dimension (LDC,NMAX)
00046 *
00047 *  LDC           (input) INTEGER
00048 *                The leading dimension of the array A.
00049 *                LDA >= max(1,NMAX).
00050 *
00051 *  CRF           (workspace) DOUBLE PRECISION array,
00052 *                dimension ((NMAX*(NMAX+1))/2).
00053 *
00054 *  A             (workspace) DOUBLE PRECISION array,
00055 *                dimension (LDA,NMAX)
00056 *
00057 *  LDA           (input) INTEGER
00058 *                The leading dimension of the array A.  LDA >= max(1,NMAX).
00059 *
00060 *  D_WORK_DLANGE (workspace) DOUBLE PRECISION array, dimension (NMAX)
00061 *
00062 *  =====================================================================
00063 *     ..
00064 *     .. Parameters ..
00065       DOUBLE PRECISION   ZERO, ONE
00066       PARAMETER          ( ZERO = 0.0D+0, ONE  = 1.0D+0 )
00067       INTEGER            NTESTS
00068       PARAMETER          ( NTESTS = 1 )
00069 *     ..
00070 *     .. Local Scalars ..
00071       CHARACTER          UPLO, CFORM, TRANS
00072       INTEGER            I, IFORM, IIK, IIN, INFO, IUPLO, J, K, N,
00073      +                   NFAIL, NRUN, IALPHA, ITRANS
00074       DOUBLE PRECISION   ALPHA, BETA, EPS, NORMA, NORMC
00075 *     ..
00076 *     .. Local Arrays ..
00077       CHARACTER          UPLOS( 2 ), FORMS( 2 ), TRANSS( 2 )
00078       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00079       DOUBLE PRECISION   RESULT( NTESTS )
00080 *     ..
00081 *     .. External Functions ..
00082       DOUBLE PRECISION   DLAMCH, DLARND, DLANGE
00083       EXTERNAL           DLAMCH, DLARND, DLANGE
00084 *     ..
00085 *     .. External Subroutines ..
00086       EXTERNAL           DSYRK, DSFRK, DTFTTR, DTRTTF
00087 *     ..
00088 *     .. Intrinsic Functions ..
00089       INTRINSIC          ABS, MAX
00090 *     ..
00091 *     .. Scalars in Common ..
00092       CHARACTER*32       SRNAMT
00093 *     ..
00094 *     .. Common blocks ..
00095       COMMON             / SRNAMC / SRNAMT
00096 *     ..
00097 *     .. Data statements ..
00098       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00099       DATA               UPLOS  / 'U', 'L' /
00100       DATA               FORMS  / 'N', 'T' /
00101       DATA               TRANSS / 'N', 'T' /
00102 *     ..
00103 *     .. Executable Statements ..
00104 *
00105 *     Initialize constants and the random number seed.
00106 *
00107       NRUN = 0
00108       NFAIL = 0
00109       INFO = 0
00110       DO 10 I = 1, 4
00111          ISEED( I ) = ISEEDY( I )
00112    10 CONTINUE
00113       EPS = DLAMCH( 'Precision' )
00114 *
00115       DO 150 IIN = 1, NN
00116 *
00117          N = NVAL( IIN )
00118 *
00119          DO 140 IIK = 1, NN
00120 *
00121             K = NVAL( IIN )
00122 *
00123             DO 130 IFORM = 1, 2
00124 *
00125                CFORM = FORMS( IFORM )
00126 *
00127                DO 120 IUPLO = 1, 2
00128 *
00129                   UPLO = UPLOS( IUPLO )
00130 *
00131                   DO 110 ITRANS = 1, 2
00132 *
00133                      TRANS = TRANSS( ITRANS )
00134 *
00135                      DO 100 IALPHA = 1, 4
00136 *
00137                         IF ( IALPHA.EQ. 1) THEN
00138                            ALPHA = ZERO
00139                            BETA = ZERO
00140                         ELSE IF ( IALPHA.EQ. 2) THEN
00141                            ALPHA = ONE
00142                            BETA = ZERO
00143                         ELSE IF ( IALPHA.EQ. 3) THEN
00144                            ALPHA = ZERO
00145                            BETA = ONE
00146                         ELSE
00147                            ALPHA = DLARND( 2, ISEED )
00148                            BETA = DLARND( 2, ISEED )
00149                         END IF
00150 *
00151 *                       All the parameters are set:
00152 *                          CFORM, UPLO, TRANS, M, N,
00153 *                          ALPHA, and BETA
00154 *                       READY TO TEST!
00155 *
00156                         NRUN = NRUN + 1
00157 *
00158                         IF ( ITRANS.EQ.1 ) THEN
00159 *
00160 *                          In this case we are NOTRANS, so A is N-by-K
00161 *
00162                            DO J = 1, K
00163                               DO I = 1, N
00164                                  A( I, J) = DLARND( 2, ISEED )
00165                               END DO
00166                            END DO
00167 *
00168                            NORMA = DLANGE( 'I', N, K, A, LDA,
00169      +                                      D_WORK_DLANGE )
00170 *
00171  
00172                         ELSE
00173 *
00174 *                          In this case we are TRANS, so A is K-by-N
00175 *
00176                            DO J = 1,N 
00177                               DO I = 1, K
00178                                  A( I, J) = DLARND( 2, ISEED )
00179                               END DO
00180                            END DO
00181 *
00182                            NORMA = DLANGE( 'I', K, N, A, LDA,
00183      +                                      D_WORK_DLANGE )
00184 *
00185                         END IF
00186 *
00187 *                       Generate C1 our N--by--N symmetric matrix. 
00188 *                       Make sure C2 has the same upper/lower part,
00189 *                       (the one that we do not touch), so
00190 *                       copy the initial C1 in C2 in it.
00191 *
00192                         DO J = 1, N
00193                            DO I = 1, N
00194                               C1( I, J) = DLARND( 2, ISEED )
00195                               C2(I,J) = C1(I,J)
00196                            END DO
00197                         END DO
00198 *
00199 *                       (See comment later on for why we use DLANGE and
00200 *                       not DLANSY for C1.)
00201 *
00202                         NORMC = DLANGE( 'I', N, N, C1, LDC,
00203      +                                      D_WORK_DLANGE )
00204 *
00205                         SRNAMT = 'DTRTTF'
00206                         CALL DTRTTF( CFORM, UPLO, N, C1, LDC, CRF,
00207      +                               INFO )
00208 *
00209 *                       call dsyrk the BLAS routine -> gives C1
00210 *
00211                         SRNAMT = 'DSYRK '
00212                         CALL DSYRK( UPLO, TRANS, N, K, ALPHA, A, LDA,
00213      +                              BETA, C1, LDC )
00214 *
00215 *                       call dsfrk the RFP routine -> gives CRF
00216 *
00217                         SRNAMT = 'DSFRK '
00218                         CALL DSFRK( CFORM, UPLO, TRANS, N, K, ALPHA, A,
00219      +                              LDA, BETA, CRF )
00220 *
00221 *                       convert CRF in full format -> gives C2
00222 *
00223                         SRNAMT = 'DTFTTR'
00224                         CALL DTFTTR( CFORM, UPLO, N, CRF, C2, LDC,
00225      +                               INFO )
00226 *
00227 *                       compare C1 and C2
00228 *
00229                         DO J = 1, N
00230                            DO I = 1, N
00231                               C1(I,J) = C1(I,J)-C2(I,J)
00232                            END DO
00233                         END DO
00234 *
00235 *                       Yes, C1 is symmetric so we could call DLANSY,
00236 *                       but we want to check the upper part that is
00237 *                       supposed to be unchanged and the diagonal that
00238 *                       is supposed to be real -> DLANGE
00239 *
00240                         RESULT(1) = DLANGE( 'I', N, N, C1, LDC,
00241      +                                      D_WORK_DLANGE )
00242                         RESULT(1) = RESULT(1) 
00243      +                              / MAX( ABS( ALPHA ) * NORMA
00244      +                                   + ABS( BETA ) , ONE )
00245      +                              / MAX( N , 1 ) / EPS
00246 *
00247                         IF( RESULT(1).GE.THRESH ) THEN
00248                            IF( NFAIL.EQ.0 ) THEN
00249                               WRITE( NOUT, * )
00250                               WRITE( NOUT, FMT = 9999 )
00251                            END IF
00252                            WRITE( NOUT, FMT = 9997 ) 'DSFRK', 
00253      +                        CFORM, UPLO, TRANS, N, K, RESULT(1)
00254                            NFAIL = NFAIL + 1
00255                         END IF
00256 *
00257   100                CONTINUE
00258   110             CONTINUE
00259   120          CONTINUE
00260   130       CONTINUE
00261   140    CONTINUE
00262   150 CONTINUE
00263 *
00264 *     Print a summary of the results.
00265 *
00266       IF ( NFAIL.EQ.0 ) THEN
00267          WRITE( NOUT, FMT = 9996 ) 'DSFRK', NRUN
00268       ELSE
00269          WRITE( NOUT, FMT = 9995 ) 'DSFRK', NFAIL, NRUN
00270       END IF
00271 *
00272  9999 FORMAT( 1X, 
00273 ' *** Error(s) or Failure(s) while testing DSFRK      +         ***')
00274  9997 FORMAT( 1X, '     Failure in ',A5,', CFORM=''',A1,''',',
00275      + ' UPLO=''',A1,''',',' TRANS=''',A1,''',', ' N=',I3,', K =', I3,
00276      + ', test=',G12.5)
00277  9996 FORMAT( 1X, 'All tests for ',A5,' auxiliary routine passed the ',
00278      +        'threshold (',I5,' tests run)')
00279  9995 FORMAT( 1X, A6, ' auxiliary routine:',I5,' out of ',I5,
00280      +        ' tests failed to pass the threshold')
00281 *
00282       RETURN
00283 *
00284 *     End of DDRVRF4
00285 *
00286       END