LAPACK: TESTING/LIN/sdrvrf3.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SDRVRF3( NOUT, NN, NVAL, THRESH, A, LDA, ARF, B1, B2,
00002      +                    S_WORK_SLANGE, S_WORK_SGEQRF, TAU )
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, NN, NOUT
00010       REAL               THRESH
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            NVAL( NN )
00014       REAL               A( LDA, * ), ARF( * ), B1( LDA, * ),
00015      +                   B2( LDA, * ), S_WORK_SGEQRF( * ),
00016      +                   S_WORK_SLANGE( * ), TAU( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  SDRVRF3 tests the LAPACK RFP routines:
00023 *      STFSM
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  NOUT          (input) INTEGER
00029 *                The unit number for output.
00030 *
00031 *  NN            (input) INTEGER
00032 *                The number of values of N contained in the vector NVAL.
00033 *
00034 *  NVAL          (input) INTEGER array, dimension (NN)
00035 *                The values of the matrix dimension N.
00036 *
00037 *  THRESH        (input) REAL
00038 *                The threshold value for the test ratios.  A result is
00039 *                included in the output file if RESULT >= THRESH.  To have
00040 *                every test ratio printed, use THRESH = 0.
00041 *
00042 *  A             (workspace) REAL array, dimension (LDA,NMAX)
00043 *
00044 *  LDA           (input) INTEGER
00045 *                The leading dimension of the array A.  LDA >= max(1,NMAX).
00046 *
00047 *  ARF           (workspace) REAL array, dimension ((NMAX*(NMAX+1))/2).
00048 *
00049 *  B1            (workspace) REAL array, dimension (LDA,NMAX)
00050 *
00051 *  B2            (workspace) REAL array, dimension (LDA,NMAX)
00052 *
00053 *  S_WORK_SLANGE (workspace) REAL array, dimension (NMAX)
00054 *
00055 *  S_WORK_SGEQRF (workspace) REAL array, dimension (NMAX)
00056 *
00057 *  TAU           (workspace) REAL array, dimension (NMAX)
00058 *
00059 *  =====================================================================
00060 *     ..
00061 *     .. Parameters ..
00062       REAL               ZERO, ONE
00063       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) ,
00064      +                     ONE  = ( 1.0E+0, 0.0E+0 ) )
00065       INTEGER            NTESTS
00066       PARAMETER          ( NTESTS = 1 )
00067 *     ..
00068 *     .. Local Scalars ..
00069       CHARACTER          UPLO, CFORM, DIAG, TRANS, SIDE
00070       INTEGER            I, IFORM, IIM, IIN, INFO, IUPLO, J, M, N, NA,
00071      +                   NFAIL, NRUN, ISIDE, IDIAG, IALPHA, ITRANS
00072       REAL               EPS, ALPHA
00073 *     ..
00074 *     .. Local Arrays ..
00075       CHARACTER          UPLOS( 2 ), FORMS( 2 ), TRANSS( 2 ),
00076      +                   DIAGS( 2 ), SIDES( 2 )
00077       INTEGER            ISEED( 4 ), ISEEDY( 4 )
00078       REAL               RESULT( NTESTS )
00079 *     ..
00080 *     .. External Functions ..
00081       REAL               SLAMCH, SLANGE, SLARND
00082       EXTERNAL           SLAMCH, SLANGE, SLARND
00083 *     ..
00084 *     .. External Subroutines ..
00085       EXTERNAL           STRTTF, SGEQRF, SGEQLF, STFSM, STRSM
00086 *     ..
00087 *     .. Intrinsic Functions ..
00088       INTRINSIC          MAX, SQRT
00089 *     ..
00090 *     .. Scalars in Common ..
00091       CHARACTER*32       SRNAMT
00092 *     ..
00093 *     .. Common blocks ..
00094       COMMON             / SRNAMC / SRNAMT
00095 *     ..
00096 *     .. Data statements ..
00097       DATA               ISEEDY / 1988, 1989, 1990, 1991 /
00098       DATA               UPLOS  / 'U', 'L' /
00099       DATA               FORMS  / 'N', 'T' /
00100       DATA               SIDES  / 'L', 'R' /
00101       DATA               TRANSS / 'N', 'T' /
00102       DATA               DIAGS  / 'N', 'U' /
00103 *     ..
00104 *     .. Executable Statements ..
00105 *
00106 *     Initialize constants and the random number seed.
00107 *
00108       NRUN = 0
00109       NFAIL = 0
00110       INFO = 0
00111       DO 10 I = 1, 4
00112          ISEED( I ) = ISEEDY( I )
00113    10 CONTINUE
00114       EPS = SLAMCH( 'Precision' )
00115 *
00116       DO 170 IIM = 1, NN
00117 *
00118          M = NVAL( IIM )
00119 *
00120          DO 160 IIN = 1, NN
00121 *
00122             N = NVAL( IIN )
00123 *
00124             DO 150 IFORM = 1, 2
00125 *
00126                CFORM = FORMS( IFORM )
00127 *
00128                DO 140 IUPLO = 1, 2
00129 *
00130                   UPLO = UPLOS( IUPLO )
00131 *
00132                   DO 130 ISIDE = 1, 2
00133 *
00134                      SIDE = SIDES( ISIDE )
00135 *
00136                      DO 120 ITRANS = 1, 2
00137 *
00138                         TRANS = TRANSS( ITRANS )
00139 *
00140                         DO 110 IDIAG = 1, 2
00141 *
00142                            DIAG = DIAGS( IDIAG )
00143 *
00144                            DO 100 IALPHA = 1, 3
00145 *
00146                               IF ( IALPHA.EQ. 1) THEN
00147                                  ALPHA = ZERO
00148                               ELSE IF ( IALPHA.EQ. 1) THEN
00149                                  ALPHA = ONE
00150                               ELSE
00151                                  ALPHA = SLARND( 2, ISEED )
00152                               END IF
00153 *
00154 *                             All the parameters are set:
00155 *                                CFORM, SIDE, UPLO, TRANS, DIAG, M, N,
00156 *                                and ALPHA
00157 *                             READY TO TEST!
00158 *
00159                               NRUN = NRUN + 1
00160 *
00161                               IF ( ISIDE.EQ.1 ) THEN
00162 *
00163 *                                The case ISIDE.EQ.1 is when SIDE.EQ.'L'
00164 *                                -> A is M-by-M ( B is M-by-N )
00165 *
00166                                  NA = M
00167 *
00168                               ELSE
00169 *
00170 *                                The case ISIDE.EQ.2 is when SIDE.EQ.'R'
00171 *                                -> A is N-by-N ( B is M-by-N )
00172 *
00173                                  NA = N
00174 *
00175                               END IF
00176 *
00177 *                             Generate A our NA--by--NA triangular
00178 *                             matrix. 
00179 *                             Our test is based on forward error so we
00180 *                             do want A to be well conditionned! To get
00181 *                             a well-conditionned triangular matrix, we
00182 *                             take the R factor of the QR/LQ factorization
00183 *                             of a random matrix. 
00184 *
00185                               DO J = 1, NA
00186                                  DO I = 1, NA
00187                                     A( I, J) = SLARND( 2, ISEED )
00188                                  END DO
00189                               END DO
00190 *
00191                               IF ( IUPLO.EQ.1 ) THEN
00192 *
00193 *                                The case IUPLO.EQ.1 is when SIDE.EQ.'U'
00194 *                                -> QR factorization.
00195 *
00196                                  SRNAMT = 'SGEQRF'
00197                                  CALL SGEQRF( NA, NA, A, LDA, TAU,
00198      +                                        S_WORK_SGEQRF, LDA,
00199      +                                        INFO )
00200                               ELSE
00201 *
00202 *                                The case IUPLO.EQ.2 is when SIDE.EQ.'L'
00203 *                                -> QL factorization.
00204 *
00205                                  SRNAMT = 'SGELQF'
00206                                  CALL SGELQF( NA, NA, A, LDA, TAU,
00207      +                                        S_WORK_SGEQRF, LDA,
00208      +                                        INFO )
00209                               END IF
00210 *
00211 *                             Store a copy of A in RFP format (in ARF).
00212 *
00213                               SRNAMT = 'STRTTF'
00214                               CALL STRTTF( CFORM, UPLO, NA, A, LDA, ARF,
00215      +                                     INFO )
00216 *
00217 *                             Generate B1 our M--by--N right-hand side
00218 *                             and store a copy in B2.
00219 *
00220                               DO J = 1, N
00221                                  DO I = 1, M
00222                                     B1( I, J) = SLARND( 2, ISEED )
00223                                     B2( I, J) = B1( I, J)
00224                                  END DO
00225                               END DO
00226 *
00227 *                             Solve op( A ) X = B or X op( A ) = B
00228 *                             with STRSM
00229 *
00230                               SRNAMT = 'STRSM'
00231                               CALL STRSM( SIDE, UPLO, TRANS, DIAG, M, N,
00232      +                               ALPHA, A, LDA, B1, LDA )
00233 *
00234 *                             Solve op( A ) X = B or X op( A ) = B
00235 *                             with STFSM
00236 *
00237                               SRNAMT = 'STFSM'
00238                               CALL STFSM( CFORM, SIDE, UPLO, TRANS,
00239      +                                    DIAG, M, N, ALPHA, ARF, B2,
00240      +                                    LDA )
00241 *
00242 *                             Check that the result agrees.
00243 *
00244                               DO J = 1, N
00245                                  DO I = 1, M
00246                                     B1( I, J) = B2( I, J ) - B1( I, J )
00247                                  END DO
00248                               END DO
00249 *
00250                               RESULT(1) = SLANGE( 'I', M, N, B1, LDA,
00251      +                                            S_WORK_SLANGE )
00252 *
00253                               RESULT(1) = RESULT(1) / SQRT( EPS )
00254      +                                    / MAX ( MAX( M, N), 1 )
00255 *
00256                               IF( RESULT(1).GE.THRESH ) THEN
00257                                  IF( NFAIL.EQ.0 ) THEN
00258                                     WRITE( NOUT, * )
00259                                     WRITE( NOUT, FMT = 9999 )
00260                                  END IF
00261                                  WRITE( NOUT, FMT = 9997 ) 'STFSM', 
00262      +                              CFORM, SIDE, UPLO, TRANS, DIAG, M,
00263      +                              N, RESULT(1)
00264                                  NFAIL = NFAIL + 1
00265                               END IF
00266 *
00267   100                      CONTINUE
00268   110                   CONTINUE
00269   120                CONTINUE
00270   130             CONTINUE
00271   140          CONTINUE
00272   150       CONTINUE
00273   160    CONTINUE
00274   170 CONTINUE
00275 *
00276 *     Print a summary of the results.
00277 *
00278       IF ( NFAIL.EQ.0 ) THEN
00279          WRITE( NOUT, FMT = 9996 ) 'STFSM', NRUN
00280       ELSE
00281          WRITE( NOUT, FMT = 9995 ) 'STFSM', NFAIL, NRUN
00282       END IF
00283 *
00284  9999 FORMAT( 1X, 
00285 ' *** Error(s) or Failure(s) while testing STFSM      +         ***')
00286  9997 FORMAT( 1X, '     Failure in ',A5,', CFORM=''',A1,''',',
00287      + ' SIDE=''',A1,''',',' UPLO=''',A1,''',',' TRANS=''',A1,''',',
00288      + ' DIAG=''',A1,''',',' M=',I3,', N =', I3,', test=',G12.5)
00289  9996 FORMAT( 1X, 'All tests for ',A5,' auxiliary routine passed the ',
00290      +        'threshold (',I5,' tests run)')
00291  9995 FORMAT( 1X, A6, ' auxiliary routine:',I5,' out of ',I5,
00292      +        ' tests failed to pass the threshold')
00293 *
00294       RETURN
00295 *
00296 *     End of SDRVRF3
00297 *
00298       END