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