1       SUBROUTINE DCHKGT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
  2      $                   A, AF, B, X, XACT, WORK, RWORK, IWORK, 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, NNS, NOUT
 11       DOUBLE PRECISION   THRESH
 12 *     ..
 13 *     .. Array Arguments ..
 14       LOGICAL            DOTYPE( * )
 15       INTEGER            IWORK( * ), NSVAL( * ), NVAL( * )
 16       DOUBLE PRECISION   A( * ), AF( * ), B( * ), RWORK( * ), WORK( * ),
 17      $                   X( * ), XACT( * )
 18 *     ..
 19 *
 20 *  Purpose
 21 *  =======
 22 *
 23 *  DCHKGT tests DGTTRF, -TRS, -RFS, and -CON
 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 *  NNS     (input) INTEGER
 40 *          The number of values of NRHS contained in the vector NSVAL.
 41 *
 42 *  NSVAL   (input) INTEGER array, dimension (NNS)
 43 *          The values of the number of right hand sides NRHS.
 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 *  A       (workspace) DOUBLE PRECISION array, dimension (NMAX*4)
 54 *
 55 *  AF      (workspace) DOUBLE PRECISION array, dimension (NMAX*4)
 56 *
 57 *  B       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
 58 *          where NSMAX is the largest entry in NSVAL.
 59 *
 60 *  X       (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
 61 *
 62 *  XACT    (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
 63 *
 64 *  WORK    (workspace) DOUBLE PRECISION array, dimension
 65 *                      (NMAX*max(3,NSMAX))
 66 *
 67 *  RWORK   (workspace) DOUBLE PRECISION array, dimension
 68 *                      (max(NMAX,2*NSMAX))
 69 *
 70 *  IWORK   (workspace) INTEGER array, dimension (2*NMAX)
 71 *
 72 *  NOUT    (input) INTEGER
 73 *          The unit number for output.
 74 *
 75 *  =====================================================================
 76 *
 77 *     .. Parameters ..
 78       DOUBLE PRECISION   ONE, ZERO
 79       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 80       INTEGER            NTYPES
 81       PARAMETER          ( NTYPES = 12 )
 82       INTEGER            NTESTS
 83       PARAMETER          ( NTESTS = 7 )
 84 *     ..
 85 *     .. Local Scalars ..
 86       LOGICAL            TRFCON, ZEROT
 87       CHARACTER          DIST, NORM, TRANS, TYPE
 88       CHARACTER*3        PATH
 89       INTEGER            I, IMAT, IN, INFO, IRHS, ITRAN, IX, IZERO, J,
 90      $                   K, KL, KOFF, KU, LDA, M, MODE, N, NERRS, NFAIL,
 91      $                   NIMAT, NRHS, NRUN
 92       DOUBLE PRECISION   AINVNM, ANORM, COND, RCOND, RCONDC, RCONDI,
 93      $                   RCONDO
 94 *     ..
 95 *     .. Local Arrays ..
 96       CHARACTER          TRANSS( 3 )
 97       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 98       DOUBLE PRECISION   RESULT( NTESTS ), Z( 3 )
 99 *     ..
100 *     .. External Functions ..
101       DOUBLE PRECISION   DASUM, DGET06, DLANGT
102       EXTERNAL           DASUM, DGET06, DLANGT
103 *     ..
104 *     .. External Subroutines ..
105       EXTERNAL           ALAERH, ALAHD, ALASUM, DCOPY, DERRGE, DGET04,
106      $                   DGTCON, DGTRFS, DGTT01, DGTT02, DGTT05, DGTTRF,
107      $                   DGTTRS, DLACPY, DLAGTM, DLARNV, DLATB4, DLATMS,
108      $                   DSCAL
109 *     ..
110 *     .. Intrinsic Functions ..
111       INTRINSIC          MAX
112 *     ..
113 *     .. Scalars in Common ..
114       LOGICAL            LERR, OK
115       CHARACTER*32       SRNAMT
116       INTEGER            INFOT, NUNIT
117 *     ..
118 *     .. Common blocks ..
119       COMMON             / INFOC / INFOT, NUNIT, OK, LERR
120       COMMON             / SRNAMC / SRNAMT
121 *     ..
122 *     .. Data statements ..
123       DATA               ISEEDY / 0001 / , TRANSS / 'N''T',
124      $                   'C' /
125 *     ..
126 *     .. Executable Statements ..
127 *
128       PATH( 11 ) = 'Double precision'
129       PATH( 23 ) = 'GT'
130       NRUN = 0
131       NFAIL = 0
132       NERRS = 0
133       DO 10 I = 14
134          ISEED( I ) = ISEEDY( I )
135    10 CONTINUE
136 *
137 *     Test the error exits
138 *
139       IF( TSTERR )
140      $   CALL DERRGE( PATH, NOUT )
141       INFOT = 0
142 *
143       DO 110 IN = 1, NN
144 *
145 *        Do for each value of N in NVAL.
146 *
147          N = NVAL( IN )
148          M = MAX( N-10 )
149          LDA = MAX1, N )
150          NIMAT = NTYPES
151          IF( N.LE.0 )
152      $      NIMAT = 1
153 *
154          DO 100 IMAT = 1, NIMAT
155 *
156 *           Do the tests only if DOTYPE( IMAT ) is true.
157 *
158             IF.NOT.DOTYPE( IMAT ) )
159      $         GO TO 100
160 *
161 *           Set up parameters with DLATB4.
162 *
163             CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
164      $                   COND, DIST )
165 *
166             ZEROT = IMAT.GE.8 .AND. IMAT.LE.10
167             IF( IMAT.LE.6 ) THEN
168 *
169 *              Types 1-6:  generate matrices of known condition number.
170 *
171                KOFF = MAX2-KU, 3-MAX1, N ) )
172                SRNAMT = 'DLATMS'
173                CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, COND,
174      $                      ANORM, KL, KU, 'Z', AF( KOFF ), 3, WORK,
175      $                      INFO )
176 *
177 *              Check the error code from DLATMS.
178 *
179                IF( INFO.NE.0 ) THEN
180                   CALL ALAERH( PATH, 'DLATMS', INFO, 0' ', N, N, KL,
181      $                         KU, -1, IMAT, NFAIL, NERRS, NOUT )
182                   GO TO 100
183                END IF
184                IZERO = 0
185 *
186                IF( N.GT.1 ) THEN
187                   CALL DCOPY( N-1, AF( 4 ), 3, A, 1 )
188                   CALL DCOPY( N-1, AF( 3 ), 3, A( N+M+1 ), 1 )
189                END IF
190                CALL DCOPY( N, AF( 2 ), 3, A( M+1 ), 1 )
191             ELSE
192 *
193 *              Types 7-12:  generate tridiagonal matrices with
194 *              unknown condition numbers.
195 *
196                IF.NOT.ZEROT .OR. .NOT.DOTYPE( 7 ) ) THEN
197 *
198 *                 Generate a matrix with elements from [-1,1].
199 *
200                   CALL DLARNV( 2, ISEED, N+2*M, A )
201                   IF( ANORM.NE.ONE )
202      $               CALL DSCAL( N+2*M, ANORM, A, 1 )
203                ELSE IF( IZERO.GT.0 ) THEN
204 *
205 *                 Reuse the last matrix by copying back the zeroed out
206 *                 elements.
207 *
208                   IF( IZERO.EQ.1 ) THEN
209                      A( N ) = Z( 2 )
210                      IF( N.GT.1 )
211      $                  A( 1 ) = Z( 3 )
212                   ELSE IF( IZERO.EQ.N ) THEN
213                      A( 3*N-2 ) = Z( 1 )
214                      A( 2*N-1 ) = Z( 2 )
215                   ELSE
216                      A( 2*N-2+IZERO ) = Z( 1 )
217                      A( N-1+IZERO ) = Z( 2 )
218                      A( IZERO ) = Z( 3 )
219                   END IF
220                END IF
221 *
222 *              If IMAT > 7, set one column of the matrix to 0.
223 *
224                IF.NOT.ZEROT ) THEN
225                   IZERO = 0
226                ELSE IF( IMAT.EQ.8 ) THEN
227                   IZERO = 1
228                   Z( 2 ) = A( N )
229                   A( N ) = ZERO
230                   IF( N.GT.1 ) THEN
231                      Z( 3 ) = A( 1 )
232                      A( 1 ) = ZERO
233                   END IF
234                ELSE IF( IMAT.EQ.9 ) THEN
235                   IZERO = N
236                   Z( 1 ) = A( 3*N-2 )
237                   Z( 2 ) = A( 2*N-1 )
238                   A( 3*N-2 ) = ZERO
239                   A( 2*N-1 ) = ZERO
240                ELSE
241                   IZERO = ( N+1 ) / 2
242                   DO 20 I = IZERO, N - 1
243                      A( 2*N-2+I ) = ZERO
244                      A( N-1+I ) = ZERO
245                      A( I ) = ZERO
246    20             CONTINUE
247                   A( 3*N-2 ) = ZERO
248                   A( 2*N-1 ) = ZERO
249                END IF
250             END IF
251 *
252 *+    TEST 1
253 *           Factor A as L*U and compute the ratio
254 *              norm(L*U - A) / (n * norm(A) * EPS )
255 *
256             CALL DCOPY( N+2*M, A, 1, AF, 1 )
257             SRNAMT = 'DGTTRF'
258             CALL DGTTRF( N, AF, AF( M+1 ), AF( N+M+1 ), AF( N+2*M+1 ),
259      $                   IWORK, INFO )
260 *
261 *           Check error code from DGTTRF.
262 *
263             IF( INFO.NE.IZERO )
264      $         CALL ALAERH( PATH, 'DGTTRF', INFO, IZERO, ' ', N, N, 1,
265      $                      1-1, IMAT, NFAIL, NERRS, NOUT )
266             TRFCON = INFO.NE.0
267 *
268             CALL DGTT01( N, A, A( M+1 ), A( N+M+1 ), AF, AF( M+1 ),
269      $                   AF( N+M+1 ), AF( N+2*M+1 ), IWORK, WORK, LDA,
270      $                   RWORK, RESULT1 ) )
271 *
272 *           Print the test ratio if it is .GE. THRESH.
273 *
274             IFRESULT1 ).GE.THRESH ) THEN
275                IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
276      $            CALL ALAHD( NOUT, PATH )
277                WRITE( NOUT, FMT = 9999 )N, IMAT, 1RESULT1 )
278                NFAIL = NFAIL + 1
279             END IF
280             NRUN = NRUN + 1
281 *
282             DO 50 ITRAN = 12
283                TRANS = TRANSS( ITRAN )
284                IF( ITRAN.EQ.1 ) THEN
285                   NORM = 'O'
286                ELSE
287                   NORM = 'I'
288                END IF
289                ANORM = DLANGT( NORM, N, A, A( M+1 ), A( N+M+1 ) )
290 *
291                IF.NOT.TRFCON ) THEN
292 *
293 *                 Use DGTTRS to solve for one column at a time of inv(A)
294 *                 or inv(A^T), computing the maximum column sum as we
295 *                 go.
296 *
297                   AINVNM = ZERO
298                   DO 40 I = 1, N
299                      DO 30 J = 1, N
300                         X( J ) = ZERO
301    30                CONTINUE
302                      X( I ) = ONE
303                      CALL DGTTRS( TRANS, N, 1, AF, AF( M+1 ),
304      $                            AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
305      $                            LDA, INFO )
306                      AINVNM = MAX( AINVNM, DASUM( N, X, 1 ) )
307    40             CONTINUE
308 *
309 *                 Compute RCONDC = 1 / (norm(A) * norm(inv(A))
310 *
311                   IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
312                      RCONDC = ONE
313                   ELSE
314                      RCONDC = ( ONE / ANORM ) / AINVNM
315                   END IF
316                   IF( ITRAN.EQ.1 ) THEN
317                      RCONDO = RCONDC
318                   ELSE
319                      RCONDI = RCONDC
320                   END IF
321                ELSE
322                   RCONDC = ZERO
323                END IF
324 *
325 *+    TEST 7
326 *              Estimate the reciprocal of the condition number of the
327 *              matrix.
328 *
329                SRNAMT = 'DGTCON'
330                CALL DGTCON( NORM, N, AF, AF( M+1 ), AF( N+M+1 ),
331      $                      AF( N+2*M+1 ), IWORK, ANORM, RCOND, WORK,
332      $                      IWORK( N+1 ), INFO )
333 *
334 *              Check error code from DGTCON.
335 *
336                IF( INFO.NE.0 )
337      $            CALL ALAERH( PATH, 'DGTCON', INFO, 0, NORM, N, N, -1,
338      $                         -1-1, IMAT, NFAIL, NERRS, NOUT )
339 *
340                RESULT7 ) = DGET06( RCOND, RCONDC )
341 *
342 *              Print the test ratio if it is .GE. THRESH.
343 *
344                IFRESULT7 ).GE.THRESH ) THEN
345                   IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
346      $               CALL ALAHD( NOUT, PATH )
347                   WRITE( NOUT, FMT = 9997 )NORM, N, IMAT, 7,
348      $               RESULT7 )
349                   NFAIL = NFAIL + 1
350                END IF
351                NRUN = NRUN + 1
352    50       CONTINUE
353 *
354 *           Skip the remaining tests if the matrix is singular.
355 *
356             IF( TRFCON )
357      $         GO TO 100
358 *
359             DO 90 IRHS = 1, NNS
360                NRHS = NSVAL( IRHS )
361 *
362 *              Generate NRHS random solution vectors.
363 *
364                IX = 1
365                DO 60 J = 1, NRHS
366                   CALL DLARNV( 2, ISEED, N, XACT( IX ) )
367                   IX = IX + LDA
368    60          CONTINUE
369 *
370                DO 80 ITRAN = 13
371                   TRANS = TRANSS( ITRAN )
372                   IF( ITRAN.EQ.1 ) THEN
373                      RCONDC = RCONDO
374                   ELSE
375                      RCONDC = RCONDI
376                   END IF
377 *
378 *                 Set the right hand side.
379 *
380                   CALL DLAGTM( TRANS, N, NRHS, ONE, A, A( M+1 ),
381      $                         A( N+M+1 ), XACT, LDA, ZERO, B, LDA )
382 *
383 *+    TEST 2
384 *                 Solve op(A) * X = B and compute the residual.
385 *
386                   CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
387                   SRNAMT = 'DGTTRS'
388                   CALL DGTTRS( TRANS, N, NRHS, AF, AF( M+1 ),
389      $                         AF( N+M+1 ), AF( N+2*M+1 ), IWORK, X,
390      $                         LDA, INFO )
391 *
392 *                 Check error code from DGTTRS.
393 *
394                   IF( INFO.NE.0 )
395      $               CALL ALAERH( PATH, 'DGTTRS', INFO, 0, TRANS, N, N,
396      $                            -1-1, NRHS, IMAT, NFAIL, NERRS,
397      $                            NOUT )
398 *
399                   CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
400                   CALL DGTT02( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
401      $                         X, LDA, WORK, LDA, RWORK, RESULT2 ) )
402 *
403 *+    TEST 3
404 *                 Check solution from generated exact solution.
405 *
406                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
407      $                         RESULT3 ) )
408 *
409 *+    TESTS 4, 5, and 6
410 *                 Use iterative refinement to improve the solution.
411 *
412                   SRNAMT = 'DGTRFS'
413                   CALL DGTRFS( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
414      $                         AF, AF( M+1 ), AF( N+M+1 ),
415      $                         AF( N+2*M+1 ), IWORK, B, LDA, X, LDA,
416      $                         RWORK, RWORK( NRHS+1 ), WORK,
417      $                         IWORK( N+1 ), INFO )
418 *
419 *                 Check error code from DGTRFS.
420 *
421                   IF( INFO.NE.0 )
422      $               CALL ALAERH( PATH, 'DGTRFS', INFO, 0, TRANS, N, N,
423      $                            -1-1, NRHS, IMAT, NFAIL, NERRS,
424      $                            NOUT )
425 *
426                   CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
427      $                         RESULT4 ) )
428                   CALL DGTT05( TRANS, N, NRHS, A, A( M+1 ), A( N+M+1 ),
429      $                         B, LDA, X, LDA, XACT, LDA, RWORK,
430      $                         RWORK( NRHS+1 ), RESULT5 ) )
431 *
432 *                 Print information about the tests that did not pass
433 *                 the threshold.
434 *
435                   DO 70 K = 26
436                      IFRESULT( K ).GE.THRESH ) THEN
437                         IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
438      $                     CALL ALAHD( NOUT, PATH )
439                         WRITE( NOUT, FMT = 9998 )TRANS, N, NRHS, IMAT,
440      $                     K, RESULT( K )
441                         NFAIL = NFAIL + 1
442                      END IF
443    70             CONTINUE
444                   NRUN = NRUN + 5
445    80          CONTINUE
446    90       CONTINUE
447 *
448   100    CONTINUE
449   110 CONTINUE
450 *
451 *     Print a summary of the results.
452 *
453       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
454 *
455  9999 FORMAT12X'N =', I5, ','10X' type ', I2, ', test(', I2,
456      $      ') = 'G12.5 )
457  9998 FORMAT' TRANS=''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
458      $      I2, ', test(', I2, ') = 'G12.5 )
459  9997 FORMAT' NORM =''', A1, ''', N =', I5, ','10X' type ', I2,
460      $      ', test(', I2, ') = 'G12.5 )
461       RETURN
462 *
463 *     End of DCHKGT
464 *
465       END