1       SUBROUTINE SDRVSY( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
  2      $                   A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
  3      $                   NOUT )
  4 *
  5 *  -- LAPACK test routine (version 3.1) --
  6 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  7 *     November 2006
  8 *
  9 *     .. Scalar Arguments ..
 10       LOGICAL            TSTERR
 11       INTEGER            NMAX, NN, NOUT, NRHS
 12       REAL               THRESH
 13 *     ..
 14 *     .. Array Arguments ..
 15       LOGICAL            DOTYPE( * )
 16       INTEGER            IWORK( * ), NVAL( * )
 17       REAL               A( * ), AFAC( * ), AINV( * ), B( * ),
 18      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
 19 *     ..
 20 *
 21 *  Purpose
 22 *  =======
 23 *
 24 *  SDRVSY tests the driver routines SSYSV and -SVX.
 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 *  NN      (input) INTEGER
 35 *          The number of values of N contained in the vector NVAL.
 36 *
 37 *  NVAL    (input) INTEGER array, dimension (NN)
 38 *          The values of the matrix dimension N.
 39 *
 40 *  NRHS    (input) INTEGER
 41 *          The number of right hand side vectors to be generated for
 42 *          each linear system.
 43 *
 44 *  THRESH  (input) REAL
 45 *          The threshold value for the test ratios.  A result is
 46 *          included in the output file if RESULT >= THRESH.  To have
 47 *          every test ratio printed, use THRESH = 0.
 48 *
 49 *  TSTERR  (input) LOGICAL
 50 *          Flag that indicates whether error exits are to be tested.
 51 *
 52 *  NMAX    (input) INTEGER
 53 *          The maximum value permitted for N, used in dimensioning the
 54 *          work arrays.
 55 *
 56 *  A       (workspace) REAL array, dimension (NMAX*NMAX)
 57 *
 58 *  AFAC    (workspace) REAL array, dimension (NMAX*NMAX)
 59 *
 60 *  AINV    (workspace) REAL array, dimension (NMAX*NMAX)
 61 *
 62 *  B       (workspace) REAL array, dimension (NMAX*NRHS)
 63 *
 64 *  X       (workspace) REAL array, dimension (NMAX*NRHS)
 65 *
 66 *  XACT    (workspace) REAL array, dimension (NMAX*NRHS)
 67 *
 68 *  WORK    (workspace) REAL array, dimension
 69 *                      (NMAX*max(2,NRHS))
 70 *
 71 *  RWORK   (workspace) REAL array, dimension (NMAX+2*NRHS)
 72 *
 73 *  IWORK   (workspace) INTEGER array, dimension (2*NMAX)
 74 *
 75 *  NOUT    (input) INTEGER
 76 *          The unit number for output.
 77 *
 78 *  =====================================================================
 79 *
 80 *     .. Parameters ..
 81       REAL               ONE, ZERO
 82       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 83       INTEGER            NTYPES, NTESTS
 84       PARAMETER          ( NTYPES = 10, NTESTS = 6 )
 85       INTEGER            NFACT
 86       PARAMETER          ( NFACT = 2 )
 87 *     ..
 88 *     .. Local Scalars ..
 89       LOGICAL            ZEROT
 90       CHARACTER          DIST, FACT, TYPE, UPLO, XTYPE
 91       CHARACTER*3        PATH
 92       INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
 93      $                   IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
 94      $                   NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT
 95       REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
 96 *     ..
 97 *     .. Local Arrays ..
 98       CHARACTER          FACTS( NFACT ), UPLOS( 2 )
 99       INTEGER            ISEED( 4 ), ISEEDY( 4 )
