1       SUBROUTINE CDRVHP( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
  2      $                   A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
  3      $                   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, NOUT, NRHS
 12       REAL               THRESH
 13 *     ..
 14 *     .. Array Arguments ..
 15       LOGICAL            DOTYPE( * )
 16       INTEGER            IWORK( * ), NVAL( * )
 17       REAL               RWORK( * )
 18       COMPLEX            A( * ), AFAC( * ), AINV( * ), B( * ),
 19      $                   WORK( * ), X( * ), XACT( * )
 20 *     ..
 21 *
 22 *  Purpose
 23 *  =======
 24 *
 25 *  CDRVHP tests the driver routines CHPSV and -SVX.
 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 *  NRHS    (input) INTEGER
 42 *          The number of right hand side vectors to be generated for
 43 *          each linear system.
 44 *
 45 *  THRESH  (input) REAL
 46 *          The threshold value for the test ratios.  A result is
 47 *          included in the output file if RESULT >= THRESH.  To have
 48 *          every test ratio printed, use THRESH = 0.
 49 *
 50 *  TSTERR  (input) LOGICAL
 51 *          Flag that indicates whether error exits are to be tested.
 52 *
 53 *  NMAX    (input) INTEGER
 54 *          The maximum value permitted for N, used in dimensioning the
 55 *          work arrays.
 56 *
 57 *  A       (workspace) COMPLEX array, dimension
 58 *                      (NMAX*(NMAX+1)/2)
 59 *
 60 *  AFAC    (workspace) COMPLEX array, dimension
 61 *                      (NMAX*(NMAX+1)/2)
 62 *
 63 *  AINV    (workspace) COMPLEX array, dimension
 64 *                      (NMAX*(NMAX+1)/2)
 65 *
 66 *  B       (workspace) COMPLEX array, dimension (NMAX*NRHS)
 67 *
 68 *  X       (workspace) COMPLEX array, dimension (NMAX*NRHS)
 69 *
 70 *  XACT    (workspace) COMPLEX array, dimension (NMAX*NRHS)
 71 *
 72 *  WORK    (workspace) COMPLEX array, dimension
 73 *                      (NMAX*max(2,NRHS))
 74 *
 75 *  RWORK   (workspace) REAL array, dimension (NMAX+2*NRHS)
 76 *
 77 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
 78 *
 79 *  NOUT    (input) INTEGER
 80 *          The unit number for output.
 81 *
 82 *  =====================================================================
 83 *
 84 *     .. Parameters ..
 85       REAL               ONE, ZERO
 86       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 87       INTEGER            NTYPES, NTESTS
 88       PARAMETER          ( NTYPES = 10, NTESTS = 6 )
 89       INTEGER            NFACT
 90       PARAMETER          ( NFACT = 2 )
 91 *     ..
 92 *     .. Local Scalars ..
 93       LOGICAL            ZEROT
 94       CHARACTER          DIST, FACT, PACKIT, TYPE, UPLO, XTYPE
 95       CHARACTER*3        PATH
 96       INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
 97      $                   IZERO, J, K, K1, KL, KU, LDA, MODE, N, NB,
 98      $                   NBMIN, NERRS, NFAIL, NIMAT, NPP, NRUN, NT
 99       REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
100 *     ..
101 *     .. Local Arrays ..
102       CHARACTER          FACTS( NFACT )
103       INTEGER            ISEED( 4 ), ISEEDY( 4 )
104       REAL               RESULT( NTESTS )
105 *     ..
106 *     .. External Functions ..
107       REAL               CLANHP, SGET06
108       EXTERNAL           CLANHP, SGET06
109 *     ..
110 *     .. External Subroutines ..
111       EXTERNAL           ALADHD, ALAERH, ALASVM, CCOPY, CERRVX, CGET04,
112      $                   CHPSV, CHPSVX, CHPT01, CHPTRF, CHPTRI, CLACPY,
113      $                   CLAIPD, CLARHS, CLASET, CLATB4, CLATMS, CPPT02,
114      $                   CPPT05, XLAENV
115 *     ..
116 *     .. Scalars in Common ..
117       LOGICAL            LERR, OK
118       CHARACTER*32       SRNAMT
119       INTEGER            INFOT, NUNIT
120 *     ..
121 *     .. Common blocks ..
122       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
123       COMMON             / SRNAMC / SRNAMT
124 *     ..
125 *     .. Intrinsic Functions ..
126       INTRINSIC          CMPLXMAXMIN
127 *     ..
128 *     .. Data statements ..
129       DATA               ISEEDY / 1988198919901991 /
130       DATA               FACTS / 'F''N' /
131 *     ..
132 *     .. Executable Statements ..
133 *
134 *     Initialize constants and the random number seed.
135 *
136       PATH( 11 ) = 'C'
137       PATH( 23 ) = 'HP'
138       NRUN = 0
139       NFAIL = 0
140       NERRS = 0
141       DO 10 I = 14
142          ISEED( I ) = ISEEDY( I )
143    10 CONTINUE
144 *
145 *     Test the error exits
146 *
147       IF( TSTERR )
148      $   CALL CERRVX( PATH, NOUT )
149       INFOT = 0
150 *
151 *     Set the block size and minimum block size for testing.
152 *
153       NB = 1
154       NBMIN = 2
155       CALL XLAENV( 1, NB )
156       CALL XLAENV( 2, NBMIN )
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          NPP = N*( N+1 ) / 2
164          XTYPE = 'N'
165          NIMAT = NTYPES
166          IF( N.LE.0 )
167      $      NIMAT = 1
168 *
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                IF( IUPLO.EQ.1 ) THEN
186                   UPLO = 'U'
187                   PACKIT = 'C'
188                ELSE
189                   UPLO = 'L'
190                   PACKIT = 'R'
191                END IF
192 *
193 *              Set up parameters with CLATB4 and generate a test matrix
194 *              with CLATMS.
195 *
196                CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
197      $                      CNDNUM, DIST )
198 *
199                SRNAMT = 'CLATMS'
200                CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
201      $                      CNDNUM, ANORM, KL, KU, PACKIT, A, LDA, WORK,
202      $                      INFO )
203 *
204 *              Check error code from CLATMS.
205 *
206                IF( INFO.NE.0 ) THEN
207                   CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1,
208      $                         -1-1, IMAT, NFAIL, NERRS, NOUT )
209                   GO TO 160
210                END IF
211 *
212 *              For types 3-6, zero one or more rows and columns of the
213 *              matrix to test that INFO is returned correctly.
214 *
215                IF( ZEROT ) THEN
216                   IF( IMAT.EQ.3 ) THEN
217                      IZERO = 1
218                   ELSE IF( IMAT.EQ.4 ) THEN
219                      IZERO = N
220                   ELSE
221                      IZERO = N / 2 + 1
222                   END IF
223 *
224                   IF( IMAT.LT.6 ) THEN
225 *
226 *                    Set row and column IZERO to zero.
227 *
228                      IF( IUPLO.EQ.1 ) THEN
229                         IOFF = ( IZERO-1 )*IZERO / 2
230                         DO 20 I = 1, IZERO - 1
231                            A( IOFF+I ) = ZERO
232    20                   CONTINUE
233                         IOFF = IOFF + IZERO
234                         DO 30 I = IZERO, N
235                            A( IOFF ) = ZERO
236                            IOFF = IOFF + I
237    30                   CONTINUE
238                      ELSE
239                         IOFF = IZERO
240                         DO 40 I = 1, IZERO - 1
241                            A( IOFF ) = ZERO
242                            IOFF = IOFF + N - I
243    40                   CONTINUE
244                         IOFF = IOFF - IZERO
245                         DO 50 I = IZERO, N
246                            A( IOFF+I ) = ZERO
247    50                   CONTINUE
248                      END IF
249                   ELSE
250                      IOFF = 0
251                      IF( IUPLO.EQ.1 ) THEN
252 *
253 *                       Set the first IZERO rows and columns to zero.
254 *
255                         DO 70 J = 1, N
256                            I2 = MIN( J, IZERO )
257                            DO 60 I = 1, I2
258                               A( IOFF+I ) = ZERO
259    60                      CONTINUE
260                            IOFF = IOFF + J
261    70                   CONTINUE
262                      ELSE
263 *
264 *                       Set the last IZERO rows and columns to zero.
265 *
266                         DO 90 J = 1, N
267                            I1 = MAX( J, IZERO )
268                            DO 80 I = I1, N
269                               A( IOFF+I ) = ZERO
270    80                      CONTINUE
271                            IOFF = IOFF + N - J
272    90                   CONTINUE
273                      END IF
274                   END IF
275                ELSE
276                   IZERO = 0
277                END IF
278 *
279 *              Set the imaginary part of the diagonals.
280 *
281                IF( IUPLO.EQ.1 ) THEN
282                   CALL CLAIPD( N, A, 21 )
283                ELSE
284                   CALL CLAIPD( N, A, N, -1 )
285                END IF
286 *
287                DO 150 IFACT = 1, NFACT
288 *
289 *                 Do first for FACT = 'F', then for other values.
290 *
291                   FACT = FACTS( IFACT )
292 *
293 *                 Compute the condition number for comparison with
294 *                 the value returned by CHPSVX.
295 *
296                   IF( ZEROT ) THEN
297                      IF( IFACT.EQ.1 )
298      $                  GO TO 150
299                      RCONDC = ZERO
300 *
301                   ELSE IF( IFACT.EQ.1 ) THEN
302 *
303 *                    Compute the 1-norm of A.
304 *
305                      ANORM = CLANHP( '1', UPLO, N, A, RWORK )
306 *
307 *                    Factor the matrix A.
308 *
309                      CALL CCOPY( NPP, A, 1, AFAC, 1 )
310                      CALL CHPTRF( UPLO, N, AFAC, IWORK, INFO )
311 *
312 *                    Compute inv(A) and take its norm.
313 *
314                      CALL CCOPY( NPP, AFAC, 1, AINV, 1 )
315                      CALL CHPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
316                      AINVNM = CLANHP( '1', UPLO, N, AINV, RWORK )
317 *
318 *                    Compute the 1-norm condition number of A.
319 *
320                      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
321                         RCONDC = ONE
322                      ELSE
323                         RCONDC = ( ONE / ANORM ) / AINVNM
324                      END IF
325                   END IF
326 *
327 *                 Form an exact solution and set the right hand side.
328 *
329                   SRNAMT = 'CLARHS'
330                   CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
331      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
332      $                         INFO )
333                   XTYPE = 'C'
334 *
335 *                 --- Test CHPSV  ---
336 *
337                   IF( IFACT.EQ.2 ) THEN
338                      CALL CCOPY( NPP, A, 1, AFAC, 1 )
339                      CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
340 *
341 *                    Factor the matrix and solve the system using CHPSV.
342 *
343                      SRNAMT = 'CHPSV '
344                      CALL CHPSV( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
345      $                           INFO )
346 *
347 *                    Adjust the expected value of INFO to account for
348 *                    pivoting.
349 *
350                      K = IZERO
351                      IF( K.GT.0 ) THEN
352   100                   CONTINUE
353                         IF( IWORK( K ).LT.0 ) THEN
354                            IF( IWORK( K ).NE.-K ) THEN
355                               K = -IWORK( K )
356                               GO TO 100
357                            END IF
358                         ELSE IF( IWORK( K ).NE.K ) THEN
359                            K = IWORK( K )
360                            GO TO 100
361                         END IF
362                      END IF
363 *
364 *                    Check error code from CHPSV .
365 *
366                      IF( INFO.NE.K ) THEN
367                         CALL ALAERH( PATH, 'CHPSV ', INFO, K, UPLO, N,
368      $                               N, -1-1, NRHS, IMAT, NFAIL,
369      $                               NERRS, NOUT )
370                         GO TO 120
371                      ELSE IF( INFO.NE.0 ) THEN
372                         GO TO 120
373                      END IF
374 *
375 *                    Reconstruct matrix from factors and compute
376 *                    residual.
377 *
378                      CALL CHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA,
379      $                            RWORK, RESULT1 ) )
380 *
381 *                    Compute residual of the computed solution.
382 *
383                      CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
384                      CALL CPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
385      $                            RWORK, RESULT2 ) )
386 *
387 *                    Check solution from generated exact solution.
388 *
389                      CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
390      $                            RESULT3 ) )
391                      NT = 3
392 *
393 *                    Print information about the tests that did not pass
394 *                    the threshold.
395 *
396                      DO 110 K = 1, NT
397                         IFRESULT( K ).GE.THRESH ) THEN
398                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
399      $                        CALL ALADHD( NOUT, PATH )
400                            WRITE( NOUT, FMT = 9999 )'CHPSV ', UPLO, N,
401      $                        IMAT, K, RESULT( K )
402                            NFAIL = NFAIL + 1
403                         END IF
404   110                CONTINUE
405                      NRUN = NRUN + NT
406   120                CONTINUE
407                   END IF
408 *
409 *                 --- Test CHPSVX ---
410 *
411                   IF( IFACT.EQ.2 .AND. NPP.GT.0 )
412      $               CALL CLASET( 'Full', NPP, 1CMPLX( ZERO ),
413      $                            CMPLX( ZERO ), AFAC, NPP )
414                   CALL CLASET( 'Full', N, NRHS, CMPLX( ZERO ),
415      $                         CMPLX( ZERO ), X, LDA )
416 *
417 *                 Solve the system and compute the condition number and
418 *                 error bounds using CHPSVX.
419 *
420                   SRNAMT = 'CHPSVX'
421                   CALL CHPSVX( FACT, UPLO, N, NRHS, A, AFAC, IWORK, B,
422      $                         LDA, X, LDA, RCOND, RWORK,
423      $                         RWORK( NRHS+1 ), WORK, RWORK( 2*NRHS+1 ),
424      $                         INFO )
425 *
426 *                 Adjust the expected value of INFO to account for
427 *                 pivoting.
428 *
429                   K = IZERO
430                   IF( K.GT.0 ) THEN
431   130                CONTINUE
432                      IF( IWORK( K ).LT.0 ) THEN
433                         IF( IWORK( K ).NE.-K ) THEN
434                            K = -IWORK( K )
435                            GO TO 130
436                         END IF
437                      ELSE IF( IWORK( K ).NE.K ) THEN
438                         K = IWORK( K )
439                         GO TO 130
440                      END IF
441                   END IF
442 *
443 *                 Check the error code from CHPSVX.
444 *
445                   IF( INFO.NE.K ) THEN
446                      CALL ALAERH( PATH, 'CHPSVX', INFO, K, FACT // UPLO,
447      $                            N, N, -1-1, NRHS, IMAT, NFAIL,
448      $                            NERRS, NOUT )
449                      GO TO 150
450                   END IF
451 *
452                   IF( INFO.EQ.0 ) THEN
453                      IF( IFACT.GE.2 ) THEN
454 *
455 *                       Reconstruct matrix from factors and compute
456 *                       residual.
457 *
458                         CALL CHPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA,
459      $                               RWORK( 2*NRHS+1 ), RESULT1 ) )
460                         K1 = 1
461                      ELSE
462                         K1 = 2
463                      END IF
464 *
465 *                    Compute residual of the computed solution.
466 *
467                      CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
468                      CALL CPPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
469      $                            RWORK( 2*NRHS+1 ), RESULT2 ) )
470 *
471 *                    Check solution from generated exact solution.
472 *
473                      CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
474      $                            RESULT3 ) )
475 *
476 *                    Check the error bounds from iterative refinement.
477 *
478                      CALL CPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA,
479      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
480      $                            RESULT4 ) )
481                   ELSE
482                      K1 = 6
483                   END IF
484 *
485 *                 Compare RCOND from CHPSVX with the computed value
486 *                 in RCONDC.
487 *
488                   RESULT6 ) = SGET06( RCOND, RCONDC )
489 *
490 *                 Print information about the tests that did not pass
491 *                 the threshold.
492 *
493                   DO 140 K = K1, 6
494                      IFRESULT( K ).GE.THRESH ) THEN
495                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
496      $                     CALL ALADHD( NOUT, PATH )
497                         WRITE( NOUT, FMT = 9998 )'CHPSVX', FACT, UPLO,
498      $                     N, IMAT, K, RESULT( K )
499                         NFAIL = NFAIL + 1
500                      END IF
501   140             CONTINUE
502                   NRUN = NRUN + 7 - K1
503 *
504   150          CONTINUE
505 *
506   160       CONTINUE
507   170    CONTINUE
508   180 CONTINUE
509 *
510 *     Print a summary of the results.
511 *
512       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
513 *
514  9999 FORMAT1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
515      $      ', test ', I2, ', ratio ='G12.5 )
516  9998 FORMAT1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N =', I5,
517      $      ', type ', I2, ', test ', I2, ', ratio ='G12.5 )
518       RETURN
519 *
520 *     End of CDRVHP
521 *
522       END