LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
schkpt.f
Go to the documentation of this file.
1 *> \brief \b SCHKPT
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 SCHKPT( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
12 * A, D, E, B, X, XACT, WORK, RWORK, NOUT )
13 *
14 * .. Scalar Arguments ..
15 * LOGICAL TSTERR
16 * INTEGER NN, NNS, NOUT
17 * REAL THRESH
18 * ..
19 * .. Array Arguments ..
20 * LOGICAL DOTYPE( * )
21 * INTEGER NSVAL( * ), NVAL( * )
22 * REAL A( * ), B( * ), D( * ), E( * ), RWORK( * ),
23 * $ WORK( * ), X( * ), XACT( * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> SCHKPT tests SPTTRF, -TRS, -RFS, and -CON
33 *> \endverbatim
34 *
35 * Arguments:
36 * ==========
37 *
38 *> \param[in] DOTYPE
39 *> \verbatim
40 *> DOTYPE is LOGICAL array, dimension (NTYPES)
41 *> The matrix types to be used for testing. Matrices of type j
42 *> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
43 *> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
44 *> \endverbatim
45 *>
46 *> \param[in] NN
47 *> \verbatim
48 *> NN is INTEGER
49 *> The number of values of N contained in the vector NVAL.
50 *> \endverbatim
51 *>
52 *> \param[in] NVAL
53 *> \verbatim
54 *> NVAL is INTEGER array, dimension (NN)
55 *> The values of the matrix dimension N.
56 *> \endverbatim
57 *>
58 *> \param[in] NNS
59 *> \verbatim
60 *> NNS is INTEGER
61 *> The number of values of NRHS contained in the vector NSVAL.
62 *> \endverbatim
63 *>
64 *> \param[in] NSVAL
65 *> \verbatim
66 *> NSVAL is INTEGER array, dimension (NNS)
67 *> The values of the number of right hand sides NRHS.
68 *> \endverbatim
69 *>
70 *> \param[in] THRESH
71 *> \verbatim
72 *> THRESH is REAL
73 *> The threshold value for the test ratios. A result is
74 *> included in the output file if RESULT >= THRESH. To have
75 *> every test ratio printed, use THRESH = 0.
76 *> \endverbatim
77 *>
78 *> \param[in] TSTERR
79 *> \verbatim
80 *> TSTERR is LOGICAL
81 *> Flag that indicates whether error exits are to be tested.
82 *> \endverbatim
83 *>
84 *> \param[out] A
85 *> \verbatim
86 *> A is REAL array, dimension (NMAX*2)
87 *> \endverbatim
88 *>
89 *> \param[out] D
90 *> \verbatim
91 *> D is REAL array, dimension (NMAX*2)
92 *> \endverbatim
93 *>
94 *> \param[out] E
95 *> \verbatim
96 *> E is REAL array, dimension (NMAX*2)
97 *> \endverbatim
98 *>
99 *> \param[out] B
100 *> \verbatim
101 *> B is REAL array, dimension (NMAX*NSMAX)
102 *> where NSMAX is the largest entry in NSVAL.
103 *> \endverbatim
104 *>
105 *> \param[out] X
106 *> \verbatim
107 *> X is REAL array, dimension (NMAX*NSMAX)
108 *> \endverbatim
109 *>
110 *> \param[out] XACT
111 *> \verbatim
112 *> XACT is REAL array, dimension (NMAX*NSMAX)
113 *> \endverbatim
114 *>
115 *> \param[out] WORK
116 *> \verbatim
117 *> WORK is REAL array, dimension
118 *> (NMAX*max(3,NSMAX))
119 *> \endverbatim
120 *>
121 *> \param[out] RWORK
122 *> \verbatim
123 *> RWORK is REAL array, dimension
124 *> (max(NMAX,2*NSMAX))
125 *> \endverbatim
126 *>
127 *> \param[in] NOUT
128 *> \verbatim
129 *> NOUT is INTEGER
130 *> The unit number for output.
131 *> \endverbatim
132 *
133 * Authors:
134 * ========
135 *
136 *> \author Univ. of Tennessee
137 *> \author Univ. of California Berkeley
138 *> \author Univ. of Colorado Denver
139 *> \author NAG Ltd.
140 *
141 *> \ingroup single_lin
142 *
143 * =====================================================================
144  SUBROUTINE schkpt( DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR,
145  $ A, D, E, B, X, XACT, WORK, RWORK, NOUT )
146 *
147 * -- LAPACK test routine --
148 * -- LAPACK is a software package provided by Univ. of Tennessee, --
149 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150 *
151 * .. Scalar Arguments ..
152  LOGICAL TSTERR
153  INTEGER NN, NNS, NOUT
154  REAL THRESH
155 * ..
156 * .. Array Arguments ..
157  LOGICAL DOTYPE( * )
158  INTEGER NSVAL( * ), NVAL( * )
159  REAL A( * ), B( * ), D( * ), E( * ), RWORK( * ),
160  $ work( * ), x( * ), xact( * )
161 * ..
162 *
163 * =====================================================================
164 *
165 * .. Parameters ..
166  REAL ONE, ZERO
167  parameter( one = 1.0e+0, zero = 0.0e+0 )
168  INTEGER NTYPES
169  parameter( ntypes = 12 )
170  INTEGER NTESTS
171  parameter( ntests = 7 )
172 * ..
173 * .. Local Scalars ..
174  LOGICAL ZEROT
175  CHARACTER DIST, TYPE
176  CHARACTER*3 PATH
177  INTEGER I, IA, IMAT, IN, INFO, IRHS, IX, IZERO, J, K,
178  $ kl, ku, lda, mode, n, nerrs, nfail, nimat,
179  $ nrhs, nrun
180  REAL AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
181 * ..
182 * .. Local Arrays ..
183  INTEGER ISEED( 4 ), ISEEDY( 4 )
184  REAL RESULT( NTESTS ), Z( 3 )
185 * ..
186 * .. External Functions ..
187  INTEGER ISAMAX
188  REAL SASUM, SGET06, SLANST
189  EXTERNAL isamax, sasum, sget06, slanst
190 * ..
191 * .. External Subroutines ..
192  EXTERNAL alaerh, alahd, alasum, scopy, serrgt, sget04,
195  $ sscal
196 * ..
197 * .. Intrinsic Functions ..
198  INTRINSIC abs, max
199 * ..
200 * .. Scalars in Common ..
201  LOGICAL LERR, OK
202  CHARACTER*32 SRNAMT
203  INTEGER INFOT, NUNIT
204 * ..
205 * .. Common blocks ..
206  COMMON / infoc / infot, nunit, ok, lerr
207  COMMON / srnamc / srnamt
208 * ..
209 * .. Data statements ..
210  DATA iseedy / 0, 0, 0, 1 /
211 * ..
212 * .. Executable Statements ..
213 *
214  path( 1: 1 ) = 'Single precision'
215  path( 2: 3 ) = 'PT'
216  nrun = 0
217  nfail = 0
218  nerrs = 0
219  DO 10 i = 1, 4
220  iseed( i ) = iseedy( i )
221  10 CONTINUE
222 *
223 * Test the error exits
224 *
225  IF( tsterr )
226  $ CALL serrgt( path, nout )
227  infot = 0
228 *
229  DO 110 in = 1, nn
230 *
231 * Do for each value of N in NVAL.
232 *
233  n = nval( in )
234  lda = max( 1, n )
235  nimat = ntypes
236  IF( n.LE.0 )
237  $ nimat = 1
238 *
239  DO 100 imat = 1, nimat
240 *
241 * Do the tests only if DOTYPE( IMAT ) is true.
242 *
243  IF( n.GT.0 .AND. .NOT.dotype( imat ) )
244  $ GO TO 100
245 *
246 * Set up parameters with SLATB4.
247 *
248  CALL slatb4( path, imat, n, n, TYPE, kl, ku, anorm, mode,
249  $ cond, dist )
250 *
251  zerot = imat.GE.8 .AND. imat.LE.10
252  IF( imat.LE.6 ) THEN
253 *
254 * Type 1-6: generate a symmetric tridiagonal matrix of
255 * known condition number in lower triangular band storage.
256 *
257  srnamt = 'SLATMS'
258  CALL slatms( n, n, dist, iseed, TYPE, rwork, mode, cond,
259  $ anorm, kl, ku, 'B', a, 2, work, info )
260 *
261 * Check the error code from SLATMS.
262 *
263  IF( info.NE.0 ) THEN
264  CALL alaerh( path, 'SLATMS', info, 0, ' ', n, n, kl,
265  $ ku, -1, imat, nfail, nerrs, nout )
266  GO TO 100
267  END IF
268  izero = 0
269 *
270 * Copy the matrix to D and E.
271 *
272  ia = 1
273  DO 20 i = 1, n - 1
274  d( i ) = a( ia )
275  e( i ) = a( ia+1 )
276  ia = ia + 2
277  20 CONTINUE
278  IF( n.GT.0 )
279  $ d( n ) = a( ia )
280  ELSE
281 *
282 * Type 7-12: generate a diagonally dominant matrix with
283 * unknown condition number in the vectors D and E.
284 *
285  IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
286 *
287 * Let D and E have values from [-1,1].
288 *
289  CALL slarnv( 2, iseed, n, d )
290  CALL slarnv( 2, iseed, n-1, e )
291 *
292 * Make the tridiagonal matrix diagonally dominant.
293 *
294  IF( n.EQ.1 ) THEN
295  d( 1 ) = abs( d( 1 ) )
296  ELSE
297  d( 1 ) = abs( d( 1 ) ) + abs( e( 1 ) )
298  d( n ) = abs( d( n ) ) + abs( e( n-1 ) )
299  DO 30 i = 2, n - 1
300  d( i ) = abs( d( i ) ) + abs( e( i ) ) +
301  $ abs( e( i-1 ) )
302  30 CONTINUE
303  END IF
304 *
305 * Scale D and E so the maximum element is ANORM.
306 *
307  ix = isamax( n, d, 1 )
308  dmax = d( ix )
309  CALL sscal( n, anorm / dmax, d, 1 )
310  CALL sscal( n-1, anorm / dmax, e, 1 )
311 *
312  ELSE IF( izero.GT.0 ) THEN
313 *
314 * Reuse the last matrix by copying back the zeroed out
315 * elements.
316 *
317  IF( izero.EQ.1 ) THEN
318  d( 1 ) = z( 2 )
319  IF( n.GT.1 )
320  $ e( 1 ) = z( 3 )
321  ELSE IF( izero.EQ.n ) THEN
322  e( n-1 ) = z( 1 )
323  d( n ) = z( 2 )
324  ELSE
325  e( izero-1 ) = z( 1 )
326  d( izero ) = z( 2 )
327  e( izero ) = z( 3 )
328  END IF
329  END IF
330 *
331 * For types 8-10, set one row and column of the matrix to
332 * zero.
333 *
334  izero = 0
335  IF( imat.EQ.8 ) THEN
336  izero = 1
337  z( 2 ) = d( 1 )
338  d( 1 ) = zero
339  IF( n.GT.1 ) THEN
340  z( 3 ) = e( 1 )
341  e( 1 ) = zero
342  END IF
343  ELSE IF( imat.EQ.9 ) THEN
344  izero = n
345  IF( n.GT.1 ) THEN
346  z( 1 ) = e( n-1 )
347  e( n-1 ) = zero
348  END IF
349  z( 2 ) = d( n )
350  d( n ) = zero
351  ELSE IF( imat.EQ.10 ) THEN
352  izero = ( n+1 ) / 2
353  IF( izero.GT.1 ) THEN
354  z( 1 ) = e( izero-1 )
355  e( izero-1 ) = zero
356  z( 3 ) = e( izero )
357  e( izero ) = zero
358  END IF
359  z( 2 ) = d( izero )
360  d( izero ) = zero
361  END IF
362  END IF
363 *
364  CALL scopy( n, d, 1, d( n+1 ), 1 )
365  IF( n.GT.1 )
366  $ CALL scopy( n-1, e, 1, e( n+1 ), 1 )
367 *
368 *+ TEST 1
369 * Factor A as L*D*L' and compute the ratio
370 * norm(L*D*L' - A) / (n * norm(A) * EPS )
371 *
372  CALL spttrf( n, d( n+1 ), e( n+1 ), info )
373 *
374 * Check error code from SPTTRF.
375 *
376  IF( info.NE.izero ) THEN
377  CALL alaerh( path, 'SPTTRF', info, izero, ' ', n, n, -1,
378  $ -1, -1, imat, nfail, nerrs, nout )
379  GO TO 100
380  END IF
381 *
382  IF( info.GT.0 ) THEN
383  rcondc = zero
384  GO TO 90
385  END IF
386 *
387  CALL sptt01( n, d, e, d( n+1 ), e( n+1 ), work,
388  $ result( 1 ) )
389 *
390 * Print the test ratio if greater than or equal to THRESH.
391 *
392  IF( result( 1 ).GE.thresh ) THEN
393  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
394  $ CALL alahd( nout, path )
395  WRITE( nout, fmt = 9999 )n, imat, 1, result( 1 )
396  nfail = nfail + 1
397  END IF
398  nrun = nrun + 1
399 *
400 * Compute RCONDC = 1 / (norm(A) * norm(inv(A))
401 *
402 * Compute norm(A).
403 *
404  anorm = slanst( '1', n, d, e )
405 *
406 * Use SPTTRS to solve for one column at a time of inv(A),
407 * computing the maximum column sum as we go.
408 *
409  ainvnm = zero
410  DO 50 i = 1, n
411  DO 40 j = 1, n
412  x( j ) = zero
413  40 CONTINUE
414  x( i ) = one
415  CALL spttrs( n, 1, d( n+1 ), e( n+1 ), x, lda, info )
416  ainvnm = max( ainvnm, sasum( n, x, 1 ) )
417  50 CONTINUE
418  rcondc = one / max( one, anorm*ainvnm )
419 *
420  DO 80 irhs = 1, nns
421  nrhs = nsval( irhs )
422 *
423 * Generate NRHS random solution vectors.
424 *
425  ix = 1
426  DO 60 j = 1, nrhs
427  CALL slarnv( 2, iseed, n, xact( ix ) )
428  ix = ix + lda
429  60 CONTINUE
430 *
431 * Set the right hand side.
432 *
433  CALL slaptm( n, nrhs, one, d, e, xact, lda, zero, b,
434  $ lda )
435 *
436 *+ TEST 2
437 * Solve A*x = b and compute the residual.
438 *
439  CALL slacpy( 'Full', n, nrhs, b, lda, x, lda )
440  CALL spttrs( n, nrhs, d( n+1 ), e( n+1 ), x, lda, info )
441 *
442 * Check error code from SPTTRS.
443 *
444  IF( info.NE.0 )
445  $ CALL alaerh( path, 'SPTTRS', info, 0, ' ', n, n, -1,
446  $ -1, nrhs, imat, nfail, nerrs, nout )
447 *
448  CALL slacpy( 'Full', n, nrhs, b, lda, work, lda )
449  CALL sptt02( n, nrhs, d, e, x, lda, work, lda,
450  $ result( 2 ) )
451 *
452 *+ TEST 3
453 * Check solution from generated exact solution.
454 *
455  CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
456  $ result( 3 ) )
457 *
458 *+ TESTS 4, 5, and 6
459 * Use iterative refinement to improve the solution.
460 *
461  srnamt = 'SPTRFS'
462  CALL sptrfs( n, nrhs, d, e, d( n+1 ), e( n+1 ), b, lda,
463  $ x, lda, rwork, rwork( nrhs+1 ), work, info )
464 *
465 * Check error code from SPTRFS.
466 *
467  IF( info.NE.0 )
468  $ CALL alaerh( path, 'SPTRFS', info, 0, ' ', n, n, -1,
469  $ -1, nrhs, imat, nfail, nerrs, nout )
470 *
471  CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
472  $ result( 4 ) )
473  CALL sptt05( n, nrhs, d, e, b, lda, x, lda, xact, lda,
474  $ rwork, rwork( nrhs+1 ), result( 5 ) )
475 *
476 * Print information about the tests that did not pass the
477 * threshold.
478 *
479  DO 70 k = 2, 6
480  IF( result( k ).GE.thresh ) THEN
481  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
482  $ CALL alahd( nout, path )
483  WRITE( nout, fmt = 9998 )n, nrhs, imat, k,
484  $ result( k )
485  nfail = nfail + 1
486  END IF
487  70 CONTINUE
488  nrun = nrun + 5
489  80 CONTINUE
490 *
491 *+ TEST 7
492 * Estimate the reciprocal of the condition number of the
493 * matrix.
494 *
495  90 CONTINUE
496  srnamt = 'SPTCON'
497  CALL sptcon( n, d( n+1 ), e( n+1 ), anorm, rcond, rwork,
498  $ info )
499 *
500 * Check error code from SPTCON.
501 *
502  IF( info.NE.0 )
503  $ CALL alaerh( path, 'SPTCON', info, 0, ' ', n, n, -1, -1,
504  $ -1, imat, nfail, nerrs, nout )
505 *
506  result( 7 ) = sget06( rcond, rcondc )
507 *
508 * Print the test ratio if greater than or equal to THRESH.
509 *
510  IF( result( 7 ).GE.thresh ) THEN
511  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
512  $ CALL alahd( nout, path )
513  WRITE( nout, fmt = 9999 )n, imat, 7, result( 7 )
514  nfail = nfail + 1
515  END IF
516  nrun = nrun + 1
517  100 CONTINUE
518  110 CONTINUE
519 *
520 * Print a summary of the results.
521 *
522  CALL alasum( path, nout, nfail, nrun, nerrs )
523 *
524  9999 FORMAT( ' N =', i5, ', type ', i2, ', test ', i2, ', ratio = ',
525  $ g12.5 )
526  9998 FORMAT( ' N =', i5, ', NRHS=', i3, ', type ', i2, ', test(', i2,
527  $ ') = ', g12.5 )
528  RETURN
529 *
530 * End of SCHKPT
531 *
532  END
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:97
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:73
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:107
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:147
subroutine spttrf(N, D, E, INFO)
SPTTRF
Definition: spttrf.f:91
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:321
subroutine sptrfs(N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, WORK, INFO)
SPTRFS
Definition: sptrfs.f:163
subroutine sptcon(N, D, E, ANORM, RCOND, WORK, INFO)
SPTCON
Definition: sptcon.f:118
subroutine spttrs(N, NRHS, D, E, B, LDB, INFO)
SPTTRS
Definition: spttrs.f:109
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine slatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
SLATB4
Definition: slatb4.f:120
subroutine serrgt(PATH, NUNIT)
SERRGT
Definition: serrgt.f:55
subroutine sptt05(N, NRHS, D, E, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
SPTT05
Definition: sptt05.f:150
subroutine sptt02(N, NRHS, D, E, X, LDX, B, LDB, RESID)
SPTT02
Definition: sptt02.f:104
subroutine sget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
SGET04
Definition: sget04.f:102
subroutine schkpt(DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, A, D, E, B, X, XACT, WORK, RWORK, NOUT)
SCHKPT
Definition: schkpt.f:146
subroutine sptt01(N, D, E, DF, EF, WORK, RESID)
SPTT01
Definition: sptt01.f:91
subroutine slaptm(N, NRHS, ALPHA, D, E, X, LDX, BETA, B, LDB)
SLAPTM
Definition: slaptm.f:116