1       SUBROUTINE SDRVRF3( NOUT, NN, NVAL, THRESH, A, LDA, ARF, B1, B2,
  2      +                    S_WORK_SLANGE, S_WORK_SGEQRF, TAU )
  3 *
  4 *  -- LAPACK test routine (version 3.2.0) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2008
  7 *
  8 *     .. Scalar Arguments ..
  9       INTEGER            LDA, NN, NOUT
 10       REAL               THRESH
 11 *     ..
 12 *     .. Array Arguments ..
 13       INTEGER            NVAL( NN )
 14       REAL               A( LDA, * ), ARF( * ), B1( LDA, * ),
 15      +                   B2( LDA, * ), S_WORK_SGEQRF( * ),
 16      +                   S_WORK_SLANGE( * ), TAU( * )
 17 *     ..
 18 *
 19 *  Purpose
 20 *  =======
 21 *
 22 *  SDRVRF3 tests the LAPACK RFP routines:
 23 *      STFSM
 24 *
 25 *  Arguments
 26 *  =========
 27 *
 28 *  NOUT          (input) INTEGER
 29 *                The unit number for output.
 30 *
 31 *  NN            (input) INTEGER
 32 *                The number of values of N contained in the vector NVAL.
 33 *
 34 *  NVAL          (input) INTEGER array, dimension (NN)
 35 *                The values of the matrix dimension N.
 36 *
 37 *  THRESH        (input) REAL
 38 *                The threshold value for the test ratios.  A result is
 39 *                included in the output file if RESULT >= THRESH.  To have
 40 *                every test ratio printed, use THRESH = 0.
 41 *
 42 *  A             (workspace) REAL array, dimension (LDA,NMAX)
 43 *
 44 *  LDA           (input) INTEGER
 45 *                The leading dimension of the array A.  LDA >= max(1,NMAX).
 46 *
 47 *  ARF           (workspace) REAL array, dimension ((NMAX*(NMAX+1))/2).
 48 *
 49 *  B1            (workspace) REAL array, dimension (LDA,NMAX)
 50 *
 51 *  B2            (workspace) REAL array, dimension (LDA,NMAX)
 52 *
 53 *  S_WORK_SLANGE (workspace) REAL array, dimension (NMAX)
 54 *
 55 *  S_WORK_SGEQRF (workspace) REAL array, dimension (NMAX)
 56 *
 57 *  TAU           (workspace) REAL array, dimension (NMAX)
 58 *
 59 *  =====================================================================
 60 *     ..
 61 *     .. Parameters ..
 62       REAL               ZERO, ONE
 63       PARAMETER          ( ZERO = ( 0.0E+00.0E+0 ) ,
 64      +                     ONE  = ( 1.0E+00.0E+0 ) )
 65       INTEGER            NTESTS
 66       PARAMETER          ( NTESTS = 1 )
 67 *     ..
 68 *     .. Local Scalars ..
 69       CHARACTER          UPLO, CFORM, DIAG, TRANS, SIDE
 70       INTEGER            I, IFORM, IIM, IIN, INFO, IUPLO, J, M, N, NA,
 71      +                   NFAIL, NRUN, ISIDE, IDIAG, IALPHA, ITRANS
 72       REAL               EPS, ALPHA
 73 *     ..
 74 *     .. Local Arrays ..
 75       CHARACTER          UPLOS( 2 ), FORMS( 2 ), TRANSS( 2 ),
 76      +                   DIAGS( 2 ), SIDES( 2 )
 77       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 78       REAL               RESULT( NTESTS )
 79 *     ..
 80 *     .. External Functions ..
 81       REAL               SLAMCH, SLANGE, SLARND
 82       EXTERNAL           SLAMCH, SLANGE, SLARND
 83 *     ..
 84 *     .. External Subroutines ..
 85       EXTERNAL           STRTTF, SGEQRF, SGEQLF, STFSM, STRSM
 86 *     ..
 87 *     .. Intrinsic Functions ..
 88       INTRINSIC          MAXSQRT
 89 *     ..
 90 *     .. Scalars in Common ..
 91       CHARACTER*32       SRNAMT
 92 *     ..
 93 *     .. Common blocks ..
 94       COMMON             / SRNAMC / SRNAMT
 95 *     ..
 96 *     .. Data statements ..
 97       DATA               ISEEDY / 1988198919901991 /
 98       DATA               UPLOS  / 'U''L' /
 99       DATA               FORMS  / 'N''T' /
