1       SUBROUTINE ZDRVSP( 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       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 *  ZDRVSP tests the driver routines ZSPSV 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
 58 *                      (NMAX*(NMAX+1)/2)
 59 *
 60 *  AFAC    (workspace) COMPLEX*16 array, dimension
 61 *                      (NMAX*(NMAX+1)/2)
 62 *
 63 *  AINV    (workspace) COMPLEX*16 array, dimension
 64 *                      (NMAX*(NMAX+1)/2)
 65 *
 66 *  B       (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
 67 *
 68 *  X       (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
 69 *
 70 *  XACT    (workspace) COMPLEX*16 array, dimension (NMAX*NRHS)
 71 *
 72 *  WORK    (workspace) COMPLEX*16 array, dimension
 73 *                      (NMAX*max(2,NRHS))
 74 *
 75 *  RWORK   (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ONE, ZERO
 86       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 87       INTEGER            NTYPES, NTESTS
 88       PARAMETER          ( NTYPES = 11, 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       DOUBLE PRECISION   AINVNM, ANORM, CNDNUM, RCOND, RCONDC
100 *     ..
101 *     .. Local Arrays ..
102       CHARACTER          FACTS( NFACT )
103       INTEGER            ISEED( 4 ), ISEEDY( 4 )
104       DOUBLE PRECISION   RESULT( NTESTS )
105 *     ..
106 *     .. External Functions ..
107       DOUBLE PRECISION   DGET06, ZLANSP
108       EXTERNAL           DGET06, ZLANSP
109 *     ..
110 *     .. External Subroutines ..
111       EXTERNAL           ALADHD, ALAERH, ALASVM, XLAENV, ZCOPY, ZERRVX,
112      $                   ZGET04, ZLACPY, ZLARHS, ZLASET, ZLATB4, ZLATMS,
113      $                   ZLATSP, ZPPT05, ZSPSV, ZSPSVX, ZSPT01, ZSPT02,
114      $                   ZSPTRF, ZSPTRI
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          DCMPLXMAXMIN
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 ) = 'Zomplex precision'
137       PATH( 23 ) = 'SP'
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 ZERRVX( 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                IF( IMAT.NE.NTYPES ) THEN
194 *
195 *                 Set up parameters with ZLATB4 and generate a test
196 *                 matrix with ZLATMS.
197 *
198                   CALL ZLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
199      $                         MODE, CNDNUM, DIST )
200 *
201                   SRNAMT = 'ZLATMS'
202                   CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
203      $                         CNDNUM, ANORM, KL, KU, PACKIT, A, LDA,
204      $                         WORK, INFO )
205 *
206 *                 Check error code from ZLATMS.
207 *
208                   IF( INFO.NE.0 ) THEN
209                      CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N,
210      $                            -1-1-1, IMAT, NFAIL, NERRS, NOUT )
211                      GO TO 160
212                   END IF
213 *
214 *                 For types 3-6, zero one or more rows and columns of
215 *                 the matrix to test that INFO is returned correctly.
216 *
217                   IF( ZEROT ) THEN
218                      IF( IMAT.EQ.3 ) THEN
219                         IZERO = 1
220                      ELSE IF( IMAT.EQ.4 ) THEN
221                         IZERO = N
222                      ELSE
223                         IZERO = N / 2 + 1
224                      END IF
225 *
226                      IF( IMAT.LT.6 ) THEN
227 *
228 *                       Set row and column IZERO to zero.
229 *
230                         IF( IUPLO.EQ.1 ) THEN
231                            IOFF = ( IZERO-1 )*IZERO / 2
232                            DO 20 I = 1, IZERO - 1
233                               A( IOFF+I ) = ZERO
234    20                      CONTINUE
235                            IOFF = IOFF + IZERO
236                            DO 30 I = IZERO, N
237                               A( IOFF ) = ZERO
238                               IOFF = IOFF + I
239    30                      CONTINUE
240                         ELSE
241                            IOFF = IZERO
242                            DO 40 I = 1, IZERO - 1
243                               A( IOFF ) = ZERO
244                               IOFF = IOFF + N - I
245    40                      CONTINUE
246                            IOFF = IOFF - IZERO
247                            DO 50 I = IZERO, N
248                               A( IOFF+I ) = ZERO
249    50                      CONTINUE
250                         END IF
251                      ELSE
252                         IF( IUPLO.EQ.1 ) THEN
253 *
254 *                          Set the first IZERO rows and columns to zero.
255 *
256                            IOFF = 0
257                            DO 70 J = 1, N
258                               I2 = MIN( J, IZERO )
259                               DO 60 I = 1, I2
260                                  A( IOFF+I ) = ZERO
261    60                         CONTINUE
262                               IOFF = IOFF + J
263    70                      CONTINUE
264                         ELSE
265 *
266 *                          Set the last IZERO rows and columns to zero.
267 *
268                            IOFF = 0
269                            DO 90 J = 1, N
270                               I1 = MAX( J, IZERO )
271                               DO 80 I = I1, N
272                                  A( IOFF+I ) = ZERO
273    80                         CONTINUE
274                               IOFF = IOFF + N - J
275    90                      CONTINUE
276                         END IF
277                      END IF
278                   ELSE
279                      IZERO = 0
280                   END IF
281                ELSE
282 *
283 *                 Use a special block diagonal matrix to test alternate
284 *                 code for the 2-by-2 blocks.
285 *
286                   CALL ZLATSP( UPLO, N, A, ISEED )
287                END IF
288 *
289                DO 150 IFACT = 1, NFACT
290 *
291 *                 Do first for FACT = 'F', then for other values.
292 *
293                   FACT = FACTS( IFACT )
294 *
295 *                 Compute the condition number for comparison with
296 *                 the value returned by ZSPSVX.
297 *
298                   IF( ZEROT ) THEN
299                      IF( IFACT.EQ.1 )
300      $                  GO TO 150
301                      RCONDC = ZERO
302 *
303                   ELSE IF( IFACT.EQ.1 ) THEN
304 *
305 *                    Compute the 1-norm of A.
306 *
307                      ANORM = ZLANSP( '1', UPLO, N, A, RWORK )
308 *
309 *                    Factor the matrix A.
310 *
311                      CALL ZCOPY( NPP, A, 1, AFAC, 1 )
312                      CALL ZSPTRF( UPLO, N, AFAC, IWORK, INFO )
313 *
314 *                    Compute inv(A) and take its norm.
315 *
316                      CALL ZCOPY( NPP, AFAC, 1, AINV, 1 )
317                      CALL ZSPTRI( UPLO, N, AINV, IWORK, WORK, INFO )
318                      AINVNM = ZLANSP( '1', UPLO, N, AINV, RWORK )
319 *
320 *                    Compute the 1-norm condition number of A.
321 *
322                      IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
323                         RCONDC = ONE
324                      ELSE
325                         RCONDC = ( ONE / ANORM ) / AINVNM
326                      END IF
327                   END IF
328 *
329 *                 Form an exact solution and set the right hand side.
330 *
331                   SRNAMT = 'ZLARHS'
332                   CALL ZLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
333      $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
334      $                         INFO )
335                   XTYPE = 'C'
336 *
337 *                 --- Test ZSPSV  ---
338 *
339                   IF( IFACT.EQ.2 ) THEN
340                      CALL ZCOPY( NPP, A, 1, AFAC, 1 )
341                      CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
342 *
343 *                    Factor the matrix and solve the system using ZSPSV.
344 *
345                      SRNAMT = 'ZSPSV '
346                      CALL ZSPSV( UPLO, N, NRHS, AFAC, IWORK, X, LDA,
347      $                           INFO )
348 *
349 *                    Adjust the expected value of INFO to account for
350 *                    pivoting.
351 *
352                      K = IZERO
353                      IF( K.GT.0 ) THEN
354   100                   CONTINUE
355                         IF( IWORK( K ).LT.0 ) THEN
356                            IF( IWORK( K ).NE.-K ) THEN
357                               K = -IWORK( K )
358                               GO TO 100
359                            END IF
360                         ELSE IF( IWORK( K ).NE.K ) THEN
361                            K = IWORK( K )
362                            GO TO 100
363                         END IF
364                      END IF
365 *
366 *                    Check error code from ZSPSV .
367 *
368                      IF( INFO.NE.K ) THEN
369                         CALL ALAERH( PATH, 'ZSPSV ', INFO, K, UPLO, N,
370      $                               N, -1-1, NRHS, IMAT, NFAIL,
371      $                               NERRS, NOUT )
372                         GO TO 120
373                      ELSE IF( INFO.NE.0 ) THEN
374                         GO TO 120
375                      END IF
376 *
377 *                    Reconstruct matrix from factors and compute
378 *                    residual.
379 *
380                      CALL ZSPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA,
381      $                            RWORK, RESULT1 ) )
382 *
383 *                    Compute residual of the computed solution.
384 *
385                      CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
386                      CALL ZSPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
387      $                            RWORK, RESULT2 ) )
388 *
389 *                    Check solution from generated exact solution.
390 *
391                      CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
392      $                            RESULT3 ) )
393                      NT = 3
394 *
395 *                    Print information about the tests that did not pass
396 *                    the threshold.
397 *
398                      DO 110 K = 1, NT
399                         IFRESULT( K ).GE.THRESH ) THEN
400                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
401      $                        CALL ALADHD( NOUT, PATH )
402                            WRITE( NOUT, FMT = 9999 )'ZSPSV ', UPLO, N,
403      $                        IMAT, K, RESULT( K )
404                            NFAIL = NFAIL + 1
405                         END IF
406   110                CONTINUE
407                      NRUN = NRUN + NT
408   120                CONTINUE
409                   END IF
410 *
411 *                 --- Test ZSPSVX ---
412 *
413                   IF( IFACT.EQ.2 .AND. NPP.GT.0 )
414      $               CALL ZLASET( 'Full', NPP, 1DCMPLX( ZERO ),
415      $                            DCMPLX( ZERO ), AFAC, NPP )
416                   CALL ZLASET( 'Full', N, NRHS, DCMPLX( ZERO ),
417      $                         DCMPLX( ZERO ), X, LDA )
418 *
419 *                 Solve the system and compute the condition number and
420 *                 error bounds using ZSPSVX.
421 *
422                   SRNAMT = 'ZSPSVX'
423                   CALL ZSPSVX( FACT, UPLO, N, NRHS, A, AFAC, IWORK, B,
424      $                         LDA, X, LDA, RCOND, RWORK,
425      $                         RWORK( NRHS+1 ), WORK, RWORK( 2*NRHS+1 ),
426      $                         INFO )
427 *
428 *                 Adjust the expected value of INFO to account for
429 *                 pivoting.
430 *
431                   K = IZERO
432                   IF( K.GT.0 ) THEN
433   130                CONTINUE
434                      IF( IWORK( K ).LT.0 ) THEN
435                         IF( IWORK( K ).NE.-K ) THEN
436                            K = -IWORK( K )
437                            GO TO 130
438                         END IF
439                      ELSE IF( IWORK( K ).NE.K ) THEN
440                         K = IWORK( K )
441                         GO TO 130
442                      END IF
443                   END IF
444 *
445 *                 Check the error code from ZSPSVX.
446 *
447                   IF( INFO.NE.K ) THEN
448                      CALL ALAERH( PATH, 'ZSPSVX', INFO, K, FACT // UPLO,
449      $                            N, N, -1-1, NRHS, IMAT, NFAIL,
450      $                            NERRS, NOUT )
451                      GO TO 150
452                   END IF
453 *
454                   IF( INFO.EQ.0 ) THEN
455                      IF( IFACT.GE.2 ) THEN
456 *
457 *                       Reconstruct matrix from factors and compute
458 *                       residual.
459 *
460                         CALL ZSPT01( UPLO, N, A, AFAC, IWORK, AINV, LDA,
461      $                               RWORK( 2*NRHS+1 ), RESULT1 ) )
462                         K1 = 1
463                      ELSE
464                         K1 = 2
465                      END IF
466 *
467 *                    Compute residual of the computed solution.
468 *
469                      CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
470                      CALL ZSPT02( UPLO, N, NRHS, A, X, LDA, WORK, LDA,
471      $                            RWORK( 2*NRHS+1 ), RESULT2 ) )
472 *
473 *                    Check solution from generated exact solution.
474 *
475                      CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
476      $                            RESULT3 ) )
477 *
478 *                    Check the error bounds from iterative refinement.
479 *
480                      CALL ZPPT05( UPLO, N, NRHS, A, B, LDA, X, LDA,
481      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
482      $                            RESULT4 ) )
483                   ELSE
484                      K1 = 6
485                   END IF
486 *
487 *                 Compare RCOND from ZSPSVX with the computed value
488 *                 in RCONDC.
489 *
490                   RESULT6 ) = DGET06( RCOND, RCONDC )
491 *
492 *                 Print information about the tests that did not pass
493 *                 the threshold.
494 *
495                   DO 140 K = K1, 6
496                      IFRESULT( K ).GE.THRESH ) THEN
497                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
498      $                     CALL ALADHD( NOUT, PATH )
499                         WRITE( NOUT, FMT = 9998 )'ZSPSVX', FACT, UPLO,
500      $                     N, IMAT, K, RESULT( K )
501                         NFAIL = NFAIL + 1
502                      END IF
503   140             CONTINUE
504                   NRUN = NRUN + 7 - K1
505 *
506   150          CONTINUE
507 *
508   160       CONTINUE
509   170    CONTINUE
510   180 CONTINUE
511 *
512 *     Print a summary of the results.
513 *
514       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
515 *
516  9999 FORMAT1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
517      $      ', test ', I2, ', ratio ='G12.5 )
518  9998 FORMAT1X, A, ', FACT=''', A1, ''', UPLO=''', A1, ''', N =', I5,
519      $      ', type ', I2, ', test ', I2, ', ratio ='G12.5 )
520       RETURN
521 *
522 *     End of ZDRVSP
523 *
524       END