LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
cchktb.f
Go to the documentation of this file.
1 *> \brief \b CCHKTB
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 CCHKTB( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
12 * NMAX, AB, AINV, B, X, XACT, WORK, RWORK, NOUT )
13 *
14 * .. Scalar Arguments ..
15 * LOGICAL TSTERR
16 * INTEGER NMAX, NN, NNS, NOUT
17 * REAL THRESH
18 * ..
19 * .. Array Arguments ..
20 * LOGICAL DOTYPE( * )
21 * INTEGER NSVAL( * ), NVAL( * )
22 * REAL RWORK( * )
23 * COMPLEX AB( * ), AINV( * ), B( * ), WORK( * ), X( * ),
24 * $ XACT( * )
25 * ..
26 *
27 *
28 *> \par Purpose:
29 * =============
30 *>
31 *> \verbatim
32 *>
33 *> CCHKTB tests CTBTRS, -RFS, and -CON, and CLATBS.
34 *> \endverbatim
35 *
36 * Arguments:
37 * ==========
38 *
39 *> \param[in] DOTYPE
40 *> \verbatim
41 *> DOTYPE is LOGICAL array, dimension (NTYPES)
42 *> The matrix types to be used for testing. Matrices of type j
43 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
44 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
45 *> \endverbatim
46 *>
47 *> \param[in] NN
48 *> \verbatim
49 *> NN is INTEGER
50 *> The number of values of N contained in the vector NVAL.
51 *> \endverbatim
52 *>
53 *> \param[in] NVAL
54 *> \verbatim
55 *> NVAL is INTEGER array, dimension (NN)
56 *> The values of the matrix column dimension N.
57 *> \endverbatim
58 *>
59 *> \param[in] NNS
60 *> \verbatim
61 *> NNS is INTEGER
62 *> The number of values of NRHS contained in the vector NSVAL.
63 *> \endverbatim
64 *>
65 *> \param[in] NSVAL
66 *> \verbatim
67 *> NSVAL is INTEGER array, dimension (NNS)
68 *> The values of the number of right hand sides NRHS.
69 *> \endverbatim
70 *>
71 *> \param[in] THRESH
72 *> \verbatim
73 *> THRESH is REAL
74 *> The threshold value for the test ratios. A result is
75 *> included in the output file if RESULT >= THRESH. To have
76 *> every test ratio printed, use THRESH = 0.
77 *> \endverbatim
78 *>
79 *> \param[in] TSTERR
80 *> \verbatim
81 *> TSTERR is LOGICAL
82 *> Flag that indicates whether error exits are to be tested.
83 *> \endverbatim
84 *>
85 *> \param[in] NMAX
86 *> \verbatim
87 *> NMAX is INTEGER
88 *> The leading dimension of the work arrays.
89 *> NMAX >= the maximum value of N in NVAL.
90 *> \endverbatim
91 *>
92 *> \param[out] AB
93 *> \verbatim
94 *> AB is COMPLEX array, dimension (NMAX*NMAX)
95 *> \endverbatim
96 *>
97 *> \param[out] AINV
98 *> \verbatim
99 *> AINV is COMPLEX array, dimension (NMAX*NMAX)
100 *> \endverbatim
101 *>
102 *> \param[out] B
103 *> \verbatim
104 *> B is COMPLEX array, dimension (NMAX*NSMAX)
105 *> where NSMAX is the largest entry in NSVAL.
106 *> \endverbatim
107 *>
108 *> \param[out] X
109 *> \verbatim
110 *> X is COMPLEX array, dimension (NMAX*NSMAX)
111 *> \endverbatim
112 *>
113 *> \param[out] XACT
114 *> \verbatim
115 *> XACT is COMPLEX array, dimension (NMAX*NSMAX)
116 *> \endverbatim
117 *>
118 *> \param[out] WORK
119 *> \verbatim
120 *> WORK is COMPLEX array, dimension
121 *> (NMAX*max(3,NSMAX))
122 *> \endverbatim
123 *>
124 *> \param[out] RWORK
125 *> \verbatim
126 *> RWORK is REAL array, dimension
127 *> (max(NMAX,2*NSMAX))
128 *> \endverbatim
129 *>
130 *> \param[in] NOUT
131 *> \verbatim
132 *> NOUT is INTEGER
133 *> The unit number for output.
134 *> \endverbatim
135 *
136 * Authors:
137 * ========
138 *
139 *> \author Univ. of Tennessee
140 *> \author Univ. of California Berkeley
141 *> \author Univ. of Colorado Denver
142 *> \author NAG Ltd.
143 *
144 *> \date November 2011
145 *
146 *> \ingroup complex_lin
147 *
148 * =====================================================================
149  SUBROUTINE cchktb( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
150  $ nmax, ab, ainv, b, x, xact, work, rwork, nout )
151 *
152 * -- LAPACK test routine (version 3.4.0) --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155 * November 2011
156 *
157 * .. Scalar Arguments ..
158  LOGICAL tsterr
159  INTEGER nmax, nn, nns, nout
160  REAL thresh
161 * ..
162 * .. Array Arguments ..
163  LOGICAL dotype( * )
164  INTEGER nsval( * ), nval( * )
165  REAL rwork( * )
166  COMPLEX ab( * ), ainv( * ), b( * ), work( * ), x( * ),
167  $ xact( * )
168 * ..
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..
173  INTEGER ntype1, ntypes
174  parameter( ntype1 = 9, ntypes = 17 )
175  INTEGER ntests
176  parameter( ntests = 8 )
177  INTEGER ntran
178  parameter( ntran = 3 )
179  REAL one, zero
180  parameter( one = 1.0e+0, zero = 0.0e+0 )
181 * ..
182 * .. Local Scalars ..
183  CHARACTER diag, norm, trans, uplo, xtype
184  CHARACTER*3 path
185  INTEGER i, idiag, ik, imat, in, info, irhs, itran,
186  $ iuplo, j, k, kd, lda, ldab, n, nerrs, nfail,
187  $ nimat, nimat2, nk, nrhs, nrun
188  REAL ainvnm, anorm, rcond, rcondc, rcondi, rcondo,
189  $ scale
190 * ..
191 * .. Local Arrays ..
192  CHARACTER transs( ntran ), uplos( 2 )
193  INTEGER iseed( 4 ), iseedy( 4 )
194  REAL result( ntests )
195 * ..
196 * .. External Functions ..
197  LOGICAL lsame
198  REAL clantb, clantr
199  EXTERNAL lsame, clantb, clantr
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL alaerh, alahd, alasum, ccopy, cerrtr, cget04,
203  $ clacpy, clarhs, claset, clatbs, clattb, ctbcon,
204  $ ctbrfs, ctbsv, ctbt02, ctbt03, ctbt05, ctbt06,
205  $ ctbtrs
206 * ..
207 * .. Scalars in Common ..
208  LOGICAL lerr, ok
209  CHARACTER*32 srnamt
210  INTEGER infot, iounit
211 * ..
212 * .. Common blocks ..
213  common / infoc / infot, iounit, ok, lerr
214  common / srnamc / srnamt
215 * ..
216 * .. Intrinsic Functions ..
217  INTRINSIC cmplx, max, min
218 * ..
219 * .. Data statements ..
220  DATA iseedy / 1988, 1989, 1990, 1991 /
221  DATA uplos / 'U', 'L' / , transs / 'N', 'T', 'C' /
222 * ..
223 * .. Executable Statements ..
224 *
225 * Initialize constants and the random number seed.
226 *
227  path( 1: 1 ) = 'Complex precision'
228  path( 2: 3 ) = 'TB'
229  nrun = 0
230  nfail = 0
231  nerrs = 0
232  DO 10 i = 1, 4
233  iseed( i ) = iseedy( i )
234  10 continue
235 *
236 * Test the error exits
237 *
238  IF( tsterr )
239  $ CALL cerrtr( path, nout )
240  infot = 0
241 *
242  DO 140 in = 1, nn
243 *
244 * Do for each value of N in NVAL
245 *
246  n = nval( in )
247  lda = max( 1, n )
248  xtype = 'N'
249  nimat = ntype1
250  nimat2 = ntypes
251  IF( n.LE.0 ) THEN
252  nimat = 1
253  nimat2 = ntype1 + 1
254  END IF
255 *
256  nk = min( n+1, 4 )
257  DO 130 ik = 1, nk
258 *
259 * Do for KD = 0, N, (3N-1)/4, and (N+1)/4. This order makes
260 * it easier to skip redundant values for small values of N.
261 *
262  IF( ik.EQ.1 ) THEN
263  kd = 0
264  ELSE IF( ik.EQ.2 ) THEN
265  kd = max( n, 0 )
266  ELSE IF( ik.EQ.3 ) THEN
267  kd = ( 3*n-1 ) / 4
268  ELSE IF( ik.EQ.4 ) THEN
269  kd = ( n+1 ) / 4
270  END IF
271  ldab = kd + 1
272 *
273  DO 90 imat = 1, nimat
274 *
275 * Do the tests only if DOTYPE( IMAT ) is true.
276 *
277  IF( .NOT.dotype( imat ) )
278  $ go to 90
279 *
280  DO 80 iuplo = 1, 2
281 *
282 * Do first for UPLO = 'U', then for UPLO = 'L'
283 *
284  uplo = uplos( iuplo )
285 *
286 * Call CLATTB to generate a triangular test matrix.
287 *
288  srnamt = 'CLATTB'
289  CALL clattb( imat, uplo, 'No transpose', diag, iseed,
290  $ n, kd, ab, ldab, x, work, rwork, info )
291 *
292 * Set IDIAG = 1 for non-unit matrices, 2 for unit.
293 *
294  IF( lsame( diag, 'N' ) ) THEN
295  idiag = 1
296  ELSE
297  idiag = 2
298  END IF
299 *
300 * Form the inverse of A so we can get a good estimate
301 * of RCONDC = 1/(norm(A) * norm(inv(A))).
302 *
303  CALL claset( 'Full', n, n, cmplx( zero ),
304  $ cmplx( one ), ainv, lda )
305  IF( lsame( uplo, 'U' ) ) THEN
306  DO 20 j = 1, n
307  CALL ctbsv( uplo, 'No transpose', diag, j, kd,
308  $ ab, ldab, ainv( ( j-1 )*lda+1 ), 1 )
309  20 continue
310  ELSE
311  DO 30 j = 1, n
312  CALL ctbsv( uplo, 'No transpose', diag, n-j+1,
313  $ kd, ab( ( j-1 )*ldab+1 ), ldab,
314  $ ainv( ( j-1 )*lda+j ), 1 )
315  30 continue
316  END IF
317 *
318 * Compute the 1-norm condition number of A.
319 *
320  anorm = clantb( '1', uplo, diag, n, kd, ab, ldab,
321  $ rwork )
322  ainvnm = clantr( '1', uplo, diag, n, n, ainv, lda,
323  $ rwork )
324  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
325  rcondo = one
326  ELSE
327  rcondo = ( one / anorm ) / ainvnm
328  END IF
329 *
330 * Compute the infinity-norm condition number of A.
331 *
332  anorm = clantb( 'I', uplo, diag, n, kd, ab, ldab,
333  $ rwork )
334  ainvnm = clantr( 'I', uplo, diag, n, n, ainv, lda,
335  $ rwork )
336  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
337  rcondi = one
338  ELSE
339  rcondi = ( one / anorm ) / ainvnm
340  END IF
341 *
342  DO 60 irhs = 1, nns
343  nrhs = nsval( irhs )
344  xtype = 'N'
345 *
346  DO 50 itran = 1, ntran
347 *
348 * Do for op(A) = A, A**T, or A**H.
349 *
350  trans = transs( itran )
351  IF( itran.EQ.1 ) THEN
352  norm = 'O'
353  rcondc = rcondo
354  ELSE
355  norm = 'I'
356  rcondc = rcondi
357  END IF
358 *
359 *+ TEST 1
360 * Solve and compute residual for op(A)*x = b.
361 *
362  srnamt = 'CLARHS'
363  CALL clarhs( path, xtype, uplo, trans, n, n, kd,
364  $ idiag, nrhs, ab, ldab, xact, lda,
365  $ b, lda, iseed, info )
366  xtype = 'C'
367  CALL clacpy( 'Full', n, nrhs, b, lda, x, lda )
368 *
369  srnamt = 'CTBTRS'
370  CALL ctbtrs( uplo, trans, diag, n, kd, nrhs, ab,
371  $ ldab, x, lda, info )
372 *
373 * Check error code from CTBTRS.
374 *
375  IF( info.NE.0 )
376  $ CALL alaerh( path, 'CTBTRS', info, 0,
377  $ uplo // trans // diag, n, n, kd,
378  $ kd, nrhs, imat, nfail, nerrs,
379  $ nout )
380 *
381  CALL ctbt02( uplo, trans, diag, n, kd, nrhs, ab,
382  $ ldab, x, lda, b, lda, work, rwork,
383  $ result( 1 ) )
384 *
385 *+ TEST 2
386 * Check solution from generated exact solution.
387 *
388  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
389  $ result( 2 ) )
390 *
391 *+ TESTS 3, 4, and 5
392 * Use iterative refinement to improve the solution
393 * and compute error bounds.
394 *
395  srnamt = 'CTBRFS'
396  CALL ctbrfs( uplo, trans, diag, n, kd, nrhs, ab,
397  $ ldab, b, lda, x, lda, rwork,
398  $ rwork( nrhs+1 ), work,
399  $ rwork( 2*nrhs+1 ), info )
400 *
401 * Check error code from CTBRFS.
402 *
403  IF( info.NE.0 )
404  $ CALL alaerh( path, 'CTBRFS', info, 0,
405  $ uplo // trans // diag, n, n, kd,
406  $ kd, nrhs, imat, nfail, nerrs,
407  $ nout )
408 *
409  CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
410  $ result( 3 ) )
411  CALL ctbt05( uplo, trans, diag, n, kd, nrhs, ab,
412  $ ldab, b, lda, x, lda, xact, lda,
413  $ rwork, rwork( nrhs+1 ),
414  $ result( 4 ) )
415 *
416 * Print information about the tests that did not
417 * pass the threshold.
418 *
419  DO 40 k = 1, 5
420  IF( result( k ).GE.thresh ) THEN
421  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
422  $ CALL alahd( nout, path )
423  WRITE( nout, fmt = 9999 )uplo, trans,
424  $ diag, n, kd, nrhs, imat, k, result( k )
425  nfail = nfail + 1
426  END IF
427  40 continue
428  nrun = nrun + 5
429  50 continue
430  60 continue
431 *
432 *+ TEST 6
433 * Get an estimate of RCOND = 1/CNDNUM.
434 *
435  DO 70 itran = 1, 2
436  IF( itran.EQ.1 ) THEN
437  norm = 'O'
438  rcondc = rcondo
439  ELSE
440  norm = 'I'
441  rcondc = rcondi
442  END IF
443  srnamt = 'CTBCON'
444  CALL ctbcon( norm, uplo, diag, n, kd, ab, ldab,
445  $ rcond, work, rwork, info )
446 *
447 * Check error code from CTBCON.
448 *
449  IF( info.NE.0 )
450  $ CALL alaerh( path, 'CTBCON', info, 0,
451  $ norm // uplo // diag, n, n, kd, kd,
452  $ -1, imat, nfail, nerrs, nout )
453 *
454  CALL ctbt06( rcond, rcondc, uplo, diag, n, kd, ab,
455  $ ldab, rwork, result( 6 ) )
456 *
457 * Print the test ratio if it is .GE. THRESH.
458 *
459  IF( result( 6 ).GE.thresh ) THEN
460  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
461  $ CALL alahd( nout, path )
462  WRITE( nout, fmt = 9998 ) 'CTBCON', norm, uplo,
463  $ diag, n, kd, imat, 6, result( 6 )
464  nfail = nfail + 1
465  END IF
466  nrun = nrun + 1
467  70 continue
468  80 continue
469  90 continue
470 *
471 * Use pathological test matrices to test CLATBS.
472 *
473  DO 120 imat = ntype1 + 1, nimat2
474 *
475 * Do the tests only if DOTYPE( IMAT ) is true.
476 *
477  IF( .NOT.dotype( imat ) )
478  $ go to 120
479 *
480  DO 110 iuplo = 1, 2
481 *
482 * Do first for UPLO = 'U', then for UPLO = 'L'
483 *
484  uplo = uplos( iuplo )
485  DO 100 itran = 1, ntran
486 *
487 * Do for op(A) = A, A**T, and A**H.
488 *
489  trans = transs( itran )
490 *
491 * Call CLATTB to generate a triangular test matrix.
492 *
493  srnamt = 'CLATTB'
494  CALL clattb( imat, uplo, trans, diag, iseed, n, kd,
495  $ ab, ldab, x, work, rwork, info )
496 *
497 *+ TEST 7
498 * Solve the system op(A)*x = b
499 *
500  srnamt = 'CLATBS'
501  CALL ccopy( n, x, 1, b, 1 )
502  CALL clatbs( uplo, trans, diag, 'N', n, kd, ab,
503  $ ldab, b, scale, rwork, info )
504 *
505 * Check error code from CLATBS.
506 *
507  IF( info.NE.0 )
508  $ CALL alaerh( path, 'CLATBS', info, 0,
509  $ uplo // trans // diag // 'N', n, n,
510  $ kd, kd, -1, imat, nfail, nerrs,
511  $ nout )
512 *
513  CALL ctbt03( uplo, trans, diag, n, kd, 1, ab, ldab,
514  $ scale, rwork, one, b, lda, x, lda,
515  $ work, result( 7 ) )
516 *
517 *+ TEST 8
518 * Solve op(A)*x = b again with NORMIN = 'Y'.
519 *
520  CALL ccopy( n, x, 1, b, 1 )
521  CALL clatbs( uplo, trans, diag, 'Y', n, kd, ab,
522  $ ldab, b, scale, rwork, info )
523 *
524 * Check error code from CLATBS.
525 *
526  IF( info.NE.0 )
527  $ CALL alaerh( path, 'CLATBS', info, 0,
528  $ uplo // trans // diag // 'Y', n, n,
529  $ kd, kd, -1, imat, nfail, nerrs,
530  $ nout )
531 *
532  CALL ctbt03( uplo, trans, diag, n, kd, 1, ab, ldab,
533  $ scale, rwork, one, b, lda, x, lda,
534  $ work, result( 8 ) )
535 *
536 * Print information about the tests that did not pass
537 * the threshold.
538 *
539  IF( result( 7 ).GE.thresh ) THEN
540  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
541  $ CALL alahd( nout, path )
542  WRITE( nout, fmt = 9997 )'CLATBS', uplo, trans,
543  $ diag, 'N', n, kd, imat, 7, result( 7 )
544  nfail = nfail + 1
545  END IF
546  IF( result( 8 ).GE.thresh ) THEN
547  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
548  $ CALL alahd( nout, path )
549  WRITE( nout, fmt = 9997 )'CLATBS', uplo, trans,
550  $ diag, 'Y', n, kd, imat, 8, result( 8 )
551  nfail = nfail + 1
552  END IF
553  nrun = nrun + 2
554  100 continue
555  110 continue
556  120 continue
557  130 continue
558  140 continue
559 *
560 * Print a summary of the results.
561 *
562  CALL alasum( path, nout, nfail, nrun, nerrs )
563 *
564  9999 format( ' UPLO=''', a1, ''', TRANS=''', a1, ''
565 ', $ DIAG=''', a1, ''', N=', i5, ', KD=', i5, ', NRHS=', i5,
566  $ ', type ', i2, ', test(', i2, ')=', g12.5 )
567  9998 format( 1x, a, '( ''', a1, ''', ''', a1, ''', ''', a1, ''',',
568  $ i5, ',', i5, ', ... ), type ', i2, ', test(', i2, ')=',
569  $ g12.5 )
570  9997 format( 1x, a, '( ''', a1, ''', ''', a1, ''', ''', a1, ''', ''',
571  $ a1, ''',', i5, ',', i5, ', ... ), type ', i2, ', test(',
572  $ i1, ')=', g12.5 )
573  return
574 *
575 * End of CCHKTB
576 *
577  END