1       SUBROUTINE SDRVPT( 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       REAL               THRESH
 12 *     ..
 13 *     .. Array Arguments ..
 14       LOGICAL            DOTYPE( * )
 15       INTEGER            NVAL( * )
 16       REAL               A( * ), B( * ), D( * ), E( * ), RWORK( * ),
 17      $                   WORK( * ), X( * ), XACT( * )
 18 *     ..
 19 *
 20 *  Purpose
 21 *  =======
 22 *
 23 *  SDRVPT tests SPTSV and -SVX.
 24 *
 25 *  Arguments
 26 *  =========
 27 *
 28 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
 29 *          The matrix types to be used for testing.  Matrices of type j
 30 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
 31 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
 32 *
 33 *  NN      (input) INTEGER
 34 *          The number of values of N contained in the vector NVAL.
 35 *
 36 *  NVAL    (input) INTEGER array, dimension (NN)
 37 *          The values of the matrix dimension N.
 38 *
 39 *  NRHS    (input) INTEGER
 40 *          The number of right hand side vectors to be generated for
 41 *          each linear system.
 42 *
 43 *  THRESH  (input) REAL
 44 *          The threshold value for the test ratios.  A result is
 45 *          included in the output file if RESULT >= THRESH.  To have
 46 *          every test ratio printed, use THRESH = 0.
 47 *
 48 *  TSTERR  (input) LOGICAL
 49 *          Flag that indicates whether error exits are to be tested.
 50 *
 51 *  A       (workspace) REAL array, dimension (NMAX*2)
 52 *
 53 *  D       (workspace) REAL array, dimension (NMAX*2)
 54 *
 55 *  E       (workspace) REAL array, dimension (NMAX*2)
 56 *
 57 *  B       (workspace) REAL array, dimension (NMAX*NRHS)
 58 *
 59 *  X       (workspace) REAL array, dimension (NMAX*NRHS)
 60 *
 61 *  XACT    (workspace) REAL array, dimension (NMAX*NRHS)
 62 *
 63 *  WORK    (workspace) REAL array, dimension
 64 *                      (NMAX*max(3,NRHS))
 65 *
 66 *  RWORK   (workspace) REAL array, dimension
 67 *                      (max(NMAX,2*NRHS))
 68 *
 69 *  NOUT    (input) INTEGER
 70 *          The unit number for output.
 71 *
 72 *  =====================================================================
 73 *
 74 *     .. Parameters ..
 75       REAL               ONE, ZERO
 76       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+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       REAL               AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
 90 *     ..
 91 *     .. Local Arrays ..
 92       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 93       REAL               RESULT( NTESTS ), Z( 3 )
 94 *     ..
 95 *     .. External Functions ..
 96       INTEGER            ISAMAX
 97       REAL               SASUM, SGET06, SLANST
 98       EXTERNAL           ISAMAX, SASUM, SGET06, SLANST
 99 *     ..
100 *     .. External Subroutines ..
101       EXTERNAL           ALADHD, ALAERH, ALASVM, SCOPY, SERRVX, SGET04,
102      $                   SLACPY, SLAPTM, SLARNV, SLASET, SLATB4, SLATMS,
103      $                   SPTSV, SPTSVX, SPTT01, SPTT02, SPTT05, SPTTRF,
104      $                   SPTTRS, SSCAL
105 *     ..
106 *     .. Intrinsic Functions ..
107       INTRINSIC          ABSMAX
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 ) = 'Single 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 SERRVX( 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 SLATB4.
156 *
157             CALL SLATB4( 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 = 'SLATMS'
167                CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND,
168      $                      ANORM, KL, KU, 'B', A, 2, WORK, INFO )
169 *
170 *              Check the error code from SLATMS.
171 *
172                IF( INFO.NE.0 ) THEN
173                   CALL ALAERH( PATH, 'SLATMS', 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 SLARNV( 2, ISEED, N, D )
199                   CALL SLARNV( 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 = ISAMAX( N, D, 1 )
217                   DMAX = D( IX )
218                   CALL SSCAL( N, ANORM / DMAX, D, 1 )
219                   IF( N.GT.1 )
220      $               CALL SSCAL( 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                      Z( 3 ) = E( IZERO )
266                      E( IZERO-1 ) = ZERO
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 SLARNV( 2, ISEED, N, XACT( IX ) )
279                IX = IX + LDA
280    40       CONTINUE
281 *
282 *           Set the right hand side.
283 *
284             CALL SLAPTM( N, NRHS, ONE, D, E, XACT, LDA, ZERO, B, LDA )
285 *
286             DO 100 IFACT = 12
287                IF( IFACT.EQ.1 ) THEN
288                   FACT = 'F'
289                ELSE
290                   FACT = 'N'
291                END IF
292 *
293 *              Compute the condition number for comparison with
294 *              the value returned by SPTSVX.
295 *
296                IF( ZEROT ) THEN
297                   IF( IFACT.EQ.1 )
298      $               GO TO 100
299                   RCONDC = ZERO
300 *
301                ELSE IF( IFACT.EQ.1 ) THEN
302 *
303 *                 Compute the 1-norm of A.
304 *
305                   ANORM = SLANST( '1', N, D, E )
306 *
307                   CALL SCOPY( N, D, 1, D( N+1 ), 1 )
308                   IF( N.GT.1 )
309      $               CALL SCOPY( N-1, E, 1, E( N+1 ), 1 )
310 *
311 *                 Factor the matrix A.
312 *
313                   CALL SPTTRF( N, D( N+1 ), E( N+1 ), INFO )
314 *
315 *                 Use SPTTRS to solve for one column at a time of
316 *                 inv(A), computing the maximum column sum as we go.
317 *
318                   AINVNM = ZERO
319                   DO 60 I = 1, N
320                      DO 50 J = 1, N
321                         X( J ) = ZERO
322    50                CONTINUE
323                      X( I ) = ONE
324                      CALL SPTTRS( N, 1, D( N+1 ), E( N+1 ), X, LDA,
325      $                            INFO )
326                      AINVNM = MAX( AINVNM, SASUM( N, X, 1 ) )
327    60             CONTINUE
328 *
329 *                 Compute the 1-norm condition number of A.
330 *
331                   IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
332                      RCONDC = ONE
333                   ELSE
334                      RCONDC = ( ONE / ANORM ) / AINVNM
335                   END IF
336                END IF
337 *
338                IF( IFACT.EQ.2 ) THEN
339 *
340 *                 --- Test SPTSV --
341 *
342                   CALL SCOPY( N, D, 1, D( N+1 ), 1 )
343                   IF( N.GT.1 )
344      $               CALL SCOPY( N-1, E, 1, E( N+1 ), 1 )
345                   CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
346 *
347 *                 Factor A as L*D*L' and solve the system A*X = B.
348 *
349                   SRNAMT = 'SPTSV '
350                   CALL SPTSV( N, NRHS, D( N+1 ), E( N+1 ), X, LDA,
351      $                        INFO )
352 *
353 *                 Check error code from SPTSV .
354 *
355                   IF( INFO.NE.IZERO )
356      $               CALL ALAERH( PATH, 'SPTSV ', INFO, IZERO, ' ', N,
357      $                            N, 11, NRHS, IMAT, NFAIL, NERRS,
358      $                            NOUT )
359                   NT = 0
360                   IF( IZERO.EQ.0 ) THEN
361 *
362 *                    Check the factorization by computing the ratio
363 *                       norm(L*D*L' - A) / (n * norm(A) * EPS )
364 *
365                      CALL SPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
366      $                            RESULT1 ) )
367 *
368 *                    Compute the residual in the solution.
369 *
370                      CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
371                      CALL SPTT02( N, NRHS, D, E, X, LDA, WORK, LDA,
372      $                            RESULT2 ) )
373 *
374 *                    Check solution from generated exact solution.
375 *
376                      CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
377      $                            RESULT3 ) )
378                      NT = 3
379                   END IF
380 *
381 *                 Print information about the tests that did not pass
382 *                 the threshold.
383 *
384                   DO 70 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 )'SPTSV ', N, IMAT, K,
389      $                     RESULT( K )
390                         NFAIL = NFAIL + 1
391                      END IF
392    70             CONTINUE
393                   NRUN = NRUN + NT
394                END IF
395 *
396 *              --- Test SPTSVX ---
397 *
398                IF( IFACT.GT.1 ) THEN
399 *
400 *                 Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
401 *
402                   DO 80 I = 1, N - 1
403                      D( N+I ) = ZERO
404                      E( N+I ) = ZERO
405    80             CONTINUE
406                   IF( N.GT.0 )
407      $               D( N+N ) = ZERO
408                END IF
409 *
410                CALL SLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDA )
411 *
412 *              Solve the system and compute the condition number and
413 *              error bounds using SPTSVX.
414 *
415                SRNAMT = 'SPTSVX'
416                CALL SPTSVX( FACT, N, NRHS, D, E, D( N+1 ), E( N+1 ), B,
417      $                      LDA, X, LDA, RCOND, RWORK, RWORK( NRHS+1 ),
418      $                      WORK, INFO )
419 *
420 *              Check the error code from SPTSVX.
421 *
422                IF( INFO.NE.IZERO )
423      $            CALL ALAERH( PATH, 'SPTSVX', INFO, IZERO, FACT, N, N,
424      $                         11, NRHS, IMAT, NFAIL, NERRS, NOUT )
425                IF( IZERO.EQ.0 ) THEN
426                   IF( IFACT.EQ.2 ) THEN
427 *
428 *                    Check the factorization by computing the ratio
429 *                       norm(L*D*L' - A) / (n * norm(A) * EPS )
430 *
431                      K1 = 1
432                      CALL SPTT01( N, D, E, D( N+1 ), E( N+1 ), WORK,
433      $                            RESULT1 ) )
434                   ELSE
435                      K1 = 2
436                   END IF
437 *
438 *                 Compute the residual in the solution.
439 *
440                   CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
441                   CALL SPTT02( N, NRHS, D, E, X, LDA, WORK, LDA,
442      $                         RESULT2 ) )
443 *
444 *                 Check solution from generated exact solution.
445 *
446                   CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
447      $                         RESULT3 ) )
448 *
449 *                 Check error bounds from iterative refinement.
450 *
451                   CALL SPTT05( N, NRHS, D, E, B, LDA, X, LDA, XACT, LDA,
452      $                         RWORK, RWORK( NRHS+1 ), RESULT4 ) )
453                ELSE
454                   K1 = 6
455                END IF
456 *
457 *              Check the reciprocal of the condition number.
458 *
459                RESULT6 ) = SGET06( RCOND, RCONDC )
460 *
461 *              Print information about the tests that did not pass
462 *              the threshold.
463 *
464                DO 90 K = K1, 6
465                   IFRESULT( K ).GE.THRESH ) THEN
466                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
467      $                  CALL ALADHD( NOUT, PATH )
468                      WRITE( NOUT, FMT = 9998 )'SPTSVX', FACT, N, IMAT,
469      $                  K, RESULT( K )
470                      NFAIL = NFAIL + 1
471                   END IF
472    90          CONTINUE
473                NRUN = NRUN + 7 - K1
474   100       CONTINUE
475   110    CONTINUE
476   120 CONTINUE
477 *
478 *     Print a summary of the results.
479 *
480       CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
481 *
482  9999 FORMAT1X, A, ', N =', I5, ', type ', I2, ', test ', I2,
483      $      ', ratio = 'G12.5 )
484  9998 FORMAT1X, A, ', FACT=''', A1, ''', N =', I5, ', type ', I2,
485      $      ', test ', I2, ', ratio = 'G12.5 )
486       RETURN
487 *
488 *     End of SDRVPT
489 *
490       END