100       DATA               SIDES  / 'L''R' /
101       DATA               TRANSS / 'N''T' /
102       DATA               DIAGS  / 'N''U' /
103 *     ..
104 *     .. Executable Statements ..
105 *
106 *     Initialize constants and the random number seed.
107 *
108       NRUN = 0
109       NFAIL = 0
110       INFO = 0
111       DO 10 I = 14
112          ISEED( I ) = ISEEDY( I )
113    10 CONTINUE
114       EPS = SLAMCH( 'Precision' )
115 *
116       DO 170 IIM = 1, NN
117 *
118          M = NVAL( IIM )
119 *
120          DO 160 IIN = 1, NN
121 *
122             N = NVAL( IIN )
123 *
124             DO 150 IFORM = 12
125 *
126                CFORM = FORMS( IFORM )
127 *
128                DO 140 IUPLO = 12
129 *
130                   UPLO = UPLOS( IUPLO )
131 *
132                   DO 130 ISIDE = 12
133 *
134                      SIDE = SIDES( ISIDE )
135 *
136                      DO 120 ITRANS = 12
137 *
138                         TRANS = TRANSS( ITRANS )
139 *
140                         DO 110 IDIAG = 12
141 *
142                            DIAG = DIAGS( IDIAG )
143 *
144                            DO 100 IALPHA = 13
145 *
146                               IF ( IALPHA.EQ. 1THEN
147                                  ALPHA = ZERO
148                               ELSE IF ( IALPHA.EQ. 1THEN
149                                  ALPHA = ONE
150                               ELSE
151                                  ALPHA = SLARND( 2, ISEED )
152                               END IF
153 *
154 *                             All the parameters are set:
155 *                                CFORM, SIDE, UPLO, TRANS, DIAG, M, N,
156 *                                and ALPHA
157 *                             READY TO TEST!
158 *
159                               NRUN = NRUN + 1
160 *
161                               IF ( ISIDE.EQ.1 ) THEN
162 *
163 *                                The case ISIDE.EQ.1 is when SIDE.EQ.'L'
164 *                                -> A is M-by-M ( B is M-by-N )
165 *
166                                  NA = M
167 *
168                               ELSE
169 *
170 *                                The case ISIDE.EQ.2 is when SIDE.EQ.'R'
171 *                                -> A is N-by-N ( B is M-by-N )
172 *
173                                  NA = N
174 *
175                               END IF
176 *
177 *                             Generate A our NA--by--NA triangular
178 *                             matrix. 
179 *                             Our test is based on forward error so we
180 *                             do want A to be well conditionned! To get
181 *                             a well-conditionned triangular matrix, we
182 *                             take the R factor of the QR/LQ factorization
183 *                             of a random matrix. 
184 *
185                               DO J = 1, NA
186                                  DO I = 1, NA
187                                     A( I, J) = SLARND( 2, ISEED )
188                                  END DO
189                               END DO
190 *
191                               IF ( IUPLO.EQ.1 ) THEN
192 *
193 *                                The case IUPLO.EQ.1 is when SIDE.EQ.'U'
194 *                                -> QR factorization.
195 *
196                                  SRNAMT = 'SGEQRF'
197                                  CALL SGEQRF( NA, NA, A, LDA, TAU,
198      +                                        S_WORK_SGEQRF, LDA,
199      +                                        INFO )
200                               ELSE
201 *
202 *                                The case IUPLO.EQ.2 is when SIDE.EQ.'L'
203 *                                -> QL factorization.
204 *
205                                  SRNAMT = 'SGELQF'
206                                  CALL SGELQF( NA, NA, A, LDA, TAU,
207      +                                        S_WORK_SGEQRF, LDA,
208      +                                        INFO )
209                               END IF
210 *
211 *                             Store a copy of A in RFP format (in ARF).
212 *
213                               SRNAMT = 'STRTTF'
214                               CALL STRTTF( CFORM, UPLO, NA, A, LDA, ARF,
215      +                                     INFO )
216 *
217 *                             Generate B1 our M--by--N right-hand side
218 *                             and store a copy in B2.
219 *
220                               DO J = 1, N
221                                  DO I = 1, M
222                                     B1( I, J) = SLARND( 2, ISEED )
223                                     B2( I, J) = B1( I, J)
224                                  END DO
225                               END DO
226 *
227 *                             Solve op( A ) X = B or X op( A ) = B
228 *                             with STRSM
229 *
230                               SRNAMT = 'STRSM'
231                               CALL STRSM( SIDE, UPLO, TRANS, DIAG, M, N,
232      +                               ALPHA, A, LDA, B1, LDA )
233 *
234 *                             Solve op( A ) X = B or X op( A ) = B
235 *                             with STFSM
236 *
237                               SRNAMT = 'STFSM'
238                               CALL STFSM( CFORM, SIDE, UPLO, TRANS,
239      +                                    DIAG, M, N, ALPHA, ARF, B2,
240      +                                    LDA )
241 *
242 *                             Check that the result agrees.
243 *
244                               DO J = 1, N
245                                  DO I = 1, M
246                                     B1( I, J) = B2( I, J ) - B1( I, J )
247                                  END DO
248                               END DO
249 *
250                               RESULT(1= SLANGE( 'I', M, N, B1, LDA,
251      +                                            S_WORK_SLANGE )
252 *
253                               RESULT(1= RESULT(1/ SQRT( EPS )
254      +                                    / MAX ( MAX( M, N), 1 )
255 *
256                               IFRESULT(1).GE.THRESH ) THEN
257                                  IF( NFAIL.EQ.0 ) THEN
258                                     WRITE( NOUT, * )
259                                     WRITE( NOUT, FMT = 9999 )
260                                  END IF
261                                  WRITE( NOUT, FMT = 9997 ) 'STFSM'
262      +                              CFORM, SIDE, UPLO, TRANS, DIAG, M,
263      +                              N, RESULT(1)
264                                  NFAIL = NFAIL + 1
265                               END IF
266 *
267   100                      CONTINUE
268   110                   CONTINUE
269   120                CONTINUE
270   130             CONTINUE
271   140          CONTINUE
272   150       CONTINUE
273   160    CONTINUE
274   170 CONTINUE
275 *
276 *     Print a summary of the results.
277 *
278       IF ( NFAIL.EQ.0 ) THEN
279          WRITE( NOUT, FMT = 9996 ) 'STFSM', NRUN
280       ELSE
281          WRITE( NOUT, FMT = 9995 ) 'STFSM', NFAIL, NRUN
282       END IF
283 *
284  9999 FORMAT1X' *** Error(s) or Failure(s) while testing STFSM 
285      +         ***')
286  9997 FORMAT1X'     Failure in ',A5,', CFORM=''',A1,''',',
287      + ' SIDE=''',A1,''',',' UPLO=''',A1,''',',' TRANS=''',A1,''',',
288      + ' DIAG=''',A1,''',',' M=',I3,', N =', I3,', test=',G12.5)
289  9996 FORMAT1X'All tests for ',A5,' auxiliary routine passed the ',
290      +        'threshold (',I5,' tests run)')
291  9995 FORMAT1X, A6, ' auxiliary routine:',I5,' out of ',I5,
292      +        ' tests failed to pass the threshold')
293 *
294       RETURN
295 *
296 *     End of SDRVRF3
297 *
298       END