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