1       SUBROUTINE ZDRVPT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D,
  2      $                   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, NOUT, NRHS
 11       DOUBLE PRECISION   THRESH
 12 *     ..
 13 *     .. Array Arguments ..
 14       LOGICAL            DOTYPE( * )
 15       INTEGER            NVAL( * )
 16       DOUBLE PRECISION   D( * ), RWORK( * )
 17       COMPLEX*16         A( * ), B( * ), E( * ), WORK( * ), X( * ),
 18      $                   XACT( * )
 19 *     ..
 20 *
 21 *  Purpose
 22 *  =======
 23 *
 24 *  ZDRVPT tests ZPTSV and -SVX.
 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 *  NRHS    (input) INTEGER
 41 *          The number of right hand side vectors to be generated for
 42 *          each linear system.
 43 *
 44 *  THRESH  (input) DOUBLE PRECISION
 45 *          The threshold value for the test ratios.  A result is
 46 *          included in the output file if RESULT >= THRESH.  To have
 47 *          every test ratio printed, use THRESH = 0.
 48 *
 49 *  TSTERR  (input) LOGICAL
 50 *          Flag that indicates whether error exits are to be tested.
 51 *
 52 *  A       (workspace) COMPLEX*16 array, dimension (NMAX*2)
 53 *
 54 *  D       (workspace) DOUBLE PRECISION array, dimension (NMAX*2)
 55 *
 56 *  E       (workspace) COMPLEX*16 array, dimension (NMAX*2)
 57 *
 58 *  B       (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
 59 *
 60 *  X       (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
 61 *
 62 *  XACT    (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
 63 *
 64 *  WORK    (workspace) COMPLEX*16 array, dimension
 65 *                      (NMAX*max(3,NRHS))
 66 *
 67 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
 68 *
 69 *  NOUT    (input) INTEGER
 70 *          The unit number for output.
 71 *
 72 *  =====================================================================
 73 *
 74 *     .. Parameters ..
 75       DOUBLE PRECISION   ONE, ZERO
 76       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 77       INTEGER            NTYPES
 78       PARAMETER          ( NTYPES = 12 )
 79       INTEGER            NTESTS
 80       PARAMETER          ( NTESTS = 6 )
 81 *     ..
 82 *     .. Local Scalars ..
 83       LOGICAL            ZEROT
 84       CHARACTER          DIST, FACT, TYPE
 85       CHARACTER*3        PATH
 86       INTEGER            I, IA, IFACT, IMAT, IN, INFO, IX, IZERO, J, K,
 87      $                   K1, KL, KU, LDA, MODE, N, NERRS, NFAIL, NIMAT,
 88      $                   NRUN, NT
 89       DOUBLE PRECISION   AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
 90 *     ..
 91 *     .. Local Arrays ..
 92       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 93       DOUBLE PRECISION   RESULT( NTESTS ), Z( 3 )
 94 *     ..
 95 *     .. External Functions ..
 96       INTEGER            IDAMAX
 97       DOUBLE PRECISION   DGET06, DZASUM, ZLANHT
 98       EXTERNAL           IDAMAX, DGET06, DZASUM, ZLANHT
 99 *     ..
100 *     .. External Subroutines ..
101       EXTERNAL           ALADHD, ALAERH, ALASVM, DCOPY, DLARNV, DSCAL,
102      $                   ZCOPY, ZDSCAL, ZERRVX, ZGET04, ZLACPY, ZLAPTM,
103      $                   ZLARNV, ZLASET, ZLATB4, ZLATMS, ZPTSV, ZPTSVX,
104      $                   ZPTT01, ZPTT02, ZPTT05, ZPTTRF, ZPTTRS
105 *     ..
106 *     .. Intrinsic Functions ..
107       INTRINSIC          ABSDCMPLXMAX
108 *     ..
109 *     .. Scalars in Common ..
110       LOGICAL            LERR, OK
111       CHARACTER*32       SRNAMT
112       INTEGER            INFOT, NUNIT
113 *     ..
114 *     .. Common blocks ..
115       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
116       COMMON             / SRNAMC / SRNAMT
117 *     ..
118 *     .. Data statements ..
119       DATA               ISEEDY / 0001 /
120 *     ..
121 *     .. Executable Statements ..
122 *
123       PATH( 11 ) = 'Zomplex precision'
124       PATH( 23 ) = 'PT'
125       NRUN = 0
126       NFAIL = 0
127       NERRS = 0
128       DO 10 I = 14
129          ISEED( I ) = ISEEDY( I )
130    10 CONTINUE
131 *
132 *     Test the error exits
133 *
134       IF( TSTERR )
135      $   CALL ZERRVX( PATH, NOUT )
136       INFOT = 0
137 *
138       DO 120 IN = 1, NN
139 *
140 *        Do for each value of N in NVAL.
141 *
142          N = NVAL( IN )
143          LDA = MAX1, N )
144          NIMAT = NTYPES
145          IF( N.LE.0 )
146      $      NIMAT = 1
147 *
148          DO 110 IMAT = 1, NIMAT
149 *
150 *           Do the tests only if DOTYPE( IMAT ) is true.
151 *
152             IF( N.GT.0 .AND. .NOT.DOTYPE( IMAT ) )
153      $         GO TO 110
154 *
155 *           Set up parameters with ZLATB4.
156 *
157             CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
158      $                   COND, DIST )
159 *
160             ZEROT = IMAT.GE.8 .AND. IMAT.LE.10
161             IF( IMAT.LE.6 ) THEN
162 *
163 *              Type 1-6:  generate a symmetric tridiagonal matrix of
164 *              known condition number in lower triangular band storage.
165 *
166                SRNAMT = 'ZLATMS'
167                CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND,
168      $                      ANORM, KL, KU, 'B', A, 2, WORK, INFO )
169 *
170 *              Check the error code from ZLATMS.
171 *
172                IF( INFO.NE.0 ) THEN
173                   CALL ALAERH( PATH, 'ZLATMS', INFO, 0' ', N, N, KL,
174      $                         KU, -1, IMAT, NFAIL, NERRS, NOUT )
175                   GO TO 110
176                END IF
177                IZERO = 0
178 *
179 *              Copy the matrix to D and E.
180 *
181                IA = 1
182                DO 20 I = 1, N - 1
183                   D( I ) = A( IA )
184                   E( I ) = A( IA+1 )
185                   IA = IA + 2
186    20          CONTINUE
187                IF( N.GT.0 )
188      $            D( N ) = A( IA )
189             ELSE
190 *
191 *              Type 7-12:  generate a diagonally dominant matrix with
192 *              unknown condition number in the vectors D and E.
193 *
194                IF.NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN
195 *
196 *                 Let D and E have values from [-1,1].
197 *
198                   CALL DLARNV( 2, ISEED, N, D )
199                   CALL ZLARNV( 2, ISEED, N-1, E )
200 *
201 *                 Make the tridiagonal matrix diagonally dominant.
202 *
203                   IF( N.EQ.1 ) THEN
204                      D( 1 ) = ABS( D( 1 ) )
205                   ELSE
206                      D( 1 ) = ABS( D( 1 ) ) + ABS( E( 1 ) )
207                      D( N ) = ABS( D( N ) ) + ABS( E( N-1 ) )
208                      DO 30 I = 2, N - 1
209                         D( I ) = ABS( D( I ) ) + ABS( E( I ) ) +
210      $                           ABS( E( I-1 ) )
211    30                CONTINUE
212                   END IF
213 *
214 *                 Scale D and E so the maximum element is ANORM.
215 *
216                   IX = IDAMAX( N, D, 1 )
217                   DMAX = D( IX )
218                   CALL DSCAL( N, ANORM / DMAX, D, 1 )
219                   IF( N.GT.1 )
220      $               CALL ZDSCAL( N-1, ANORM / DMAX, E, 1 )
221 *
222                ELSE IF( IZERO.GT.0 ) THEN
223 *
224 *                 Reuse the last matrix by copying back the zeroed out
225 *                 elements.
226 *
227                   IF( IZERO.EQ.1 ) THEN
228                      D( 1 ) = Z( 2 )
229                      IF( N.GT.1 )
230      $                  E( 1 ) = Z( 3 )
231                   ELSE IF( IZERO.EQ.N ) THEN
232                      E( N-1 ) = Z( 1 )
233                      D( N ) = Z( 2 )
234                   ELSE
235                      E( IZERO-1 ) = Z( 1 )
236                      D( IZERO ) = Z( 2 )
237                      E( IZERO ) = Z( 3 )
238                   END IF
239                END IF
240 *
241 *              For types 8-10, set one row and column of the matrix to
242 *              zero.
243 *
244                IZERO = 0
245                IF( IMAT.EQ.8 ) THEN
246                   IZERO = 1
247                   Z( 2 ) = D( 1 )
248                   D( 1 ) = ZERO
249                   IF( N.GT.1 ) THEN
250                      Z( 3 ) = E( 1 )
251                      E( 1 ) = ZERO
252                   END IF
253                ELSE IF( IMAT.EQ.9 ) THEN
254                   IZERO = N
255                   IF( N.GT.1 ) THEN
256                      Z( 1 ) = E( N-1 )
257                      E( N-1 ) = ZERO
258                   END IF
259                   Z( 2 ) = D( N )
260                   D( N ) = ZERO
261                ELSE IF( IMAT.EQ.10 ) THEN
262                   IZERO = ( N+1 ) / 2
263                   IF( IZERO.GT.1 ) THEN
264                      Z( 1 ) = E( IZERO-1 )
265                      E( IZERO-1 ) = ZERO
266                      Z( 3 ) = E( IZERO )
267                      E( IZERO ) = ZERO
268                   END IF
269                   Z( 2 ) = D( IZERO )
270                   D( IZERO ) = ZERO
271                END IF
272             END IF
273 *
274 *           Generate NRHS random solution vectors.
275 *
276             IX = 1
277             DO 40 J = 1, NRHS
278                CALL ZLARNV( 2, ISEED, N, XACT( IX ) )
279                IX = IX + LDA
280    40       CONTINUE
281 *
282 *           Set the right hand side.
283 *
284             CALL ZLAPTM( 'Lower', N, NRHS, ONE, D, E, XACT, LDA, ZERO,
285      $                   B, LDA )
286 *
287             DO 100 IFACT = 12
288                IF( IFACT.EQ.1 ) THEN
289                   FACT = 'F'
290                ELSE
291                   FACT = 'N'
292                END IF
293 *
294 *              Compute the condition number for comparison with
295 *              the value returned by ZPTSVX.
296 *
297                IF( ZEROT ) THEN
298                   IF( IFACT.EQ.1 )
299      $               GO TO 100
300                   RCONDC = ZERO
301 *
302                ELSE IF( IFACT.EQ.1 ) THEN
303 *
304 *                 Compute the 1-norm of A.
305 *
306                   ANORM = ZLANHT( '1', N, D, E )
307 *
308                   CALL DCOPY( N, D, 1, D( N+1 ), 1 )
309                   IF( N.GT.1 )
310      $               CALL ZCOPY( N-1, E, 1, E( N+1 ), 1 )
311 *
312 *                 Factor the matrix A.
313 *
314                   CALL ZPTTRF( N, D( N+1 ), E( N+1 ), INFO )
315 *
316 *                 Use ZPTTRS to solve for one column at a time of
317 *                 inv(A), computing the maximum column sum as we go.
318 *
319                   AINVNM = ZERO
320                   DO 60 I = 1, N
321                      DO 50 J = 1, N
322                         X( J ) = ZERO
323    50                CONTINUE
324                      X( I ) = ONE
325                      CALL ZPTTRS( 'Lower', N, 1, D( N+1 ), E( N+1 ), X,
326      $                            LDA, INFO )
327                      AINVNM = MAX( AINVNM, DZASUM( N, X, 1 ) )
328    60             CONTINUE
329 *
330 *                 Compute the 1-norm condition number of A.
331 *
332                   IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
333                      RCONDC = ONE
334                   ELSE
335                      RCONDC = ( ONE / ANORM ) / AINVNM
336                   END IF
337                END IF
338 *
339                IF( IFACT.EQ.2 ) THEN
340 *
341 *                 --- Test ZPTSV --
342 *
343                   CALL DCOPY( N, D, 1, D( N+1 ), 1 )
344                   IF( N.GT.1 )
345      $               CALL ZCOPY( N-1, E, 1, E( N+1 ), 1 )
346                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
347 *
348 *                 Factor A as L*D*L' and solve the system A*X = B.
349 *
350                   SRNAMT = 'ZPTSV '
351                   CALL ZPTSV( N, NRHS, D( N+1 ), E( N+1 ), X, LDA,
352      $                        INFO )
353 *
354 *                 Check error code from ZPTSV .
355 *
356                   IF( INFO.NE.IZERO )
357      $               CALL ALAERH( PATH, 'ZPTSV ', INFO, IZERO, ' ', N,
358      $                            N, 11, NRHS, IMAT, NFAIL, NERRS,
359      $                            NOUT )
360                   NT = 0
361                   IF( IZERO.EQ.0 ) THEN
362 *
363 *                    Check the factorization by computing the ratio
364 *                       norm(L*D*L' - A) / (n * norm(A) * EPS )
365 *
366                      CALL ZPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
367      $                            RESULT1 ) )
368 *
369 *                    Compute the residual in the solution.
370 *
371                      CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
372                      CALL ZPTT02( 'Lower', N, NRHS, D, E, X, LDA, WORK,
373      $                            LDA, RESULT2 ) )
374 *
375 *                    Check solution from generated exact solution.
376 *
377                      CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
378      $                            RESULT3 ) )
379                      NT = 3
380                   END IF
381 *
382 *                 Print information about the tests that did not pass
383 *                 the threshold.
384 *
385                   DO 70 K = 1, NT
386                      IFRESULT( K ).GE.THRESH ) THEN
387                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
388      $                     CALL ALADHD( NOUT, PATH )
389                         WRITE( NOUT, FMT = 9999 )'ZPTSV ', N, IMAT, K,
390      $                     RESULT( K )
391                         NFAIL = NFAIL + 1
392                      END IF
393    70             CONTINUE
394                   NRUN = NRUN + NT
395                END IF
396 *
397 *              --- Test ZPTSVX ---
398 *
399                IF( IFACT.GT.1 ) THEN
400 *
401 *                 Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
402 *
403                   DO 80 I = 1, N - 1
404                      D( N+I ) = ZERO
405                      E( N+I ) = ZERO
406    80             CONTINUE
407                   IF( N.GT.0 )
408      $               D( N+N ) = ZERO
409                END IF
410 *
411                CALL ZLASET( 'Full', N, NRHS, DCMPLX( ZERO ),
412      $                      DCMPLX( ZERO ), X, LDA )
413 *
414 *              Solve the system and compute the condition number and
415 *              error bounds using ZPTSVX.
416 *
417                SRNAMT = 'ZPTSVX'
418                CALL ZPTSVX( FACT, N, NRHS, D, E, D( N+1 ), E( N+1 ), B,
419      $                      LDA, X, LDA, RCOND, RWORK, RWORK( NRHS+1 ),
420      $                      WORK, RWORK( 2*NRHS+1 ), INFO )
421 *
422 *              Check the error code from ZPTSVX.
423 *
424                IF( INFO.NE.IZERO )
425      $            CALL ALAERH( PATH, 'ZPTSVX', INFO, IZERO, FACT, N, N,
426      $                         11, NRHS, IMAT, NFAIL, NERRS, NOUT )
427                IF( IZERO.EQ.0 ) THEN
428                   IF( IFACT.EQ.2 ) THEN
429 *
430 *                    Check the factorization by computing the ratio
431 *                       norm(L*D*L' - A) / (n * norm(A) * EPS )
432 *
433                      K1 = 1
434                      CALL ZPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
435      $                            RESULT1 ) )
436                   ELSE
437                      K1 = 2
438                   END IF
439 *
440 *                 Compute the residual in the solution.
441 *
442                   CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
443                   CALL ZPTT02( 'Lower', N, NRHS, D, E, X, LDA, WORK,
444      $                         LDA, RESULT2 ) )
445 *
446 *                 Check solution from generated exact solution.
447 *
448                   CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
449      $                         RESULT3 ) )
450 *
451 *                 Check error bounds from iterative refinement.
452 *
453                   CALL ZPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA,
454      $                         RWORK, RWORK( NRHS+1 ), RESULT4 ) )
455                ELSE
456                   K1 = 6
457                END IF
458 *
459 *              Check the reciprocal of the condition number.
460 *
461                RESULT6 ) = DGET06( RCOND, RCONDC )
462 *
463 *              Print information about the tests that did not pass
464 *              the threshold.
465 *
466                DO 90 K = K1, 6
467                   IFRESULT( K ).GE.THRESH ) THEN
468                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
469      $                  CALL ALADHD( NOUT, PATH )
470                      WRITE( NOUT, FMT = 9998 )'ZPTSVX', FACT, N, IMAT,
471      $                  K, RESULT( K )
472                      NFAIL = NFAIL + 1
473                   END IF
474    90          CONTINUE
475                NRUN = NRUN + 7 - K1
476   100       CONTINUE
477   110    CONTINUE
478   120 CONTINUE
479 *
480 *     Print a summary of the results.
481 *
482       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
483 *
484  9999 FORMAT1X, A, ', N =', I5, ', type ', I2, ', test ', I2,
485      $      ', ratio = 'G12.5 )
486  9998 FORMAT1X, A, ', FACT=''', A1, ''', N =', I5, ', type ', I2,
487      $      ', test ', I2, ', ratio = 'G12.5 )
488       RETURN
489 *
490 *     End of ZDRVPT
491 *
492       END