1*> \brief \b DCHKGE
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       SUBROUTINE DCHKGE( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS,
12*                          NSVAL, THRESH, TSTERR, NMAX, A, AFAC, AINV, B,
13*                          X, XACT, WORK, RWORK, IWORK, NOUT )
14*
15*       .. Scalar Arguments ..
16*       LOGICAL            TSTERR
17*       INTEGER            NM, NMAX, NN, NNB, NNS, NOUT
18*       DOUBLE PRECISION   THRESH
19*       ..
20*       .. Array Arguments ..
21*       LOGICAL            DOTYPE( * )
22*       INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
23*      $                   NVAL( * )
24*       DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
25*      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
26*       ..
27*
28*
29*> \par Purpose:
30*  =============
31*>
32*> \verbatim
33*>
34*> DCHKGE tests DGETRF, -TRI, -TRS, -RFS, and -CON.
35*> \endverbatim
36*
37*  Arguments:
38*  ==========
39*
40*> \param[in] DOTYPE
41*> \verbatim
42*>          DOTYPE is LOGICAL array, dimension (NTYPES)
43*>          The matrix types to be used for testing.  Matrices of type j
44*>          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
45*>          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
46*> \endverbatim
47*>
48*> \param[in] NM
49*> \verbatim
50*>          NM is INTEGER
51*>          The number of values of M contained in the vector MVAL.
52*> \endverbatim
53*>
54*> \param[in] MVAL
55*> \verbatim
56*>          MVAL is INTEGER array, dimension (NM)
57*>          The values of the matrix row dimension M.
58*> \endverbatim
59*>
60*> \param[in] NN
61*> \verbatim
62*>          NN is INTEGER
63*>          The number of values of N contained in the vector NVAL.
64*> \endverbatim
65*>
66*> \param[in] NVAL
67*> \verbatim
68*>          NVAL is INTEGER array, dimension (NN)
69*>          The values of the matrix column dimension N.
70*> \endverbatim
71*>
72*> \param[in] NNB
73*> \verbatim
74*>          NNB is INTEGER
75*>          The number of values of NB contained in the vector NBVAL.
76*> \endverbatim
77*>
78*> \param[in] NBVAL
79*> \verbatim
80*>          NBVAL is INTEGER array, dimension (NNB)
81*>          The values of the blocksize NB.
82*> \endverbatim
83*>
84*> \param[in] NNS
85*> \verbatim
86*>          NNS is INTEGER
87*>          The number of values of NRHS contained in the vector NSVAL.
88*> \endverbatim
89*>
90*> \param[in] NSVAL
91*> \verbatim
92*>          NSVAL is INTEGER array, dimension (NNS)
93*>          The values of the number of right hand sides NRHS.
94*> \endverbatim
95*>
96*> \param[in] THRESH
97*> \verbatim
98*>          THRESH is DOUBLE PRECISION
99*>          The threshold value for the test ratios.  A result is
100*>          included in the output file if RESULT >= THRESH.  To have
101*>          every test ratio printed, use THRESH = 0.
102*> \endverbatim
103*>
104*> \param[in] TSTERR
105*> \verbatim
106*>          TSTERR is LOGICAL
107*>          Flag that indicates whether error exits are to be tested.
108*> \endverbatim
109*>
110*> \param[in] NMAX
111*> \verbatim
112*>          NMAX is INTEGER
113*>          The maximum value permitted for M or N, used in dimensioning
114*>          the work arrays.
115*> \endverbatim
116*>
117*> \param[out] A
118*> \verbatim
119*>          A is DOUBLE PRECISION array, dimension (NMAX*NMAX)
120*> \endverbatim
121*>
122*> \param[out] AFAC
123*> \verbatim
124*>          AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX)
125*> \endverbatim
126*>
127*> \param[out] AINV
128*> \verbatim
129*>          AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
130*> \endverbatim
131*>
132*> \param[out] B
133*> \verbatim
134*>          B is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
135*>          where NSMAX is the largest entry in NSVAL.
136*> \endverbatim
137*>
138*> \param[out] X
139*> \verbatim
140*>          X is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
141*> \endverbatim
142*>
143*> \param[out] XACT
144*> \verbatim
145*>          XACT is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
146*> \endverbatim
147*>
148*> \param[out] WORK
149*> \verbatim
150*>          WORK is DOUBLE PRECISION array, dimension
151*>                      (NMAX*max(3,NSMAX))
152*> \endverbatim
153*>
154*> \param[out] RWORK
155*> \verbatim
156*>          RWORK is DOUBLE PRECISION array, dimension
157*>                      (max(2*NMAX,2*NSMAX+NWORK))
158*> \endverbatim
159*>
160*> \param[out] IWORK
161*> \verbatim
162*>          IWORK is INTEGER array, dimension (2*NMAX)
163*> \endverbatim
164*>
165*> \param[in] NOUT
166*> \verbatim
167*>          NOUT is INTEGER
168*>          The unit number for output.
169*> \endverbatim
170*
171*  Authors:
172*  ========
173*
174*> \author Univ. of Tennessee
175*> \author Univ. of California Berkeley
176*> \author Univ. of Colorado Denver
177*> \author NAG Ltd.
178*
179*> \ingroup double_lin
180*
181*  =====================================================================
182      SUBROUTINE DCHKGE( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS,
183     $                   NSVAL, THRESH, TSTERR, NMAX, A, AFAC, AINV, B,
184     $                   X, XACT, WORK, RWORK, IWORK, NOUT )
185*
186*  -- LAPACK test routine --
187*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
188*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190*     .. Scalar Arguments ..
191      LOGICAL            TSTERR
192      INTEGER            NM, NMAX, NN, NNB, NNS, NOUT
193      DOUBLE PRECISION   THRESH
194*     ..
195*     .. Array Arguments ..
196      LOGICAL            DOTYPE( * )
197      INTEGER            IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
198     $                   NVAL( * )
199      DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
200     $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
201*     ..
202*
203*  =====================================================================
204*
205*     .. Parameters ..
206      DOUBLE PRECISION   ONE, ZERO
207      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
208      INTEGER            NTYPES
209      PARAMETER          ( NTYPES = 11 )
210      INTEGER            NTESTS
211      PARAMETER          ( NTESTS = 8 )
212      INTEGER            NTRAN
213      PARAMETER          ( NTRAN = 3 )
214*     ..
215*     .. Local Scalars ..
216      LOGICAL            TRFCON, ZEROT
217      CHARACTER          DIST, NORM, TRANS, TYPE, XTYPE
218      CHARACTER*3        PATH
219      INTEGER            I, IM, IMAT, IN, INB, INFO, IOFF, IRHS, ITRAN,
220     $                   IZERO, K, KL, KU, LDA, LWORK, M, MODE, N, NB,
221     $                   NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
222      DOUBLE PRECISION   AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, DUMMY,
223     $                   RCOND, RCONDC, RCONDI, RCONDO
224*     ..
225*     .. Local Arrays ..
226      CHARACTER          TRANSS( NTRAN )
227      INTEGER            ISEED( 4 ), ISEEDY( 4 )
228      DOUBLE PRECISION   RESULT( NTESTS )
229*     ..
230*     .. External Functions ..
231      DOUBLE PRECISION   DGET06, DLANGE
232      EXTERNAL           DGET06, DLANGE
233*     ..
234*     .. External Subroutines ..
235      EXTERNAL           ALAERH, ALAHD, ALASUM, DERRGE, DGECON, DGERFS,
236     $                   DGET01, DGET02, DGET03, DGET04, DGET07, DGETRF,
237     $                   DGETRI, DGETRS, DLACPY, DLARHS, DLASET, DLATB4,
238     $                   DLATMS, XLAENV
239*     ..
240*     .. Intrinsic Functions ..
241      INTRINSIC          MAX, MIN
242*     ..
243*     .. Scalars in Common ..
244      LOGICAL            LERR, OK
245      CHARACTER*32       SRNAMT
246      INTEGER            INFOT, NUNIT
247*     ..
248*     .. Common blocks ..
249      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
250      COMMON             / SRNAMC / SRNAMT
251*     ..
252*     .. Data statements ..
253      DATA               ISEEDY / 1988, 1989, 1990, 1991 / ,
254     $                   TRANSS / 'N', 'T', 'C' /
255*     ..
256*     .. Executable Statements ..
257*
258*     Initialize constants and the random number seed.
259*
260      PATH( 1: 1 ) = 'Double precision'
261      PATH( 2: 3 ) = 'GE'
262      NRUN = 0
263      NFAIL = 0
264      NERRS = 0
265      DO 10 I = 1, 4
266         ISEED( I ) = ISEEDY( I )
267   10 CONTINUE
268*
269*     Test the error exits
270*
271      CALL XLAENV( 1, 1 )
272      IF( TSTERR )
273     $   CALL DERRGE( PATH, NOUT )
274      INFOT = 0
275      CALL XLAENV( 2, 2 )
276*
277*     Do for each value of M in MVAL
278*
279      DO 120 IM = 1, NM
280         M = MVAL( IM )
281         LDA = MAX( 1, M )
282*
283*        Do for each value of N in NVAL
284*
285         DO 110 IN = 1, NN
286            N = NVAL( IN )
287            XTYPE = 'N'
288            NIMAT = NTYPES
289            IF( M.LE.0 .OR. N.LE.0 )
290     $         NIMAT = 1
291*
292            DO 100 IMAT = 1, NIMAT
293*
294*              Do the tests only if DOTYPE( IMAT ) is true.
295*
296               IF( .NOT.DOTYPE( IMAT ) )
297     $            GO TO 100
298*
299*              Skip types 5, 6, or 7 if the matrix size is too small.
300*
301               ZEROT = IMAT.GE.5 .AND. IMAT.LE.7
302               IF( ZEROT .AND. N.LT.IMAT-4 )
303     $            GO TO 100
304*
305*              Set up parameters with DLATB4 and generate a test matrix
306*              with DLATMS.
307*
308               CALL DLATB4( PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE,
309     $                      CNDNUM, DIST )
310*
311               SRNAMT = 'DLATMS'
312               CALL DLATMS( M, N, DIST, ISEED, TYPE, RWORK, MODE,
313     $                      CNDNUM, ANORM, KL, KU, 'No packing', A, LDA,
314     $                      WORK, INFO )
315*
316*              Check error code from DLATMS.
317*
318               IF( INFO.NE.0 ) THEN
319                  CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', M, N, -1,
320     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
321                  GO TO 100
322               END IF
323*
324*              For types 5-7, zero one or more columns of the matrix to
325*              test that INFO is returned correctly.
326*
327               IF( ZEROT ) THEN
328                  IF( IMAT.EQ.5 ) THEN
329                     IZERO = 1
330                  ELSE IF( IMAT.EQ.6 ) THEN
331                     IZERO = MIN( M, N )
332                  ELSE
333                     IZERO = MIN( M, N ) / 2 + 1
334                  END IF
335                  IOFF = ( IZERO-1 )*LDA
336                  IF( IMAT.LT.7 ) THEN
337                     DO 20 I = 1, M
338                        A( IOFF+I ) = ZERO
339   20                CONTINUE
340                  ELSE
341                     CALL DLASET( 'Full', M, N-IZERO+1, ZERO, ZERO,
342     $                            A( IOFF+1 ), LDA )
343                  END IF
344               ELSE
345                  IZERO = 0
346               END IF
347*
348*              These lines, if used in place of the calls in the DO 60
349*              loop, cause the code to bomb on a Sun SPARCstation.
350*
351*               ANORMO = DLANGE( 'O', M, N, A, LDA, RWORK )
352*               ANORMI = DLANGE( 'I', M, N, A, LDA, RWORK )
353*
354*              Do for each blocksize in NBVAL
355*
356               DO 90 INB = 1, NNB
357                  NB = NBVAL( INB )
358                  CALL XLAENV( 1, NB )
359*
360*                 Compute the LU factorization of the matrix.
361*
362                  CALL DLACPY( 'Full', M, N, A, LDA, AFAC, LDA )
363                  SRNAMT = 'DGETRF'
364                  CALL DGETRF( M, N, AFAC, LDA, IWORK, INFO )
365*
366*                 Check error code from DGETRF.
367*
368                  IF( INFO.NE.IZERO )
369     $               CALL ALAERH( PATH, 'DGETRF', INFO, IZERO, ' ', M,
370     $                            N, -1, -1, NB, IMAT, NFAIL, NERRS,
371     $                            NOUT )
372                  TRFCON = .FALSE.
373*
374*+    TEST 1
375*                 Reconstruct matrix from factors and compute residual.
376*
377                  CALL DLACPY( 'Full', M, N, AFAC, LDA, AINV, LDA )
378                  CALL DGET01( M, N, A, LDA, AINV, LDA, IWORK, RWORK,
379     $                         RESULT( 1 ) )
380                  NT = 1
381*
382*+    TEST 2
383*                 Form the inverse if the factorization was successful
384*                 and compute the residual.
385*
386                  IF( M.EQ.N .AND. INFO.EQ.0 ) THEN
387                     CALL DLACPY( 'Full', N, N, AFAC, LDA, AINV, LDA )
388                     SRNAMT = 'DGETRI'
389                     NRHS = NSVAL( 1 )
390                     LWORK = NMAX*MAX( 3, NRHS )
391                     CALL DGETRI( N, AINV, LDA, IWORK, WORK, LWORK,
392     $                            INFO )
393*
394*                    Check error code from DGETRI.
395*
396                     IF( INFO.NE.0 )
397     $                  CALL ALAERH( PATH, 'DGETRI', INFO, 0, ' ', N, N,
398     $                               -1, -1, NB, IMAT, NFAIL, NERRS,
399     $                               NOUT )
400*
401*                    Compute the residual for the matrix times its
402*                    inverse.  Also compute the 1-norm condition number
403*                    of A.
404*
405                     CALL DGET03( N, A, LDA, AINV, LDA, WORK, LDA,
406     $                            RWORK, RCONDO, RESULT( 2 ) )
407                     ANORMO = DLANGE( 'O', M, N, A, LDA, RWORK )
408*
409*                    Compute the infinity-norm condition number of A.
410*
411                     ANORMI = DLANGE( 'I', M, N, A, LDA, RWORK )
412                     AINVNM = DLANGE( 'I', N, N, AINV, LDA, RWORK )
413                     IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
414                        RCONDI = ONE
415                     ELSE
416                        RCONDI = ( ONE / ANORMI ) / AINVNM
417                     END IF
418                     NT = 2
419                  ELSE
420*
421*                    Do only the condition estimate if INFO > 0.
422*
423                     TRFCON = .TRUE.
424                     ANORMO = DLANGE( 'O', M, N, A, LDA, RWORK )
425                     ANORMI = DLANGE( 'I', M, N, A, LDA, RWORK )
426                     RCONDO = ZERO
427                     RCONDI = ZERO
428                  END IF
429*
430*                 Print information about the tests so far that did not
431*                 pass the threshold.
432*
433                  DO 30 K = 1, NT
434                     IF( RESULT( K ).GE.THRESH ) THEN
435                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
436     $                     CALL ALAHD( NOUT, PATH )
437                        WRITE( NOUT, FMT = 9999 )M, N, NB, IMAT, K,
438     $                     RESULT( K )
439                        NFAIL = NFAIL + 1
440                     END IF
441   30             CONTINUE
442                  NRUN = NRUN + NT
443*
444*                 Skip the remaining tests if this is not the first
445*                 block size or if M .ne. N.  Skip the solve tests if
446*                 the matrix is singular.
447*
448                  IF( INB.GT.1 .OR. M.NE.N )
449     $               GO TO 90
450                  IF( TRFCON )
451     $               GO TO 70
452*
453                  DO 60 IRHS = 1, NNS
454                     NRHS = NSVAL( IRHS )
455                     XTYPE = 'N'
456*
457                     DO 50 ITRAN = 1, NTRAN
458                        TRANS = TRANSS( ITRAN )
459                        IF( ITRAN.EQ.1 ) THEN
460                           RCONDC = RCONDO
461                        ELSE
462                           RCONDC = RCONDI
463                        END IF
464*
465*+    TEST 3
466*                       Solve and compute residual for A * X = B.
467*
468                        SRNAMT = 'DLARHS'
469                        CALL DLARHS( PATH, XTYPE, ' ', TRANS, N, N, KL,
470     $                               KU, NRHS, A, LDA, XACT, LDA, B,
471     $                               LDA, ISEED, INFO )
472                        XTYPE = 'C'
473*
474                        CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
475                        SRNAMT = 'DGETRS'
476                        CALL DGETRS( TRANS, N, NRHS, AFAC, LDA, IWORK,
477     $                               X, LDA, INFO )
478*
479*                       Check error code from DGETRS.
480*
481                        IF( INFO.NE.0 )
482     $                     CALL ALAERH( PATH, 'DGETRS', INFO, 0, TRANS,
483     $                                  N, N, -1, -1, NRHS, IMAT, NFAIL,
484     $                                  NERRS, NOUT )
485*
486                        CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK,
487     $                               LDA )
488                        CALL DGET02( TRANS, N, N, NRHS, A, LDA, X, LDA,
489     $                               WORK, LDA, RWORK, RESULT( 3 ) )
490*
491*+    TEST 4
492*                       Check solution from generated exact solution.
493*
494                        CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
495     $                               RESULT( 4 ) )
496*
497*+    TESTS 5, 6, and 7
498*                       Use iterative refinement to improve the
499*                       solution.
500*
501                        SRNAMT = 'DGERFS'
502                        CALL DGERFS( TRANS, N, NRHS, A, LDA, AFAC, LDA,
503     $                               IWORK, B, LDA, X, LDA, RWORK,
504     $                               RWORK( NRHS+1 ), WORK,
505     $                               IWORK( N+1 ), INFO )
506*
507*                       Check error code from DGERFS.
508*
509                        IF( INFO.NE.0 )
510     $                     CALL ALAERH( PATH, 'DGERFS', INFO, 0, TRANS,
511     $                                  N, N, -1, -1, NRHS, IMAT, NFAIL,
512     $                                  NERRS, NOUT )
513*
514                        CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
515     $                               RESULT( 5 ) )
516                        CALL DGET07( TRANS, N, NRHS, A, LDA, B, LDA, X,
517     $                               LDA, XACT, LDA, RWORK, .TRUE.,
518     $                               RWORK( NRHS+1 ), RESULT( 6 ) )
519*
520*                       Print information about the tests that did not
521*                       pass the threshold.
522*
523                        DO 40 K = 3, 7
524                           IF( RESULT( K ).GE.THRESH ) THEN
525                              IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
526     $                           CALL ALAHD( NOUT, PATH )
527                              WRITE( NOUT, FMT = 9998 )TRANS, N, NRHS,
528     $                           IMAT, K, RESULT( K )
529                              NFAIL = NFAIL + 1
530                           END IF
531   40                   CONTINUE
532                        NRUN = NRUN + 5
533   50                CONTINUE
534   60             CONTINUE
535*
536*+    TEST 8
537*                    Get an estimate of RCOND = 1/CNDNUM.
538*
539   70             CONTINUE
540                  DO 80 ITRAN = 1, 2
541                     IF( ITRAN.EQ.1 ) THEN
542                        ANORM = ANORMO
543                        RCONDC = RCONDO
544                        NORM = 'O'
545                     ELSE
546                        ANORM = ANORMI
547                        RCONDC = RCONDI
548                        NORM = 'I'
549                     END IF
550                     SRNAMT = 'DGECON'
551                     CALL DGECON( NORM, N, AFAC, LDA, ANORM, RCOND,
552     $                            WORK, IWORK( N+1 ), INFO )
553*
554*                       Check error code from DGECON.
555*
556                     IF( INFO.NE.0 )
557     $                  CALL ALAERH( PATH, 'DGECON', INFO, 0, NORM, N,
558     $                               N, -1, -1, -1, IMAT, NFAIL, NERRS,
559     $                               NOUT )
560*
561*                       This line is needed on a Sun SPARCstation.
562*
563                     DUMMY = RCOND
564*
565                     RESULT( 8 ) = DGET06( RCOND, RCONDC )
566*
567*                    Print information about the tests that did not pass
568*                    the threshold.
569*
570                     IF( RESULT( 8 ).GE.THRESH ) THEN
571                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
572     $                     CALL ALAHD( NOUT, PATH )
573                        WRITE( NOUT, FMT = 9997 )NORM, N, IMAT, 8,
574     $                     RESULT( 8 )
575                        NFAIL = NFAIL + 1
576                     END IF
577                     NRUN = NRUN + 1
578   80             CONTINUE
579   90          CONTINUE
580  100       CONTINUE
581  110    CONTINUE
582  120 CONTINUE
583*
584*     Print a summary of the results.
585*
586      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
587*
588 9999 FORMAT( ' M = ', I5, ', N =', I5, ', NB =', I4, ', type ', I2,
589     $      ', test(', I2, ') =', G12.5 )
590 9998 FORMAT( ' TRANS=''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
591     $      I2, ', test(', I2, ') =', G12.5 )
592 9997 FORMAT( ' NORM =''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
593     $      ', test(', I2, ') =', G12.5 )
594      RETURN
595*
596*     End of DCHKGE
597*
598      END
599