1       SUBROUTINE CCHKPT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
  2      $                   A, D, E, B, X, XACT, WORK, RWORK, NOUT )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       LOGICAL            TSTERR
 10       INTEGER            NN, NNS, NOUT
 11       REAL               THRESH
 12 *     ..
 13 *     .. Array Arguments ..
 14       LOGICAL            DOTYPE( * )
 15       INTEGER            NSVAL( * ), NVAL( * )
 16       REAL               D( * ), RWORK( * )
 17       COMPLEX            A( * ), B( * ), E( * ), WORK( * ), X( * ),
 18      $                   XACT( * )
 19 *     ..
 20 *
 21 *  Purpose
 22 *  =======
 23 *
 24 *  CCHKPT tests CPTTRF, -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 *  A       (workspace) COMPLEX array, dimension (NMAX*2)
 55 *
 56 *  D       (workspace) REAL array, dimension (NMAX*2)
 57 *
 58 *  E       (workspace) COMPLEX array, dimension (NMAX*2)
 59 *
 60 *  B       (workspace) COMPLEX array, dimension (NMAX*NSMAX)
 61 *          where NSMAX is the largest entry in NSVAL.
 62 *
 63 *  X       (workspace) COMPLEX array, dimension (NMAX*NSMAX)
 64 *
 65 *  XACT    (workspace) COMPLEX array, dimension (NMAX*NSMAX)
 66 *
 67 *  WORK    (workspace) COMPLEX array, dimension
 68 *                      (NMAX*max(3,NSMAX))
 69 *
 70 *  RWORK   (workspace) REAL array, dimension
 71 *                      (max(NMAX,2*NSMAX))
 72 *
 73 *  NOUT    (input) INTEGER
 74 *          The unit number for output.
 75 *
 76 *  =====================================================================
 77 *
 78 *     .. Parameters ..
 79       REAL               ONE, ZERO
 80       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 81       INTEGER            NTYPES
 82       PARAMETER          ( NTYPES = 12 )
 83       INTEGER            NTESTS
 84       PARAMETER          ( NTESTS = 7 )
 85 *     ..
 86 *     .. Local Scalars ..
 87       LOGICAL            ZEROT
 88       CHARACTER          DIST, TYPE, UPLO
 89       CHARACTER*3        PATH
 90       INTEGER            I, IA, IMAT, IN, INFO, IRHS, IUPLO, IX, IZERO,
 91      $                   J, K, KL, KU, LDA, MODE, N, NERRS, NFAIL,
 92      $                   NIMAT, NRHS, NRUN
 93       REAL               AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
 94 *     ..
 95 *     .. Local Arrays ..
 96       CHARACTER          UPLOS( 2 )
 97       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 98       REAL               RESULT( NTESTS )
 99       COMPLEX            Z( 3 )
