LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dchkge.f
Go to the documentation of this file.
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 (NBVAL)
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 *> \date November 2011
180 *
181 *> \ingroup double_lin
182 *
183 * =====================================================================
184  SUBROUTINE dchkge( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS,
185  $ nsval, thresh, tsterr, nmax, a, afac, ainv, b,
186  $ x, xact, work, rwork, iwork, nout )
187 *
188 * -- LAPACK test routine (version 3.4.0) --
189 * -- LAPACK is a software package provided by Univ. of Tennessee, --
190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191 * November 2011
192 *
193 * .. Scalar Arguments ..
194  LOGICAL tsterr
195  INTEGER nm, nmax, nn, nnb, nns, nout
196  DOUBLE PRECISION thresh
197 * ..
198 * .. Array Arguments ..
199  LOGICAL dotype( * )
200  INTEGER iwork( * ), mval( * ), nbval( * ), nsval( * ),
201  $ nval( * )
202  DOUBLE PRECISION a( * ), afac( * ), ainv( * ), b( * ),
203  $ rwork( * ), work( * ), x( * ), xact( * )
204 * ..
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209  DOUBLE PRECISION one, zero
210  parameter( one = 1.0d+0, zero = 0.0d+0 )
211  INTEGER ntypes
212  parameter( ntypes = 11 )
213  INTEGER ntests
214  parameter( ntests = 8 )
215  INTEGER ntran
216  parameter( ntran = 3 )
217 * ..
218 * .. Local Scalars ..
219  LOGICAL trfcon, zerot
220  CHARACTER dist, norm, trans, type, xtype
221  CHARACTER*3 path
222  INTEGER i, im, imat, in, inb, info, ioff, irhs, itran,
223  $ izero, k, kl, ku, lda, lwork, m, mode, n, nb,
224  $ nerrs, nfail, nimat, nrhs, nrun, nt
225  DOUBLE PRECISION ainvnm, anorm, anormi, anormo, cndnum, dummy,
226  $ rcond, rcondc, rcondi, rcondo
227 * ..
228 * .. Local Arrays ..
229  CHARACTER transs( ntran )
230  INTEGER iseed( 4 ), iseedy( 4 )
231  DOUBLE PRECISION result( ntests )
232 * ..
233 * .. External Functions ..
234  DOUBLE PRECISION dget06, dlange
235  EXTERNAL dget06, dlange
236 * ..
237 * .. External Subroutines ..
238  EXTERNAL alaerh, alahd, alasum, derrge, dgecon, dgerfs,
239  $ dget01, dget02, dget03, dget04, dget07, dgetrf,
240  $ dgetri, dgetrs, dlacpy, dlarhs, dlaset, dlatb4,
241  $ dlatms, xlaenv
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC max, min
245 * ..
246 * .. Scalars in Common ..
247  LOGICAL lerr, ok
248  CHARACTER*32 srnamt
249  INTEGER infot, nunit
250 * ..
251 * .. Common blocks ..
252  common / infoc / infot, nunit, ok, lerr
253  common / srnamc / srnamt
254 * ..
255 * .. Data statements ..
256  DATA iseedy / 1988, 1989, 1990, 1991 / ,
257  $ transs / 'N', 'T', 'C' /
258 * ..
259 * .. Executable Statements ..
260 *
261 * Initialize constants and the random number seed.
262 *
263  path( 1: 1 ) = 'Double precision'
264  path( 2: 3 ) = 'GE'
265  nrun = 0
266  nfail = 0
267  nerrs = 0
268  DO 10 i = 1, 4
269  iseed( i ) = iseedy( i )
270  10 continue
271 *
272 * Test the error exits
273 *
274  CALL xlaenv( 1, 1 )
275  IF( tsterr )
276  $ CALL derrge( path, nout )
277  infot = 0
278  CALL xlaenv( 2, 2 )
279 *
280 * Do for each value of M in MVAL
281 *
282  DO 120 im = 1, nm
283  m = mval( im )
284  lda = max( 1, m )
285 *
286 * Do for each value of N in NVAL
287 *
288  DO 110 in = 1, nn
289  n = nval( in )
290  xtype = 'N'
291  nimat = ntypes
292  IF( m.LE.0 .OR. n.LE.0 )
293  $ nimat = 1
294 *
295  DO 100 imat = 1, nimat
296 *
297 * Do the tests only if DOTYPE( IMAT ) is true.
298 *
299  IF( .NOT.dotype( imat ) )
300  $ go to 100
301 *
302 * Skip types 5, 6, or 7 if the matrix size is too small.
303 *
304  zerot = imat.GE.5 .AND. imat.LE.7
305  IF( zerot .AND. n.LT.imat-4 )
306  $ go to 100
307 *
308 * Set up parameters with DLATB4 and generate a test matrix
309 * with DLATMS.
310 *
311  CALL dlatb4( path, imat, m, n, type, kl, ku, anorm, mode,
312  $ cndnum, dist )
313 *
314  srnamt = 'DLATMS'
315  CALL dlatms( m, n, dist, iseed, type, rwork, mode,
316  $ cndnum, anorm, kl, ku, 'No packing', a, lda,
317  $ work, info )
318 *
319 * Check error code from DLATMS.
320 *
321  IF( info.NE.0 ) THEN
322  CALL alaerh( path, 'DLATMS', info, 0, ' ', m, n, -1,
323  $ -1, -1, imat, nfail, nerrs, nout )
324  go to 100
325  END IF
326 *
327 * For types 5-7, zero one or more columns of the matrix to
328 * test that INFO is returned correctly.
329 *
330  IF( zerot ) THEN
331  IF( imat.EQ.5 ) THEN
332  izero = 1
333  ELSE IF( imat.EQ.6 ) THEN
334  izero = min( m, n )
335  ELSE
336  izero = min( m, n ) / 2 + 1
337  END IF
338  ioff = ( izero-1 )*lda
339  IF( imat.LT.7 ) THEN
340  DO 20 i = 1, m
341  a( ioff+i ) = zero
342  20 continue
343  ELSE
344  CALL dlaset( 'Full', m, n-izero+1, zero, zero,
345  $ a( ioff+1 ), lda )
346  END IF
347  ELSE
348  izero = 0
349  END IF
350 *
351 * These lines, if used in place of the calls in the DO 60
352 * loop, cause the code to bomb on a Sun SPARCstation.
353 *
354 * ANORMO = DLANGE( 'O', M, N, A, LDA, RWORK )
355 * ANORMI = DLANGE( 'I', M, N, A, LDA, RWORK )
356 *
357 * Do for each blocksize in NBVAL
358 *
359  DO 90 inb = 1, nnb
360  nb = nbval( inb )
361  CALL xlaenv( 1, nb )
362 *
363 * Compute the LU factorization of the matrix.
364 *
365  CALL dlacpy( 'Full', m, n, a, lda, afac, lda )
366  srnamt = 'DGETRF'
367  CALL dgetrf( m, n, afac, lda, iwork, info )
368 *
369 * Check error code from DGETRF.
370 *
371  IF( info.NE.izero )
372  $ CALL alaerh( path, 'DGETRF', info, izero, ' ', m,
373  $ n, -1, -1, nb, imat, nfail, nerrs,
374  $ nout )
375  trfcon = .false.
376 *
377 *+ TEST 1
378 * Reconstruct matrix from factors and compute residual.
379 *
380  CALL dlacpy( 'Full', m, n, afac, lda, ainv, lda )
381  CALL dget01( m, n, a, lda, ainv, lda, iwork, rwork,
382  $ result( 1 ) )
383  nt = 1
384 *
385 *+ TEST 2
386 * Form the inverse if the factorization was successful
387 * and compute the residual.
388 *
389  IF( m.EQ.n .AND. info.EQ.0 ) THEN
390  CALL dlacpy( 'Full', n, n, afac, lda, ainv, lda )
391  srnamt = 'DGETRI'
392  nrhs = nsval( 1 )
393  lwork = nmax*max( 3, nrhs )
394  CALL dgetri( n, ainv, lda, iwork, work, lwork,
395  $ info )
396 *
397 * Check error code from DGETRI.
398 *
399  IF( info.NE.0 )
400  $ CALL alaerh( path, 'DGETRI', info, 0, ' ', n, n,
401  $ -1, -1, nb, imat, nfail, nerrs,
402  $ nout )
403 *
404 * Compute the residual for the matrix times its
405 * inverse. Also compute the 1-norm condition number
406 * of A.
407 *
408  CALL dget03( n, a, lda, ainv, lda, work, lda,
409  $ rwork, rcondo, result( 2 ) )
410  anormo = dlange( 'O', m, n, a, lda, rwork )
411 *
412 * Compute the infinity-norm condition number of A.
413 *
414  anormi = dlange( 'I', m, n, a, lda, rwork )
415  ainvnm = dlange( 'I', n, n, ainv, lda, rwork )
416  IF( anormi.LE.zero .OR. ainvnm.LE.zero ) THEN
417  rcondi = one
418  ELSE
419  rcondi = ( one / anormi ) / ainvnm
420  END IF
421  nt = 2
422  ELSE
423 *
424 * Do only the condition estimate if INFO > 0.
425 *
426  trfcon = .true.
427  anormo = dlange( 'O', m, n, a, lda, rwork )
428  anormi = dlange( 'I', m, n, a, lda, rwork )
429  rcondo = zero
430  rcondi = zero
431  END IF
432 *
433 * Print information about the tests so far that did not
434 * pass the threshold.
435 *
436  DO 30 k = 1, nt
437  IF( result( k ).GE.thresh ) THEN
438  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
439  $ CALL alahd( nout, path )
440  WRITE( nout, fmt = 9999 )m, n, nb, imat, k,
441  $ result( k )
442  nfail = nfail + 1
443  END IF
444  30 continue
445  nrun = nrun + nt
446 *
447 * Skip the remaining tests if this is not the first
448 * block size or if M .ne. N. Skip the solve tests if
449 * the matrix is singular.
450 *
451  IF( inb.GT.1 .OR. m.NE.n )
452  $ go to 90
453  IF( trfcon )
454  $ go to 70
455 *
456  DO 60 irhs = 1, nns
457  nrhs = nsval( irhs )
458  xtype = 'N'
459 *
460  DO 50 itran = 1, ntran
461  trans = transs( itran )
462  IF( itran.EQ.1 ) THEN
463  rcondc = rcondo
464  ELSE
465  rcondc = rcondi
466  END IF
467 *
468 *+ TEST 3
469 * Solve and compute residual for A * X = B.
470 *
471  srnamt = 'DLARHS'
472  CALL dlarhs( path, xtype, ' ', trans, n, n, kl,
473  $ ku, nrhs, a, lda, xact, lda, b,
474  $ lda, iseed, info )
475  xtype = 'C'
476 *
477  CALL dlacpy( 'Full', n, nrhs, b, lda, x, lda )
478  srnamt = 'DGETRS'
479  CALL dgetrs( trans, n, nrhs, afac, lda, iwork,
480  $ x, lda, info )
481 *
482 * Check error code from DGETRS.
483 *
484  IF( info.NE.0 )
485  $ CALL alaerh( path, 'DGETRS', info, 0, trans,
486  $ n, n, -1, -1, nrhs, imat, nfail,
487  $ nerrs, nout )
488 *
489  CALL dlacpy( 'Full', n, nrhs, b, lda, work,
490  $ lda )
491  CALL dget02( trans, n, n, nrhs, a, lda, x, lda,
492  $ work, lda, rwork, result( 3 ) )
493 *
494 *+ TEST 4
495 * Check solution from generated exact solution.
496 *
497  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
498  $ result( 4 ) )
499 *
500 *+ TESTS 5, 6, and 7
501 * Use iterative refinement to improve the
502 * solution.
503 *
504  srnamt = 'DGERFS'
505  CALL dgerfs( trans, n, nrhs, a, lda, afac, lda,
506  $ iwork, b, lda, x, lda, rwork,
507  $ rwork( nrhs+1 ), work,
508  $ iwork( n+1 ), info )
509 *
510 * Check error code from DGERFS.
511 *
512  IF( info.NE.0 )
513  $ CALL alaerh( path, 'DGERFS', info, 0, trans,
514  $ n, n, -1, -1, nrhs, imat, nfail,
515  $ nerrs, nout )
516 *
517  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
518  $ result( 5 ) )
519  CALL dget07( trans, n, nrhs, a, lda, b, lda, x,
520  $ lda, xact, lda, rwork, .true.,
521  $ rwork( nrhs+1 ), result( 6 ) )
522 *
523 * Print information about the tests that did not
524 * pass the threshold.
525 *
526  DO 40 k = 3, 7
527  IF( result( k ).GE.thresh ) THEN
528  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
529  $ CALL alahd( nout, path )
530  WRITE( nout, fmt = 9998 )trans, n, nrhs,
531  $ imat, k, result( k )
532  nfail = nfail + 1
533  END IF
534  40 continue
535  nrun = nrun + 5
536  50 continue
537  60 continue
538 *
539 *+ TEST 8
540 * Get an estimate of RCOND = 1/CNDNUM.
541 *
542  70 continue
543  DO 80 itran = 1, 2
544  IF( itran.EQ.1 ) THEN
545  anorm = anormo
546  rcondc = rcondo
547  norm = 'O'
548  ELSE
549  anorm = anormi
550  rcondc = rcondi
551  norm = 'I'
552  END IF
553  srnamt = 'DGECON'
554  CALL dgecon( norm, n, afac, lda, anorm, rcond,
555  $ work, iwork( n+1 ), info )
556 *
557 * Check error code from DGECON.
558 *
559  IF( info.NE.0 )
560  $ CALL alaerh( path, 'DGECON', info, 0, norm, n,
561  $ n, -1, -1, -1, imat, nfail, nerrs,
562  $ nout )
563 *
564 * This line is needed on a Sun SPARCstation.
565 *
566  dummy = rcond
567 *
568  result( 8 ) = dget06( rcond, rcondc )
569 *
570 * Print information about the tests that did not pass
571 * the threshold.
572 *
573  IF( result( 8 ).GE.thresh ) THEN
574  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
575  $ CALL alahd( nout, path )
576  WRITE( nout, fmt = 9997 )norm, n, imat, 8,
577  $ result( 8 )
578  nfail = nfail + 1
579  END IF
580  nrun = nrun + 1
581  80 continue
582  90 continue
583  100 continue
584  110 continue
585  120 continue
586 *
587 * Print a summary of the results.
588 *
589  CALL alasum( path, nout, nfail, nrun, nerrs )
590 *
591  9999 format( ' M = ', i5, ', N =', i5, ', NB =', i4, ', type ', i2,
592  $ ', test(', i2, ') =', g12.5 )
593  9998 format( ' TRANS=''', a1, ''', N =', i5, ', NRHS=', i3, ', type ',
594  $ i2, ', test(', i2, ') =', g12.5 )
595  9997 format( ' NORM =''', a1, ''', N =', i5, ',', 10x, ' type ', i2,
596  $ ', test(', i2, ') =', g12.5 )
597  return
598 *
599 * End of DCHKGE
600 *
601  END