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