100 *     ..
101 *     .. External Functions ..
102       INTEGER            ISAMAX
103       REAL               CLANHT, SCASUM, SGET06
104       EXTERNAL           ISAMAX, CLANHT, SCASUM, SGET06
105 *     ..
106 *     .. External Subroutines ..
107       EXTERNAL           ALAERH, ALAHD, ALASUM, CCOPY, CERRGT, CGET04,
108      $                   CLACPY, CLAPTM, CLARNV, CLATB4, CLATMS, CPTCON,
109      $                   CPTRFS, CPTT01, CPTT02, CPTT05, CPTTRF, CPTTRS,
110      $                   CSSCAL, SCOPY, SLARNV, SSCAL
111 *     ..
112 *     .. Intrinsic Functions ..
113       INTRINSIC          ABSMAX, REAL
114 *     ..
115 *     .. Scalars in Common ..
116       LOGICAL            LERR, OK
117       CHARACTER*32       SRNAMT
118       INTEGER            INFOT, NUNIT
119 *     ..
120 *     .. Common blocks ..
121       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
122       COMMON             / SRNAMC / SRNAMT
123 *     ..
124 *     .. Data statements ..
125       DATA               ISEEDY / 0001 / , UPLOS / 'U''L' /
126 *     ..
127 *     .. Executable Statements ..
128 *
129       PATH( 11 ) = 'Complex precision'
130       PATH( 23 ) = 'PT'
131       NRUN = 0
132       NFAIL = 0
133       NERRS = 0
134       DO 10 I = 14
135          ISEED( I ) = ISEEDY( I )
136    10 CONTINUE
137 *
138 *     Test the error exits
139 *
140       IF( TSTERR )
141      $   CALL CERRGT( PATH, NOUT )
142       INFOT = 0
143 *
144       DO 120 IN = 1, NN
145 *
146 *        Do for each value of N in NVAL.
147 *
148          N = NVAL( IN )
149          LDA = MAX1, N )
150          NIMAT = NTYPES
151          IF( N.LE.0 )
152      $      NIMAT = 1
153 *
154          DO 110 IMAT = 1, NIMAT
155 *
156 *           Do the tests only if DOTYPE( IMAT ) is true.
157 *
158             IF( N.GT.0 .AND. .NOT.DOTYPE( IMAT ) )
159      $         GO TO 110
160 *
161 *           Set up parameters with CLATB4.
162 *
163             CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
164      $                   COND, DIST )
165 *
166             ZEROT = IMAT.GE.8 .AND. IMAT.LE.10
167             IF( IMAT.LE.6 ) THEN
168 *
169 *              Type 1-6:  generate a Hermitian tridiagonal matrix of
170 *              known condition number in lower triangular band storage.
171 *
172                SRNAMT = 'CLATMS'
173                CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND,
174      $                      ANORM, KL, KU, 'B', A, 2, WORK, INFO )
175 *
176 *              Check the error code from CLATMS.
177 *
178                IF( INFO.NE.0 ) THEN
179                   CALL ALAERH( PATH, 'CLATMS', INFO, 0' ', N, N, KL,
180      $                         KU, -1, IMAT, NFAIL, NERRS, NOUT )
181                   GO TO 110
182                END IF
183                IZERO = 0
184 *
185 *              Copy the matrix to D and E.
186 *
187                IA = 1
188                DO 20 I = 1, N - 1
189                   D( I ) = REAL( A( IA ) )
190                   E( I ) = A( IA+1 )
191                   IA = IA + 2
192    20          CONTINUE
193                IF( N.GT.0 )
194      $            D( N ) = REAL( A( IA ) )
195             ELSE
196 *
197 *              Type 7-12:  generate a diagonally dominant matrix with
198 *              unknown condition number in the vectors D and E.
199 *
200                IF.NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN
201 *
202 *                 Let E be complex, D real, with values from [-1,1].
203 *
204                   CALL SLARNV( 2, ISEED, N, D )
205                   CALL CLARNV( 2, ISEED, N-1, E )
206 *
207 *                 Make the tridiagonal matrix diagonally dominant.
208 *
209                   IF( N.EQ.1 ) THEN
210                      D( 1 ) = ABS( D( 1 ) )
211                   ELSE
212                      D( 1 ) = ABS( D( 1 ) ) + ABS( E( 1 ) )
213                      D( N ) = ABS( D( N ) ) + ABS( E( N-1 ) )
214                      DO 30 I = 2, N - 1
215                         D( I ) = ABS( D( I ) ) + ABS( E( I ) ) +
216      $                           ABS( E( I-1 ) )
217    30                CONTINUE
218                   END IF
219 *
220 *                 Scale D and E so the maximum element is ANORM.
221 *
222                   IX = ISAMAX( N, D, 1 )
223                   DMAX = D( IX )
224                   CALL SSCAL( N, ANORM / DMAX, D, 1 )
225                   CALL CSSCAL( N-1, ANORM / DMAX, E, 1 )
226 *
227                ELSE IF( IZERO.GT.0 ) THEN
228 *
229 *                 Reuse the last matrix by copying back the zeroed out
230 *                 elements.
231 *
232                   IF( IZERO.EQ.1 ) THEN
233                      D( 1 ) = Z( 2 )
234                      IF( N.GT.1 )
235      $                  E( 1 ) = Z( 3 )
236                   ELSE IF( IZERO.EQ.N ) THEN
237                      E( N-1 ) = Z( 1 )
238                      D( N ) = Z( 2 )
239                   ELSE
240                      E( IZERO-1 ) = Z( 1 )
241                      D( IZERO ) = Z( 2 )
242                      E( IZERO ) = Z( 3 )
243                   END IF
244                END IF
245 *
246 *              For types 8-10, set one row and column of the matrix to
247 *              zero.
248 *
249                IZERO = 0
250                IF( IMAT.EQ.8 ) THEN
251                   IZERO = 1
252                   Z( 2 ) = D( 1 )
253                   D( 1 ) = ZERO
254                   IF( N.GT.1 ) THEN
255                      Z( 3 ) = E( 1 )
256                      E( 1 ) = ZERO
257                   END IF
258                ELSE IF( IMAT.EQ.9 ) THEN
259                   IZERO = N
260                   IF( N.GT.1 ) THEN
261                      Z( 1 ) = E( N-1 )
262                      E( N-1 ) = ZERO
263                   END IF
264                   Z( 2 ) = D( N )
265                   D( N ) = ZERO
266                ELSE IF( IMAT.EQ.10 ) THEN
267                   IZERO = ( N+1 ) / 2
268                   IF( IZERO.GT.1 ) THEN
269                      Z( 1 ) = E( IZERO-1 )
270                      Z( 3 ) = E( IZERO )
271                      E( IZERO-1 ) = ZERO
272                      E( IZERO ) = ZERO
273                   END IF
274                   Z( 2 ) = D( IZERO )
275                   D( IZERO ) = ZERO
276                END IF
277             END IF
278 *
279             CALL SCOPY( N, D, 1, D( N+1 ), 1 )
280             IF( N.GT.1 )
281      $         CALL CCOPY( N-1, E, 1, E( N+1 ), 1 )
282 *
283 *+    TEST 1
284 *           Factor A as L*D*L' and compute the ratio
285 *              norm(L*D*L' - A) / (n * norm(A) * EPS )
286 *
287             CALL CPTTRF( N, D( N+1 ), E( N+1 ), INFO )
288 *
289 *           Check error code from CPTTRF.
290 *
291             IF( INFO.NE.IZERO ) THEN
292                CALL ALAERH( PATH, 'CPTTRF', INFO, IZERO, ' ', N, N, -1,
293      $                      -1-1, IMAT, NFAIL, NERRS, NOUT )
294                GO TO 110
295             END IF
296 *
297             IF( INFO.GT.0 ) THEN
298                RCONDC = ZERO
299                GO TO 100
300             END IF
301 *
302             CALL CPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
303      $                   RESULT1 ) )
304 *
305 *           Print the test ratio if greater than or equal to THRESH.
306 *
307             IFRESULT1 ).GE.THRESH ) THEN
308                IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
309      $            CALL ALAHD( NOUT, PATH )
310                WRITE( NOUT, FMT = 9999 )N, IMAT, 1RESULT1 )
311                NFAIL = NFAIL + 1
312             END IF
313             NRUN = NRUN + 1
314 *
315 *           Compute RCONDC = 1 / (norm(A) * norm(inv(A))
316 *
317 *           Compute norm(A).
318 *
319             ANORM = CLANHT( '1', N, D, E )
320 *
321 *           Use CPTTRS to solve for one column at a time of inv(A),
322 *           computing the maximum column sum as we go.
323 *
324             AINVNM = ZERO
325             DO 50 I = 1, N
326                DO 40 J = 1, N
327                   X( J ) = ZERO
328    40          CONTINUE
329                X( I ) = ONE
330                CALL CPTTRS( 'Lower', N, 1, D( N+1 ), E( N+1 ), X, LDA,
331      $                      INFO )
332                AINVNM = MAX( AINVNM, SCASUM( N, X, 1 ) )
333    50       CONTINUE
334             RCONDC = ONE / MAX( ONE, ANORM*AINVNM )
335 *
336             DO 90 IRHS = 1, NNS
337                NRHS = NSVAL( IRHS )
338 *
339 *           Generate NRHS random solution vectors.
340 *
341                IX = 1
342                DO 60 J = 1, NRHS
343                   CALL CLARNV( 2, ISEED, N, XACT( IX ) )
344                   IX = IX + LDA
345    60          CONTINUE
346 *
347                DO 80 IUPLO = 12
348 *
349 *              Do first for UPLO = 'U', then for UPLO = 'L'.
350 *
351                   UPLO = UPLOS( IUPLO )
352 *
353 *              Set the right hand side.
354 *
355                   CALL CLAPTM( UPLO, N, NRHS, ONE, D, E, XACT, LDA,
356      $                         ZERO, B, LDA )
357 *
358 *+    TEST 2
359 *              Solve A*x = b and compute the residual.
360 *
361                   CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
362                   CALL CPTTRS( UPLO, N, NRHS, D( N+1 ), E( N+1 ), X,
363      $                         LDA, INFO )
364 *
365 *              Check error code from CPTTRS.
366 *
367                   IF( INFO.NE.0 )
368      $               CALL ALAERH( PATH, 'CPTTRS', INFO, 0, UPLO, N, N,
369      $                            -1-1, NRHS, IMAT, NFAIL, NERRS,
370      $                            NOUT )
371 *
372                   CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
373                   CALL CPTT02( UPLO, N, NRHS, D, E, X, LDA, WORK, LDA,
374      $                         RESULT2 ) )
375 *
376 *+    TEST 3
377 *              Check solution from generated exact solution.
378 *
379                   CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
380      $                         RESULT3 ) )
381 *
382 *+    TESTS 4, 5, and 6
383 *              Use iterative refinement to improve the solution.
384 *
385                   SRNAMT = 'CPTRFS'
386                   CALL CPTRFS( UPLO, N, NRHS, D, E, D( N+1 ), E( N+1 ),
387      $                         B, LDA, X, LDA, RWORK, RWORK( NRHS+1 ),
388      $                         WORK, RWORK( 2*NRHS+1 ), INFO )
389 *
390 *              Check error code from CPTRFS.
391 *
392                   IF( INFO.NE.0 )
393      $               CALL ALAERH( PATH, 'CPTRFS', INFO, 0, UPLO, N, N,
394      $                            -1-1, NRHS, IMAT, NFAIL, NERRS,
395      $                            NOUT )
396 *
397                   CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
398      $                         RESULT4 ) )
399                   CALL CPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA,
400      $                         RWORK, RWORK( NRHS+1 ), RESULT5 ) )
401 *
402 *              Print information about the tests that did not pass the
403 *              threshold.
404 *
405                   DO 70 K = 26
406                      IFRESULT( K ).GE.THRESH ) THEN
407                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
408      $                     CALL ALAHD( NOUT, PATH )
409                         WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, IMAT,
410      $                     K, RESULT( K )
411                         NFAIL = NFAIL + 1
412                      END IF
413    70             CONTINUE
414                   NRUN = NRUN + 5
415 *
416    80          CONTINUE
417    90       CONTINUE
418 *
419 *+    TEST 7
420 *           Estimate the reciprocal of the condition number of the
421 *           matrix.
422 *
423   100       CONTINUE
424             SRNAMT = 'CPTCON'
425             CALL CPTCON( N, D( N+1 ), E( N+1 ), ANORM, RCOND, RWORK,
426      $                   INFO )
427 *
428 *           Check error code from CPTCON.
429 *
430             IF( INFO.NE.0 )
431      $         CALL ALAERH( PATH, 'CPTCON', INFO, 0' ', N, N, -1-1,
432      $                      -1, IMAT, NFAIL, NERRS, NOUT )
433 *
434             RESULT7 ) = SGET06( RCOND, RCONDC )
435 *
436 *           Print the test ratio if greater than or equal to THRESH.
437 *
438             IFRESULT7 ).GE.THRESH ) THEN
439                IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
440      $            CALL ALAHD( NOUT, PATH )
441                WRITE( NOUT, FMT = 9999 )N, IMAT, 7RESULT7 )
442                NFAIL = NFAIL + 1
443             END IF
444             NRUN = NRUN + 1
445   110    CONTINUE
446   120 CONTINUE
447 *
448 *     Print a summary of the results.
449 *
450       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
451 *
452  9999 FORMAT' N =', I5, ', type ', I2, ', test ', I2, ', ratio = ',
453      $      G12.5 )
454  9998 FORMAT' UPLO = ''', A1, ''', N =', I5, ', NRHS =', I3,
455      $        ', type ', I2, ', test ', I2, ', ratio = 'G12.5 )
456       RETURN
457 *
458 *     End of CCHKPT
459 *
460       END