100       REAL               RESULT( NTESTS )
101 *     ..
102 *     .. External Functions ..
103       REAL               SGET06, SLANSY
104       EXTERNAL           SGET06, SLANSY
105 *     ..
106 *     .. External Subroutines ..
107       EXTERNAL           ALADHD, ALAERH, ALASVM, SERRVX, SGET04, SLACPY,
108      $                   SLARHS, SLASET, SLATB4, SLATMS, SPOT02, SPOT05,
109      $                   SSYSV, SSYSVX, SSYT01, SSYTRF, SSYTRI2, XLAENV
110 *     ..
111 *     .. Scalars in Common ..
112       LOGICAL            LERR, OK
113       CHARACTER*32       SRNAMT
114       INTEGER            INFOT, NUNIT
115 *     ..
116 *     .. Common blocks ..
117       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
118       COMMON             / SRNAMC / SRNAMT
119 *     ..
120 *     .. Intrinsic Functions ..
121       INTRINSIC          MAXMIN
122 *     ..
123 *     .. Data statements ..
124       DATA               ISEEDY / 1988198919901991 /
125       DATA               UPLOS / 'U''L' / , FACTS / 'F''N' /
126 *     ..
127 *     .. Executable Statements ..
128 *
129 *     Initialize constants and the random number seed.
130 *
131       PATH( 11 ) = 'Single precision'
132       PATH( 23 ) = 'SY'
133       NRUN = 0
134       NFAIL = 0
135       NERRS = 0
136       DO 10 I = 14
137          ISEED( I ) = ISEEDY( I )
138    10 CONTINUE
139       LWORK = MAX2*NMAX, NMAX*NRHS )
140 *
141 *     Test the error exits
142 *
143       IF( TSTERR )
144      $   CALL SERRVX( PATH, NOUT )
145       INFOT = 0
146 *
147 *     Set the block size and minimum block size for testing.
148 *
149       NB = 1
150       NBMIN = 2
151       CALL XLAENV( 1, NB )
152       CALL XLAENV( 2, NBMIN )
153 *
154 *     Do for each value of N in NVAL
155 *
156       DO 180 IN = 1, NN
157          N = NVAL( IN )
158          LDA = MAX( N, 1 )
159          XTYPE = 'N'
160          NIMAT = NTYPES
161          IF( N.LE.0 )
162      $      NIMAT = 1
163 *
164          DO 170 IMAT = 1, NIMAT
165 *
166 *           Do the tests only if DOTYPE( IMAT ) is true.
167 *
168             IF.NOT.DOTYPE( IMAT ) )
169      $         GO TO 170
170 *
171 *           Skip types 3, 4, 5, or 6 if the matrix size is too small.
172 *
173             ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
174             IF( ZEROT .AND. N.LT.IMAT-2 )
175      $         GO TO 170
176 *
177 *           Do first for UPLO = 'U', then for UPLO = 'L'
178 *
179             DO 160 IUPLO = 12
180                UPLO = UPLOS( IUPLO )
181 *
182 *              Set up parameters with SLATB4 and generate a test matrix
183 *              with SLATMS.
184 *
185                CALL SLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
186      $                      CNDNUM, DIST )
187 *
188                SRNAMT = 'SLATMS'
189                CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
190      $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
191      $                      INFO )
192 *
193 *              Check error code from SLATMS.
194 *
195                IF( INFO.NE.0 ) THEN
196                   CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, N, -1,
197      $                         -1-1, IMAT, NFAIL, NERRS, NOUT )
198                   GO TO 160
199                END IF
200 *
201 *              For types 3-6, zero one or more rows and columns of the
202 *              matrix to test that INFO is returned correctly.
203 *
204                IF( ZEROT ) THEN
205                   IF( IMAT.EQ.3 ) THEN
206                      IZERO = 1
207                   ELSE IF( IMAT.EQ.4 ) THEN
208                      IZERO = N
209                   ELSE
210                      IZERO = N / 2 + 1
211                   END IF
212 *
213                   IF( IMAT.LT.6 ) THEN
214 *
215 *                    Set row and column IZERO to zero.
216 *
217                      IF( IUPLO.EQ.1 ) THEN
218                         IOFF = ( IZERO-1 )*LDA
219                         DO 20 I = 1, IZERO - 1
220                            A( IOFF+I ) = ZERO
221    20                   CONTINUE
222                         IOFF = IOFF + IZERO
223                         DO 30 I = IZERO, N
224                            A( IOFF ) = ZERO
225                            IOFF = IOFF + LDA
226    30                   CONTINUE
227                      ELSE
228                         IOFF = IZERO
229                         DO 40 I = 1, IZERO - 1
230                            A( IOFF ) = ZERO
231                            IOFF = IOFF + LDA
232    40                   CONTINUE
233                         IOFF = IOFF - IZERO
234                         DO 50 I = IZERO, N
235                            A( IOFF+I ) = ZERO
236    50                   CONTINUE
237                      END IF
238                   ELSE
239                      IOFF = 0
240                      IF( IUPLO.EQ.1 ) THEN
241 *
242 *                       Set the first IZERO rows and columns to zero.
243 *
244                         DO 70 J = 1, N
245                            I2 = MIN( J, IZERO )
246                            DO 60 I = 1, I2
247                               A( IOFF+I ) = ZERO
248    60                      CONTINUE
249                            IOFF = IOFF + LDA
250    70                   CONTINUE
251                      ELSE
252 *
253 *                       Set the last IZERO rows and columns to zero.
254 *
255                         DO 90 J = 1, N
256                            I1 = MAX( J, IZERO )
257                            DO 80 I = I1, N
258                               A( IOFF+I ) = ZERO
259    80                      CONTINUE
260                            IOFF = IOFF + LDA
261    90                   CONTINUE
262                      END IF
263                   END IF
264                ELSE
265                   IZERO = 0
266                END IF
267 *
268                DO 150 IFACT = 1, NFACT
269 *
270 *                 Do first for FACT = 'F', then for other values.
271 *
272                   FACT = FACTS( IFACT )
273 *
274 *                 Compute the condition number for comparison with
275 *                 the value returned by SSYSVX.
276 *
277                   IF( ZEROT ) THEN
278                      IF( IFACT.EQ.1 )
279      $                  GO TO 150
280                      RCONDC = ZERO
281 *
282                   ELSE IF( IFACT.EQ.1 ) THEN
283 *
284 *                    Compute the 1-norm of A.
285 *
286                      ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
287 *
288 *                    Factor the matrix A.
289 *
290                      CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
291                      CALL SSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK,
292      $                            LWORK, INFO )
293 *
294 *                    Compute inv(A) and take its norm.
295 *
296                      CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
297                      LWORK = (N+NB+1)*(NB+3)
298                      CALL SSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
299      $                            LWORK, INFO )
300                      AINVNM = SLANSY( '1', UPLO, N, AINV, LDA, RWORK )
301 *
302 *                    Compute the 1-norm condition number of A.
303 *
304                      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
305                         RCONDC = ONE
306                      ELSE
307                         RCONDC = ( ONE / ANORM ) / AINVNM
308                      END IF
309                   END IF
310 *
311 *                 Form an exact solution and set the right hand side.
312 *
313                   SRNAMT = 'SLARHS'
314                   CALL SLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
315      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
316      $                         INFO )
317                   XTYPE = 'C'
318 *
319 *                 --- Test SSYSV  ---
320 *
321                   IF( IFACT.EQ.2 ) THEN
322                      CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
323                      CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
324 *
325 *                    Factor the matrix and solve the system using SSYSV.
326 *
327                      SRNAMT = 'SSYSV '
328                      CALL SSYSV( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
329      $                           LDA, WORK, LWORK, INFO )
330 *
331 *                    Adjust the expected value of INFO to account for
332 *                    pivoting.
333 *
334                      K = IZERO
335                      IF( K.GT.0 ) THEN
336   100                   CONTINUE
337                         IF( IWORK( K ).LT.0 ) THEN
338                            IF( IWORK( K ).NE.-K ) THEN
339                               K = -IWORK( K )
340                               GO TO 100
341                            END IF
342                         ELSE IF( IWORK( K ).NE.K ) THEN
343                            K = IWORK( K )
344                            GO TO 100
345                         END IF
346                      END IF
347 *
348 *                    Check error code from SSYSV .
349 *
350                      IF( INFO.NE.K ) THEN
351                         CALL ALAERH( PATH, 'SSYSV ', INFO, K, UPLO, N,
352      $                               N, -1-1, NRHS, IMAT, NFAIL,
353      $                               NERRS, NOUT )
354                         GO TO 120
355                      ELSE IF( INFO.NE.0 ) THEN
356                         GO TO 120
357                      END IF
358 *
359 *                    Reconstruct matrix from factors and compute
360 *                    residual.
361 *
362                      CALL SSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
363      $                            AINV, LDA, RWORK, RESULT1 ) )
364 *
365 *                    Compute residual of the computed solution.
366 *
367                      CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
368                      CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
369      $                            LDA, RWORK, RESULT2 ) )
370 *
371 *                    Check solution from generated exact solution.
372 *
373                      CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
374      $                            RESULT3 ) )
375                      NT = 3
376 *
377 *                    Print information about the tests that did not pass
378 *                    the threshold.
379 *
380                      DO 110 K = 1, NT
381                         IFRESULT( K ).GE.THRESH ) THEN
382                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
383      $                        CALL ALADHD( NOUT, PATH )
384                            WRITE( NOUT, FMT = 9999 )'SSYSV ', UPLO, N,
385      $                        IMAT, K, RESULT( K )
386                            NFAIL = NFAIL + 1
387                         END IF
388   110                CONTINUE
389                      NRUN = NRUN + NT
390   120                CONTINUE
391                   END IF
392 *
393 *                 --- Test SSYSVX ---
394 *
395                   IF( IFACT.EQ.2 )
396      $               CALL SLASET( UPLO, N, N, ZERO, ZERO, AFAC, LDA )
397                   CALL SLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
398 *
399 *                 Solve the system and compute the condition number and
400 *                 error bounds using SSYSVX.
401 *
402                   SRNAMT = 'SSYSVX'
403                   CALL SSYSVX( FACT, UPLO, N, NRHS, A, LDA, AFAC, LDA,
404      $                         IWORK, B, LDA, X, LDA, RCOND, RWORK,
405      $                         RWORK( NRHS+1 ), WORK, LWORK,
406      $                         IWORK( N+1 ), INFO )
407 *
408 *                 Adjust the expected value of INFO to account for
409 *                 pivoting.
410 *
411                   K = IZERO
412                   IF( K.GT.0 ) THEN
413   130                CONTINUE
414                      IF( IWORK( K ).LT.0 ) THEN
415                         IF( IWORK( K ).NE.-K ) THEN
416                            K = -IWORK( K )
417                            GO TO 130
418                         END IF
419                      ELSE IF( IWORK( K ).NE.K ) THEN
420                         K = IWORK( K )
421                         GO TO 130
422                      END IF
423                   END IF
424 *
425 *                 Check the error code from SSYSVX.
426 *
427                   IF( INFO.NE.K ) THEN
428                      CALL ALAERH( PATH, 'SSYSVX', INFO, K, FACT // UPLO,
429      $                            N, N, -1-1, NRHS, IMAT, NFAIL,
430      $                            NERRS, NOUT )
431                      GO TO 150
432                   END IF
433 *
434                   IF( INFO.EQ.0 ) THEN
435                      IF( IFACT.GE.2 ) THEN
436 *
437 *                       Reconstruct matrix from factors and compute
438 *                       residual.
439 *
440                         CALL SSYT01( UPLO, N, A, LDA, AFAC, LDA, IWORK,
441      $                               AINV, LDA, RWORK( 2*NRHS+1 ),
442      $                               RESULT1 ) )
443                         K1 = 1
444                      ELSE
445                         K1 = 2
446                      END IF
447 *
448 *                    Compute residual of the computed solution.
449 *
450                      CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
451                      CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
452      $                            LDA, RWORK( 2*NRHS+1 ), RESULT2 ) )
453 *
454 *                    Check solution from generated exact solution.
455 *
456                      CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
457      $                            RESULT3 ) )
458 *
459 *                    Check the error bounds from iterative refinement.
460 *
461                      CALL SPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
462      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
463      $                            RESULT4 ) )
464                   ELSE
465                      K1 = 6
466                   END IF
467 *
468 *                 Compare RCOND from SSYSVX with the computed value
469 *                 in RCONDC.
470 *
471                   RESULT6 ) = SGET06( RCOND, RCONDC )
472 *
473 *                 Print information about the tests that did not pass
474 *                 the threshold.
475 *
476                   DO 140 K = K1, 6
477                      IFRESULT( K ).GE.THRESH ) THEN
478                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
479      $                     CALL ALADHD( NOUT, PATH )
480                         WRITE( NOUT, FMT = 9998 )'SSYSVX', FACT, UPLO,
481      $                     N, IMAT, K, RESULT( K )
482                         NFAIL = NFAIL + 1
483                      END IF
484   140             CONTINUE
485                   NRUN = NRUN + 7 - K1
486 *
487   150          CONTINUE
488 *
489   160       CONTINUE
490   170    CONTINUE
491   180 CONTINUE
492 *
493 *     Print a summary of the results.
494 *
495       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
496 *
497  9999 FORMAT1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
498      $      ', test ', I2, ', ratio ='G12.5 )
499  9998 FORMAT1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N =', I5,
500      $      ', type ', I2, ', test ', I2, ', ratio ='G12.5 )
501       RETURN
502 *
503 *     End of SDRVSY
504 *
505       END