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