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