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