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