1       SUBROUTINE ZCHKHP( 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       DOUBLE PRECISION   THRESH
 13 *     ..
 14 *     .. Array Arguments ..
 15       LOGICAL            DOTYPE( * )
 16       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
 17       DOUBLE PRECISION   RWORK( * )
 18       COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
 19      $                   WORK( * ), X( * ), XACT( * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  ZCHKHP tests ZHPTRF, -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) DOUBLE PRECISION
 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*16 array, dimension
 60 *                      (NMAX*(NMAX+1)/2)
 61 *
 62 *  AFAC    (workspace) COMPLEX*16 array, dimension
 63 *                      (NMAX*(NMAX+1)/2)
 64 *
 65 *  AINV    (workspace) COMPLEX*16 array, dimension
 66 *                      (NMAX*(NMAX+1)/2)
 67 *
 68 *  B       (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX)
 69 *          where NSMAX is the largest entry in NSVAL.
 70 *
 71 *  X       (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX)
 72 *
 73 *  XACT    (workspace) COMPLEX*16 array, dimension (NMAX*NSMAX)
 74 *
 75 *  WORK    (workspace) COMPLEX*16 array, dimension
 76 *                      (NMAX*max(2,NSMAX))
 77 *
 78 *  RWORK   (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ZERO
 90       PARAMETER          ( ZERO = 0.0D+0 )
 91       INTEGER            NTYPES
 92       PARAMETER          ( NTYPES = 10 )
 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       DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
104 *     ..
105 *     .. Local Arrays ..
106       CHARACTER          UPLOS( 2 )
107       INTEGER            ISEED( 4 ), ISEEDY( 4 )
108       DOUBLE PRECISION   RESULT( NTESTS )
109 *     ..
110 *     .. External Functions ..
111       LOGICAL            LSAME
112       DOUBLE PRECISION   DGET06, ZLANHP
113       EXTERNAL           LSAME, DGET06, ZLANHP
114 *     ..
115 *     .. External Subroutines ..
116       EXTERNAL           ALAERH, ALAHD, ALASUM, ZCOPY, ZERRSY, ZGET04,
117      $                   ZHPCON, ZHPRFS, ZHPT01, ZHPTRF, ZHPTRI, ZHPTRS,
118      $                   ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS, ZPPT02,
119      $                   ZPPT03, ZPPT05
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 ) = 'Zomplex precision'
142       PATH( 23 ) = 'HP'
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 ZERRSY( 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          IZERO = 0
167          DO 160 IMAT = 1, NIMAT
168 *
169 *           Do the tests only if DOTYPE( IMAT ) is true.
170 *
171             IF.NOT.DOTYPE( IMAT ) )
172      $         GO TO 160
173 *
174 *           Skip types 3, 4, 5, or 6 if the matrix size is too small.
175 *
176             ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
177             IF( ZEROT .AND. N.LT.IMAT-2 )
178      $         GO TO 160
179 *
180 *           Do first for UPLO = 'U', then for UPLO = 'L'
181 *
182             DO 150 IUPLO = 12
183                UPLO = UPLOS( IUPLO )
184                IF( LSAME( UPLO, 'U' ) ) THEN
185                   PACKIT = 'C'
186                ELSE
187                   PACKIT = 'R'
188                END IF
189 *
190 *              Set up parameters with ZLATB4 and generate a test matrix
191 *              with ZLATMS.
192 *
193                CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
194      $                      CNDNUM, DIST )
195 *
196                SRNAMT = 'ZLATMS'
197                CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
198      $                      CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK,
199      $                      INFO )
200 *
201 *              Check error code from ZLATMS.
202 *
203                IF( INFO.NE.0 ) THEN
204                   CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1,
205      $                         -1-1, IMAT, NFAIL, NERRS, NOUT )
206                   GO TO 150
207                END IF
208 *
209 *              For types 3-6, zero one or more rows and columns of
210 *              the matrix to test that INFO is returned correctly.
211 *
212                IF( ZEROT ) THEN
213                   IF( IMAT.EQ.3 ) THEN
214                      IZERO = 1
215                   ELSE IF( IMAT.EQ.4 ) THEN
216                      IZERO = N
217                   ELSE
218                      IZERO = N / 2 + 1
219                   END IF
220 *
221                   IF( IMAT.LT.6 ) THEN
222 *
223 *                    Set row and column IZERO to zero.
224 *
225                      IF( IUPLO.EQ.1 ) THEN
226                         IOFF = ( IZERO-1 )*IZERO / 2
227                         DO 20 I = 1, IZERO - 1
228                            A( IOFF+I ) = ZERO
229    20                   CONTINUE
230                         IOFF = IOFF + IZERO
231                         DO 30 I = IZERO, N
232                            A( IOFF ) = ZERO
233                            IOFF = IOFF + I
234    30                   CONTINUE
235                      ELSE
236                         IOFF = IZERO
237                         DO 40 I = 1, IZERO - 1
238                            A( IOFF ) = ZERO
239                            IOFF = IOFF + N - I
240    40                   CONTINUE
241                         IOFF = IOFF - IZERO
242                         DO 50 I = IZERO, N
243                            A( IOFF+I ) = ZERO
244    50                   CONTINUE
245                      END IF
246                   ELSE
247                      IOFF = 0
248                      IF( IUPLO.EQ.1 ) THEN
249 *
250 *                       Set the first IZERO rows and columns to zero.
251 *
252                         DO 70 J = 1, N
253                            I2 = MIN( J, IZERO )
254                            DO 60 I = 1, I2
255                               A( IOFF+I ) = ZERO
256    60                      CONTINUE
257                            IOFF = IOFF + J
258    70                   CONTINUE
259                      ELSE
260 *
261 *                       Set the last IZERO rows and columns to zero.
262 *
263                         DO 90 J = 1, N
264                            I1 = MAX( J, IZERO )
265                            DO 80 I = I1, N
266                               A( IOFF+I ) = ZERO
267    80                      CONTINUE
268                            IOFF = IOFF + N - J
269    90                   CONTINUE
270                      END IF
271                   END IF
272                ELSE
273                   IZERO = 0
274                END IF
275 *
276 *              Set the imaginary part of the diagonals.
277 *
278                IF( IUPLO.EQ.1 ) THEN
279                   CALL ZLAIPD( N, A, 21 )
280                ELSE
281                   CALL ZLAIPD( N, A, N, -1 )
282                END IF
283 *
284 *              Compute the L*D*L' or U*D*U' factorization of the matrix.
285 *
286                NPP = N*( N+1 ) / 2
287                CALL ZCOPY( NPP, A, 1, AFAC, 1 )
288                SRNAMT = 'ZHPTRF'
289                CALL ZHPTRF( UPLO, N, AFAC, IWORK, INFO )
290 *
291 *              Adjust the expected value of INFO to account for
292 *              pivoting.
293 *
294                K = IZERO
295                IF( K.GT.0 ) THEN
296   100             CONTINUE
297                   IF( IWORK( K ).LT.0 ) THEN
298                      IF( IWORK( K ).NE.-K ) THEN
299                         K = -IWORK( K )
300                         GO TO 100
301                      END IF
302                   ELSE IF( IWORK( K ).NE.K ) THEN
303                      K = IWORK( K )
304                      GO TO 100
305                   END IF
306                END IF
307 *
308 *              Check error code from ZHPTRF.
309 *
310                IF( INFO.NE.K )
311      $            CALL ALAERH( PATH, 'ZHPTRF', INFO, K, UPLO, N, N, -1,
312      $                         -1-1, IMAT, NFAIL, NERRS, NOUT )
313                IF( INFO.NE.0 ) THEN
314                   TRFCON = .TRUE.
315                ELSE
316                   TRFCON = .FALSE.
317                END IF
318 *
319 *+    TEST 1
320 *              Reconstruct matrix from factors and compute residual.
321 *
322                CALL ZHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA, RWORK,
323      $                      RESULT1 ) )
324                NT = 1
325 *
326 *+    TEST 2
327 *              Form the inverse and compute the residual.
328 *
329                IF.NOT.TRFCON ) THEN
330                   CALL ZCOPY( NPP, AFAC, 1, AINV, 1 )
331                   SRNAMT = 'ZHPTRI'
332                   CALL ZHPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
333 *
334 *              Check error code from ZHPTRI.
335 *
336                   IF( INFO.NE.0 )
337      $               CALL ALAERH( PATH, 'ZHPTRI', INFO, 0, UPLO, N, N,
338      $                            -1-1-1, IMAT, NFAIL, NERRS, NOUT )
339 *
340                   CALL ZPPT03( UPLO, N, A, AINV, WORK, LDA, RWORK,
341      $                         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, IMAT, K,
353      $                  RESULT( K )
354                      NFAIL = NFAIL + 1
355                   END IF
356   110          CONTINUE
357                NRUN = NRUN + NT
358 *
359 *              Do only the condition estimate if INFO is not 0.
360 *
361                IF( TRFCON ) THEN
362                   RCONDC = ZERO
363                   GO TO 140
364                END IF
365 *
366                DO 130 IRHS = 1, NNS
367                   NRHS = NSVAL( IRHS )
368 *
369 *+    TEST 3
370 *              Solve and compute residual for  A * X = B.
371 *
372                   SRNAMT = 'ZLARHS'
373                   CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
374      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
375      $                         INFO )
376                   XTYPE = 'C'
377                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
378 *
379                   SRNAMT = 'ZHPTRS'
380                   CALL ZHPTRS( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
381      $                         INFO )
382 *
383 *              Check error code from ZHPTRS.
384 *
385                   IF( INFO.NE.0 )
386      $               CALL ALAERH( PATH, 'ZHPTRS', INFO, 0, UPLO, N, N,
387      $                            -1-1, NRHS, IMAT, NFAIL, NERRS,
388      $                            NOUT )
389 *
390                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
391                   CALL ZPPT02( 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 ZGET04( 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 = 'ZHPRFS'
404                   CALL ZHPRFS( 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 ZHPRFS.
409 *
410                   IF( INFO.NE.0 )
411      $               CALL ALAERH( PATH, 'ZHPRFS', INFO, 0, UPLO, N, N,
412      $                            -1-1, NRHS, IMAT, NFAIL, NERRS,
413      $                            NOUT )
414 *
415                   CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
416      $                         RESULT5 ) )
417                   CALL ZPPT05( 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 = ZLANHP( '1', UPLO, N, A, RWORK )
441                SRNAMT = 'ZHPCON'
442                CALL ZHPCON( UPLO, N, AFAC, IWORK, ANORM, RCOND, WORK,
443      $                      INFO )
444 *
445 *              Check error code from ZHPCON.
446 *
447                IF( INFO.NE.0 )
448      $            CALL ALAERH( PATH, 'ZHPCON', INFO, 0, UPLO, N, N, -1,
449      $                         -1-1, IMAT, NFAIL, NERRS, NOUT )
450 *
451                RESULT8 ) = DGET06( 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 ZCHKHP
478 *
479       END