1       SUBROUTINE DDRVAB( DOTYPE, NM, MVAL, NNS,
  2      $                   NSVAL, THRESH, NMAX, A, AFAC, B,
  3      $                   X, WORK, RWORK, SWORK, IWORK, NOUT )
  4       IMPLICIT NONE
  5 *
  6 *  -- LAPACK test routine (version 3.1) --
  7 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  8 *     June 2010
  9 *
 10 *     .. Scalar Arguments ..
 11       INTEGER            NM, NMAX, NNS, NOUT
 12       DOUBLE PRECISION   THRESH
 13 *     ..
 14 *     .. Array Arguments ..
 15       LOGICAL            DOTYPE( * )
 16       INTEGER            MVAL( * ), NSVAL( * ), IWORK( * )
 17       REAL               SWORK(*)
 18       DOUBLE PRECISION   A( * ), AFAC( * ), B( * ),
 19      $                   RWORK( * ), WORK( * ), X( * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  DDRVAB tests DSGESV
 26 *
 27 *  Arguments
 28 *  =========
 29 *
 30 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
 31 *          The matrix types to be used for testing.  Matrices of type j
 32 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
 33 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
 34 *
 35 *  NM      (input) INTEGER
 36 *          The number of values of M contained in the vector MVAL.
 37 *
 38 *  MVAL    (input) INTEGER array, dimension (NM)
 39 *          The values of the matrix row dimension M.
 40 *
 41 *  NNS     (input) INTEGER
 42 *          The number of values of NRHS contained in the vector NSVAL.
 43 *
 44 *  NSVAL   (input) INTEGER array, dimension (NNS)
 45 *          The values of the number of right hand sides NRHS.
 46 *
 47 *  THRESH  (input) DOUBLE PRECISION
 48 *          The threshold value for the test ratios.  A result is
 49 *          included in the output file if RESULT >= THRESH.  To have
 50 *          every test ratio printed, use THRESH = 0.
 51 *
 52 *  NMAX    (input) INTEGER
 53 *          The maximum value permitted for M or N, used in dimensioning
 54 *          the work arrays.
 55 *
 56 *  A       (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
 57 *
 58 *  AFAC    (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
 59 *
 60 *  B       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
 61 *          where NSMAX is the largest entry in NSVAL.
 62 *
 63 *  X       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
 64 *
 65 *  WORK    (workspace) DOUBLE PRECISION array, dimension
 66 *                      (NMAX*max(3,NSMAX))
 67 *
 68 *  RWORK   (workspace) DOUBLE PRECISION array, dimension
 69 *                      (max(2*NMAX,2*NSMAX+NWORK))
 70 *
 71 *  SWORK   (workspace) REAL array, dimension
 72 *                      (NMAX*(NSMAX+NMAX))
 73 *
 74 *  IWORK   (workspace) INTEGER array, dimension
 75 *                      NMAX
 76 *
 77 *  NOUT    (input) INTEGER
 78 *          The unit number for output.
 79 *
 80 *  =====================================================================
 81 *
 82 *     .. Parameters ..
 83       DOUBLE PRECISION   ONE, ZERO
 84       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 85       INTEGER            NTYPES
 86       PARAMETER          ( NTYPES = 11 )
 87       INTEGER            NTESTS
 88       PARAMETER          ( NTESTS = 1 )
 89 *     ..
 90 *     .. Local Scalars ..
 91       LOGICAL            ZEROT
 92       CHARACTER          DIST, TRANS, TYPE, XTYPE
 93       CHARACTER*3        PATH
 94       INTEGER            I, IM, IMAT, INFO, IOFF, IRHS,
 95      $                   IZERO, KL, KU, LDA, M, MODE, N,
 96      $                   NERRS, NFAIL, NIMAT, NRHS, NRUN
 97       DOUBLE PRECISION   ANORM, CNDNUM
 98 *     ..
 99 *     .. Local Arrays ..
100       INTEGER            ISEED( 4 ), ISEEDY( 4 )
101       DOUBLE PRECISION   RESULT( NTESTS )
102 *     ..
103 *     .. Local Variables ..
104       INTEGER            ITER, KASE
105 *     ..
106 *     .. External Subroutines ..
107       EXTERNAL           ALAERH, ALAHD, DGET08, DLACPY, DLARHS, DLASET,
108      $                   DLATB4, DLATMS
109 *     ..
110 *     .. Intrinsic Functions ..
111       INTRINSIC          DBLEMAXMINSQRT
112 *     ..
113 *     .. Scalars in Common ..
114       LOGICAL            LERR, OK
115       CHARACTER*32       SRNAMT
116       INTEGER            INFOT, NUNIT
117 *     ..
118 *     .. Common blocks ..
119       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
120       COMMON             / SRNAMC / SRNAMT
121 *     ..
122 *     .. Data statements ..
123       DATA               ISEEDY / 2006200720082009 / 
124 *     ..
125 *     .. Executable Statements ..
126 *
127 *     Initialize constants and the random number seed.
128 *
129       KASE = 0
130       PATH( 11 ) = 'Double precision'
131       PATH( 23 ) = 'GE'
132       NRUN = 0
133       NFAIL = 0
134       NERRS = 0
135       DO 10 I = 14
136          ISEED( I ) = ISEEDY( I )
137    10 CONTINUE
138 *
139       INFOT = 0
140 *
141 *     Do for each value of M in MVAL
142 *
143       DO 120 IM = 1, NM
144          M = MVAL( IM )
145          LDA = MAX1, M )
146 *
147          N = M
148          NIMAT = NTYPES
149          IF( M.LE.0 .OR. N.LE.0 )
150      $      NIMAT = 1
151 *
152          DO 100 IMAT = 1, NIMAT
153 *
154 *           Do the tests only if DOTYPE( IMAT ) is true.
155 *
156             IF.NOT.DOTYPE( IMAT ) )
157      $         GO TO 100
158 *
159 *           Skip types 5, 6, or 7 if the matrix size is too small.
160 *
161             ZEROT = IMAT.GE.5 .AND. IMAT.LE.7
162             IF( ZEROT .AND. N.LT.IMAT-4 )
163      $         GO TO 100
164 *
165 *           Set up parameters with DLATB4 and generate a test matrix
166 *           with DLATMS.
167 *
168             CALL DLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE,
169      $                   CNDNUM, DIST )
170 *
171             SRNAMT = 'DLATMS'
172             CALL DLATMS( M, N, DIST, ISEED, TYPE, RWORK, MODE,
173      $                   CNDNUM, ANORM, KL, KU, 'No packing', A, LDA,
174      $                   WORK, INFO )
175 *
176 *           Check error code from DLATMS.
177 *
178             IF( INFO.NE.0 ) THEN
179                CALL ALAERH( PATH, 'DLATMS', INFO, 0' ', M, N, -1,
180      $                      -1-1, IMAT, NFAIL, NERRS, NOUT )
181                GO TO 100
182             END IF
183 *
184 *           For types 5-7, zero one or more columns of the matrix to
185 *           test that INFO is returned correctly.
186 *
187             IF( ZEROT ) THEN
188                IF( IMAT.EQ.5 ) THEN
189                   IZERO = 1
190                ELSE IF( IMAT.EQ.6 ) THEN
191                   IZERO = MIN( M, N )
192                ELSE
193                   IZERO = MIN( M, N ) / 2 + 1
194                END IF
195                IOFF = ( IZERO-1 )*LDA
196                IF( IMAT.LT.7 ) THEN
197                   DO 20 I = 1, M
198                      A( IOFF+I ) = ZERO
199    20             CONTINUE
200                ELSE
201                   CALL DLASET( 'Full', M, N-IZERO+1, ZERO, ZERO,
202      $                         A( IOFF+1 ), LDA )
203                END IF
204             ELSE
205                IZERO = 0
206             END IF
207 *
208             DO 60 IRHS = 1, NNS
209                NRHS = NSVAL( IRHS )
210                XTYPE = 'N'
211                TRANS = 'N'
212 *
213                SRNAMT = 'DLARHS'
214                CALL DLARHS( PATH, XTYPE, ' ', TRANS, N, N, KL,
215      $                      KU, NRHS, A, LDA, X, LDA, B,
216      $                      LDA, ISEED, INFO )
217 *
218                SRNAMT = 'DSGESV'
219 *
220                KASE = KASE + 1
221 *
222                CALL DLACPY( 'Full', M, N, A, LDA, AFAC, LDA )
223 *
224                CALL DSGESV( N, NRHS, A, LDA, IWORK, B, LDA, X, LDA,
225      $                      WORK, SWORK, ITER, INFO)
226 *
227                IF (ITER.LT.0THEN
228                    CALL DLACPY( 'Full', M, N, AFAC, LDA, A, LDA )
229                ENDIF
230 *
231 *              Check error code from DSGESV. This should be the same as 
232 *              the one of DGETRF.
233 *
234                IF( INFO.NE.IZERO ) THEN
235 *
236                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
237      $               CALL ALAHD( NOUT, PATH )
238                   NERRS = NERRS + 1
239 *
240                   IF( INFO.NE.IZERO .AND. IZERO.NE.0 ) THEN
241                      WRITE( NOUT, FMT = 9988 )'DSGESV',INFO,
242      $                         IZERO,M,IMAT
243                   ELSE
244                      WRITE( NOUT, FMT = 9975 )'DSGESV',INFO,
245      $                         M, IMAT
246                   END IF
247                END IF
248 *
249 *              Skip the remaining test if the matrix is singular.
250 *
251                IF( INFO.NE.0 )
252      $            GO TO 100
253 *
254 *              Check the quality of the solution
255 *
256                CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
257 *
258                CALL DGET08( TRANS, N, N, NRHS, A, LDA, X, LDA, WORK,
259      $                      LDA, RWORK, RESULT1 ) )
260 *
261 *              Check if the test passes the tesing.
262 *              Print information about the tests that did not
263 *              pass the testing.
264 *
265 *              If iterative refinement has been used and claimed to 
266 *              be successful (ITER>0), we want
267 *                NORMI(B - A*X)/(NORMI(A)*NORMI(X)*EPS*SRQT(N)) < 1
268 *
269 *              If double precision has been used (ITER<0), we want
270 *                NORMI(B - A*X)/(NORMI(A)*NORMI(X)*EPS) < THRES
271 *              (Cf. the linear solver testing routines)
272 *
273                IF ((THRESH.LE.0.0E+00)
274      $            .OR.((ITER.GE.0).AND.(N.GT.0)
275      $                 .AND.(RESULT(1).GE.SQRT(DBLE(N))))
276      $            .OR.((ITER.LT.0).AND.(RESULT(1).GE.THRESH))) THEN
277 *
278                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) THEN
279                      WRITE( NOUT, FMT = 8999 )'DGE'
280                      WRITE( NOUT, FMT = '( '' Matrix types:'' )' )
281                      WRITE( NOUT, FMT = 8979 )
282                      WRITE( NOUT, FMT = '( '' Test ratios:'' )' )
283                      WRITE( NOUT, FMT = 8960 )1
284                      WRITE( NOUT, FMT = '( '' Messages:'' )' )
285                   END IF
286 *
287                   WRITE( NOUT, FMT = 9998 )TRANS, N, NRHS,
288      $               IMAT, 1RESULT1 )
289                   NFAIL = NFAIL + 1
290                END IF
291                NRUN = NRUN + 1
292    60       CONTINUE
293   100    CONTINUE
294   120 CONTINUE
295 *
296 *     Print a summary of the results.
297 *
298       IF( NFAIL.GT.0 ) THEN
299          WRITE( NOUT, FMT = 9996 )'DSGESV', NFAIL, NRUN
300       ELSE
301          WRITE( NOUT, FMT = 9995 )'DSGESV', NRUN
302       END IF
303       IF( NERRS.GT.0 ) THEN
304          WRITE( NOUT, FMT = 9994 )NERRS
305       END IF
306 *
307  9998 FORMAT' TRANS=''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
308      $      I2, ', test(', I2, ') ='G12.5 )
309  9996 FORMAT1X, A6, ': ', I6, ' out of ', I6,
310      $      ' tests failed to pass the threshold' )
311  9995 FORMAT/1X'All tests for ', A6,
312      $      ' routines passed the threshold (', I6, ' tests run)' )
313  9994 FORMAT6X, I6, ' error messages recorded' )
314 *
315 *     SUBNAM, INFO, INFOE, M, IMAT
316 *
317  9988 FORMAT' *** ', A6, ' returned with INFO =', I5, ' instead of ',
318      $      I5, / ' ==> M =', I5, ', type ',
319      $      I2 )
320 *
321 *     SUBNAM, INFO, M, IMAT
322 *
323  9975 FORMAT' *** Error code from ', A6, '=', I5, ' for M=', I5,
324      $      ', type ', I2 )
325  8999 FORMAT/ 1X, A3, ':  General dense matrices' )
326  8979 FORMAT4X'1. Diagonal'24X'7. Last n/2 columns zero'/ 4X,
327      $      '2. Upper triangular'16X,
328      $      '8. Random, CNDNUM = sqrt(0.1/EPS)'/ 4X,
329      $      '3. Lower triangular'16X'9. Random, CNDNUM = 0.1/EPS',
330      $      / 4X'4. Random, CNDNUM = 2'13X,
331      $      '10. Scaled near underflow'/ 4X'5. First column zero',
332      $      14X'11. Scaled near overflow'/ 4X,
333      $      '6. Last column zero' )
334  8960 FORMAT3X, I2, ': norm_1( B - A * X )  / ',
335      $      '( norm_1(A) * norm_1(X) * EPS * SQRT(N) ) > 1 if ITERREF'
336      $      / 4x'or norm_1( B - A * X )  / ',
337      $      '( norm_1(A) * norm_1(X) * EPS ) > THRES if DGETRF' )
338       RETURN
339 *
340 *     End of DDRVAB
341 *
342       END