MAGMA  1.2.0
MatrixAlgebraonGPUandMulticoreArchitectures
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
dchkgb.f
Go to the documentation of this file.
1  SUBROUTINE dchkgb( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS,
2  $ nsval, thresh, tsterr, a, la, afac, lafac, b,
3  $ x, xact, work, rwork, iwork, 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 la, lafac, nm, nn, nnb, nns, nout
12  DOUBLE PRECISION thresh
13 * ..
14 * .. Array Arguments ..
15  LOGICAL dotype( * )
16  INTEGER iwork( * ), mval( * ), nbval( * ), nsval( * ),
17  $ nval( * )
18  DOUBLE PRECISION a( * ), afac( * ), b( * ), rwork( * ),
19  $ work( * ), x( * ), xact( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DCHKGB tests DGBTRF, -TRS, -RFS, and -CON
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 * NM (input) INTEGER
36 * The number of values of M contained in the vector MVAL.
37 *
38 * MVAL (input) INTEGER array, dimension (NM)
39 * The values of the matrix row dimension M.
40 *
41 * NN (input) INTEGER
42 * The number of values of N contained in the vector NVAL.
43 *
44 * NVAL (input) INTEGER array, dimension (NN)
45 * The values of the matrix column dimension N.
46 *
47 * NNB (input) INTEGER
48 * The number of values of NB contained in the vector NBVAL.
49 *
50 * NBVAL (input) INTEGER array, dimension (NNB)
51 * The values of the blocksize NB.
52 *
53 * NNS (input) INTEGER
54 * The number of values of NRHS contained in the vector NSVAL.
55 *
56 * NSVAL (input) INTEGER array, dimension (NNS)
57 * The values of the number of right hand sides NRHS.
58 *
59 * THRESH (input) DOUBLE PRECISION
60 * The threshold value for the test ratios. A result is
61 * included in the output file if RESULT >= THRESH. To have
62 * every test ratio printed, use THRESH = 0.
63 *
64 * TSTERR (input) LOGICAL
65 * Flag that indicates whether error exits are to be tested.
66 *
67 * A (workspace) DOUBLE PRECISION array, dimension (LA)
68 *
69 * LA (input) INTEGER
70 * The length of the array A. LA >= (KLMAX+KUMAX+1)*NMAX
71 * where KLMAX is the largest entry in the local array KLVAL,
72 * KUMAX is the largest entry in the local array KUVAL and
73 * NMAX is the largest entry in the input array NVAL.
74 *
75 * AFAC (workspace) DOUBLE PRECISION array, dimension (LAFAC)
76 *
77 * LAFAC (input) INTEGER
78 * The length of the array AFAC. LAFAC >= (2*KLMAX+KUMAX+1)*NMAX
79 * where KLMAX is the largest entry in the local array KLVAL,
80 * KUMAX is the largest entry in the local array KUVAL and
81 * NMAX is the largest entry in the input array NVAL.
82 *
83 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
84 * where NSMAX is the largest entry in NSVAL.
85 *
86 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
87 *
88 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*NSMAX)
89 *
90 * WORK (workspace) DOUBLE PRECISION array, dimension
91 * (NMAX*max(3,NSMAX,NMAX))
92 *
93 * RWORK (workspace) DOUBLE PRECISION array, dimension
94 * (max(NMAX,2*NSMAX))
95 *
96 * IWORK (workspace) INTEGER array, dimension (2*NMAX)
97 *
98 * NOUT (input) INTEGER
99 * The unit number for output.
100 *
101 * =====================================================================
102 *
103 * .. Parameters ..
104  DOUBLE PRECISION one, zero
105  parameter( one = 1.0d+0, zero = 0.0d+0 )
106  INTEGER ntypes, ntests
107  parameter( ntypes = 8, ntests = 7 )
108  INTEGER nbw, ntran
109  parameter( nbw = 4, ntran = 3 )
110 * ..
111 * .. Local Scalars ..
112  LOGICAL trfcon, zerot
113  CHARACTER dist, norm, trans, type, xtype
114  CHARACTER*3 path
115  INTEGER i, i1, i2, ikl, iku, im, imat, in, inb, info,
116  $ ioff, irhs, itran, izero, j, k, kl, koff, ku,
117  $ lda, ldafac, ldb, m, mode, n, nb, nerrs, nfail,
118  $ nimat, nkl, nku, nrhs, nrun
119  DOUBLE PRECISION ainvnm, anorm, anormi, anormo, cndnum, rcond,
120  $ rcondc, rcondi, rcondo
121 * ..
122 * .. Local Arrays ..
123  CHARACTER transs( ntran )
124  INTEGER iseed( 4 ), iseedy( 4 ), klval( nbw ),
125  $ kuval( nbw )
126  DOUBLE PRECISION result( ntests )
127 * ..
128 * .. External Functions ..
129  DOUBLE PRECISION dget06, dlangb, dlange
130  EXTERNAL dget06, dlangb, dlange
131 * ..
132 * .. External Subroutines ..
133  EXTERNAL alaerh, alahd, alasum, dcopy, derrge, dgbcon,
134  $ dgbrfs, dgbt01, dgbt02, dgbt05, dgbtrf, dgbtrs,
135  $ dget04, dlacpy, dlarhs, dlaset, dlatb4, dlatms,
136  $ xlaenv
137 * ..
138 * .. Intrinsic Functions ..
139  INTRINSIC max, min
140 * ..
141 * .. Scalars in Common ..
142  LOGICAL lerr, ok
143  CHARACTER*32 srnamt
144  INTEGER infot, nunit
145 * ..
146 * .. Common blocks ..
147  common / infoc / infot, nunit, ok, lerr
148  common / srnamc / srnamt
149 * ..
150 * .. Data statements ..
151  DATA iseedy / 1988, 1989, 1990, 1991 / ,
152  $ transs / 'N', 'T', 'C' /
153 * ..
154 * .. Executable Statements ..
155 *
156 * Initialize constants and the random number seed.
157 *
158  path( 1: 1 ) = 'Double precision'
159  path( 2: 3 ) = 'GB'
160  nrun = 0
161  nfail = 0
162  nerrs = 0
163  DO 10 i = 1, 4
164  iseed( i ) = iseedy( i )
165  10 continue
166 *
167 * Test the error exits
168 *
169  IF( tsterr )
170  $ CALL derrge( path, nout )
171  infot = 0
172  CALL xlaenv( 2, 2 )
173 *
174 * Initialize the first value for the lower and upper bandwidths.
175 *
176  klval( 1 ) = 0
177  kuval( 1 ) = 0
178 *
179 * Do for each value of M in MVAL
180 *
181  DO 160 im = 1, nm
182  m = mval( im )
183 *
184 * Set values to use for the lower bandwidth.
185 *
186  klval( 2 ) = m + ( m+1 ) / 4
187 *
188 * KLVAL( 2 ) = MAX( M-1, 0 )
189 *
190  klval( 3 ) = ( 3*m-1 ) / 4
191  klval( 4 ) = ( m+1 ) / 4
192 *
193 * Do for each value of N in NVAL
194 *
195  DO 150 in = 1, nn
196  n = nval( in )
197  xtype = 'N'
198 *
199 * Set values to use for the upper bandwidth.
200 *
201  kuval( 2 ) = n + ( n+1 ) / 4
202 *
203 * KUVAL( 2 ) = MAX( N-1, 0 )
204 *
205  kuval( 3 ) = ( 3*n-1 ) / 4
206  kuval( 4 ) = ( n+1 ) / 4
207 *
208 * Set limits on the number of loop iterations.
209 *
210  nkl = min( m+1, 4 )
211  IF( n.EQ.0 )
212  $ nkl = 2
213  nku = min( n+1, 4 )
214  IF( m.EQ.0 )
215  $ nku = 2
216  nimat = ntypes
217  IF( m.LE.0 .OR. n.LE.0 )
218  $ nimat = 1
219 *
220  DO 140 ikl = 1, nkl
221 *
222 * Do for KL = 0, (5*M+1)/4, (3M-1)/4, and (M+1)/4. This
223 * order makes it easier to skip redundant values for small
224 * values of M.
225 *
226  kl = klval( ikl )
227  DO 130 iku = 1, nku
228 *
229 * Do for KU = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This
230 * order makes it easier to skip redundant values for
231 * small values of N.
232 *
233  ku = kuval( iku )
234 *
235 * Check that A and AFAC are big enough to generate this
236 * matrix.
237 *
238  lda = kl + ku + 1
239  ldafac = 2*kl + ku + 1
240  IF( ( lda*n ).GT.la .OR. ( ldafac*n ).GT.lafac ) THEN
241  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
242  $ CALL alahd( nout, path )
243  IF( n*( kl+ku+1 ).GT.la ) THEN
244  WRITE( nout, fmt = 9999 )la, m, n, kl, ku,
245  $ n*( kl+ku+1 )
246  nerrs = nerrs + 1
247  END IF
248  IF( n*( 2*kl+ku+1 ).GT.lafac ) THEN
249  WRITE( nout, fmt = 9998 )lafac, m, n, kl, ku,
250  $ n*( 2*kl+ku+1 )
251  nerrs = nerrs + 1
252  END IF
253  go to 130
254  END IF
255 *
256  DO 120 imat = 1, nimat
257 *
258 * Do the tests only if DOTYPE( IMAT ) is true.
259 *
260  IF( .NOT.dotype( imat ) )
261  $ go to 120
262 *
263 * Skip types 2, 3, or 4 if the matrix size is too
264 * small.
265 *
266  zerot = imat.GE.2 .AND. imat.LE.4
267  IF( zerot .AND. n.LT.imat-1 )
268  $ go to 120
269 *
270  IF( .NOT.zerot .OR. .NOT.dotype( 1 ) ) THEN
271 *
272 * Set up parameters with DLATB4 and generate a
273 * test matrix with DLATMS.
274 *
275  CALL dlatb4( path, imat, m, n, type, kl, ku,
276  $ anorm, mode, cndnum, dist )
277 *
278  koff = max( 1, ku+2-n )
279  DO 20 i = 1, koff - 1
280  a( i ) = zero
281  20 continue
282  srnamt = 'DLATMS'
283  CALL dlatms( m, n, dist, iseed, type, rwork,
284  $ mode, cndnum, anorm, kl, ku, 'Z',
285  $ a( koff ), lda, work, info )
286 *
287 * Check the error code from DLATMS.
288 *
289  IF( info.NE.0 ) THEN
290  CALL alaerh( path, 'DLATMS', info, 0, ' ', m,
291  $ n, kl, ku, -1, imat, nfail,
292  $ nerrs, nout )
293  go to 120
294  END IF
295  ELSE IF( izero.GT.0 ) THEN
296 *
297 * Use the same matrix for types 3 and 4 as for
298 * type 2 by copying back the zeroed out column.
299 *
300  CALL dcopy( i2-i1+1, b, 1, a( ioff+i1 ), 1 )
301  END IF
302 *
303 * For types 2, 3, and 4, zero one or more columns of
304 * the matrix to test that INFO is returned correctly.
305 *
306  izero = 0
307  IF( zerot ) THEN
308  IF( imat.EQ.2 ) THEN
309  izero = 1
310  ELSE IF( imat.EQ.3 ) THEN
311  izero = min( m, n )
312  ELSE
313  izero = min( m, n ) / 2 + 1
314  END IF
315  ioff = ( izero-1 )*lda
316  IF( imat.LT.4 ) THEN
317 *
318 * Store the column to be zeroed out in B.
319 *
320  i1 = max( 1, ku+2-izero )
321  i2 = min( kl+ku+1, ku+1+( m-izero ) )
322  CALL dcopy( i2-i1+1, a( ioff+i1 ), 1, b, 1 )
323 *
324  DO 30 i = i1, i2
325  a( ioff+i ) = zero
326  30 continue
327  ELSE
328  DO 50 j = izero, n
329  DO 40 i = max( 1, ku+2-j ),
330  $ min( kl+ku+1, ku+1+( m-j ) )
331  a( ioff+i ) = zero
332  40 continue
333  ioff = ioff + lda
334  50 continue
335  END IF
336  END IF
337 *
338 * These lines, if used in place of the calls in the
339 * loop over INB, cause the code to bomb on a Sun
340 * SPARCstation.
341 *
342 * ANORMO = DLANGB( 'O', N, KL, KU, A, LDA, RWORK )
343 * ANORMI = DLANGB( 'I', N, KL, KU, A, LDA, RWORK )
344 *
345 * Do for each blocksize in NBVAL
346 *
347  DO 110 inb = 1, nnb
348  nb = nbval( inb )
349  CALL xlaenv( 1, nb )
350 *
351 * Compute the LU factorization of the band matrix.
352 *
353  IF( m.GT.0 .AND. n.GT.0 )
354  $ CALL dlacpy( 'Full', kl+ku+1, n, a, lda,
355  $ afac( kl+1 ), ldafac )
356  srnamt = 'DGBTRF'
357  CALL dgbtrf( m, n, kl, ku, afac, ldafac, iwork,
358  $ info )
359 *
360 * Check error code from DGBTRF.
361 *
362  IF( info.NE.izero )
363  $ CALL alaerh( path, 'DGBTRF', info, izero,
364  $ ' ', m, n, kl, ku, nb, imat,
365  $ nfail, nerrs, nout )
366  trfcon = .false.
367 *
368 *+ TEST 1
369 * Reconstruct matrix from factors and compute
370 * residual.
371 *
372  CALL dgbt01( m, n, kl, ku, a, lda, afac, ldafac,
373  $ iwork, work, result( 1 ) )
374 *
375 * Print information about the tests so far that
376 * did not pass the threshold.
377 *
378  IF( result( 1 ).GE.thresh ) THEN
379  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
380  $ CALL alahd( nout, path )
381  WRITE( nout, fmt = 9997 )m, n, kl, ku, nb,
382  $ imat, 1, result( 1 )
383  nfail = nfail + 1
384  END IF
385  nrun = nrun + 1
386 *
387 * Skip the remaining tests if this is not the
388 * first block size or if M .ne. N.
389 *
390  IF( inb.GT.1 .OR. m.NE.n )
391  $ go to 110
392 *
393  anormo = dlangb( 'O', n, kl, ku, a, lda, rwork )
394  anormi = dlangb( 'I', n, kl, ku, a, lda, rwork )
395 *
396  IF( info.EQ.0 ) THEN
397 *
398 * Form the inverse of A so we can get a good
399 * estimate of CNDNUM = norm(A) * norm(inv(A)).
400 *
401  ldb = max( 1, n )
402  CALL dlaset( 'Full', n, n, zero, one, work,
403  $ ldb )
404  srnamt = 'DGBTRS'
405  CALL dgbtrs( 'No transpose', n, kl, ku, n,
406  $ afac, ldafac, iwork, work, ldb,
407  $ info )
408 *
409 * Compute the 1-norm condition number of A.
410 *
411  ainvnm = dlange( 'O', n, n, work, ldb,
412  $ rwork )
413  IF( anormo.LE.zero .OR. ainvnm.LE.zero ) THEN
414  rcondo = one
415  ELSE
416  rcondo = ( one / anormo ) / ainvnm
417  END IF
418 *
419 * Compute the infinity-norm condition number of
420 * A.
421 *
422  ainvnm = dlange( 'I', n, n, work, ldb,
423  $ rwork )
424  IF( anormi.LE.zero .OR. ainvnm.LE.zero ) THEN
425  rcondi = one
426  ELSE
427  rcondi = ( one / anormi ) / ainvnm
428  END IF
429  ELSE
430 *
431 * Do only the condition estimate if INFO.NE.0.
432 *
433  trfcon = .true.
434  rcondo = zero
435  rcondi = zero
436  END IF
437 *
438 * Skip the solve tests if the matrix is singular.
439 *
440  IF( trfcon )
441  $ go to 90
442 *
443  DO 80 irhs = 1, nns
444  nrhs = nsval( irhs )
445  xtype = 'N'
446 *
447  DO 70 itran = 1, ntran
448  trans = transs( itran )
449  IF( itran.EQ.1 ) THEN
450  rcondc = rcondo
451  norm = 'O'
452  ELSE
453  rcondc = rcondi
454  norm = 'I'
455  END IF
456 *
457 *+ TEST 2:
458 * Solve and compute residual for A * X = B.
459 *
460  srnamt = 'DLARHS'
461  CALL dlarhs( path, xtype, ' ', trans, n,
462  $ n, kl, ku, nrhs, a, lda,
463  $ xact, ldb, b, ldb, iseed,
464  $ info )
465  xtype = 'C'
466  CALL dlacpy( 'Full', n, nrhs, b, ldb, x,
467  $ ldb )
468 *
469  srnamt = 'DGBTRS'
470  CALL dgbtrs( trans, n, kl, ku, nrhs, afac,
471  $ ldafac, iwork, x, ldb, info )
472 *
473 * Check error code from DGBTRS.
474 *
475  IF( info.NE.0 )
476  $ CALL alaerh( path, 'DGBTRS', info, 0,
477  $ trans, n, n, kl, ku, -1,
478  $ imat, nfail, nerrs, nout )
479 *
480  CALL dlacpy( 'Full', n, nrhs, b, ldb,
481  $ work, ldb )
482  CALL dgbt02( trans, m, n, kl, ku, nrhs, a,
483  $ lda, x, ldb, work, ldb,
484  $ result( 2 ) )
485 *
486 *+ TEST 3:
487 * Check solution from generated exact
488 * solution.
489 *
490  CALL dget04( n, nrhs, x, ldb, xact, ldb,
491  $ rcondc, result( 3 ) )
492 *
493 *+ TESTS 4, 5, 6:
494 * Use iterative refinement to improve the
495 * solution.
496 *
497  srnamt = 'DGBRFS'
498  CALL dgbrfs( trans, n, kl, ku, nrhs, a,
499  $ lda, afac, ldafac, iwork, b,
500  $ ldb, x, ldb, rwork,
501  $ rwork( nrhs+1 ), work,
502  $ iwork( n+1 ), info )
503 *
504 * Check error code from DGBRFS.
505 *
506  IF( info.NE.0 )
507  $ CALL alaerh( path, 'DGBRFS', info, 0,
508  $ trans, n, n, kl, ku, nrhs,
509  $ imat, nfail, nerrs, nout )
510 *
511  CALL dget04( n, nrhs, x, ldb, xact, ldb,
512  $ rcondc, result( 4 ) )
513  CALL dgbt05( trans, n, kl, ku, nrhs, a,
514  $ lda, b, ldb, x, ldb, xact,
515  $ ldb, rwork, rwork( nrhs+1 ),
516  $ result( 5 ) )
517  DO 60 k = 2, 6
518  IF( result( k ).GE.thresh ) THEN
519  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
520  $ CALL alahd( nout, path )
521  WRITE( nout, fmt = 9996 )trans, n,
522  $ kl, ku, nrhs, imat, k,
523  $ result( k )
524  nfail = nfail + 1
525  END IF
526  60 continue
527  nrun = nrun + 5
528  70 continue
529  80 continue
530 *
531 *+ TEST 7:
532 * Get an estimate of RCOND = 1/CNDNUM.
533 *
534  90 continue
535  DO 100 itran = 1, 2
536  IF( itran.EQ.1 ) THEN
537  anorm = anormo
538  rcondc = rcondo
539  norm = 'O'
540  ELSE
541  anorm = anormi
542  rcondc = rcondi
543  norm = 'I'
544  END IF
545  srnamt = 'DGBCON'
546  CALL dgbcon( norm, n, kl, ku, afac, ldafac,
547  $ iwork, anorm, rcond, work,
548  $ iwork( n+1 ), info )
549 *
550 * Check error code from DGBCON.
551 *
552  IF( info.NE.0 )
553  $ CALL alaerh( path, 'DGBCON', info, 0,
554  $ norm, n, n, kl, ku, -1, imat,
555  $ nfail, nerrs, nout )
556 *
557  result( 7 ) = dget06( rcond, rcondc )
558 *
559 * Print information about the tests that did
560 * not pass the threshold.
561 *
562  IF( result( 7 ).GE.thresh ) THEN
563  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
564  $ CALL alahd( nout, path )
565  WRITE( nout, fmt = 9995 )norm, n, kl, ku,
566  $ imat, 7, result( 7 )
567  nfail = nfail + 1
568  END IF
569  nrun = nrun + 1
570  100 continue
571 *
572  110 continue
573  120 continue
574  130 continue
575  140 continue
576  150 continue
577  160 continue
578 *
579 * Print a summary of the results.
580 *
581  CALL alasum( path, nout, nfail, nrun, nerrs )
582 *
583  9999 format( ' *** In DCHKGB, LA=', i5, ' is too small for M=', i5,
584  $ ', N=', i5, ', KL=', i4, ', KU=', i4,
585  $ / ' ==> Increase LA to at least ', i5 )
586  9998 format( ' *** In DCHKGB, LAFAC=', i5, ' is too small for M=', i5,
587  $ ', N=', i5, ', KL=', i4, ', KU=', i4,
588  $ / ' ==> Increase LAFAC to at least ', i5 )
589  9997 format( ' M =', i5, ', N =', i5, ', KL=', i5, ', KU=', i5,
590  $ ', NB =', i4, ', type ', i1, ', test(', i1, ')=', g12.5 )
591  9996 format( ' TRANS=''', a1, ''', N=', i5, ', KL=', i5, ', KU=', i5,
592  $ ', NRHS=', i3, ', type ', i1, ', test(', i1, ')=', g12.5 )
593  9995 format( ' NORM =''', a1, ''', N=', i5, ', KL=', i5, ', KU=', i5,
594  $ ',', 10x, ' type ', i1, ', test(', i1, ')=', g12.5 )
595 *
596  return
597 *
598 * End of DCHKGB
599 *
600  END