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