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