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