LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dchkpt.f
Go to the documentation of this file.
1 *> \brief \b DCHKPT
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 DCHKPT( 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 * DOUBLE PRECISION THRESH
18 * ..
19 * .. Array Arguments ..
20 * LOGICAL DOTYPE( * )
21 * INTEGER NSVAL( * ), NVAL( * )
22 * DOUBLE PRECISION A( * ), B( * ), D( * ), E( * ), RWORK( * ),
23 * $ WORK( * ), X( * ), XACT( * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> DCHKPT tests DPTTRF, -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 DOUBLE PRECISION
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 DOUBLE PRECISION array, dimension (NMAX*2)
87 *> \endverbatim
88 *>
89 *> \param[out] D
90 *> \verbatim
91 *> D is DOUBLE PRECISION array, dimension (NMAX*2)
92 *> \endverbatim
93 *>
94 *> \param[out] E
95 *> \verbatim
96 *> E is DOUBLE PRECISION array, dimension (NMAX*2)
97 *> \endverbatim
98 *>
99 *> \param[out] B
100 *> \verbatim
101 *> B is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (NMAX*NSMAX)
108 *> \endverbatim
109 *>
110 *> \param[out] XACT
111 *> \verbatim
112 *> XACT is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
113 *> \endverbatim
114 *>
115 *> \param[out] WORK
116 *> \verbatim
117 *> WORK is DOUBLE PRECISION array, dimension
118 *> (NMAX*max(3,NSMAX))
119 *> \endverbatim
120 *>
121 *> \param[out] RWORK
122 *> \verbatim
123 *> RWORK is DOUBLE PRECISION 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 double_lin
142 *
143 * =====================================================================
144  SUBROUTINE dchkpt( 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  DOUBLE PRECISION THRESH
155 * ..
156 * .. Array Arguments ..
157  LOGICAL DOTYPE( * )
158  INTEGER NSVAL( * ), NVAL( * )
159  DOUBLE PRECISION A( * ), B( * ), D( * ), E( * ), RWORK( * ),
160  $ work( * ), x( * ), xact( * )
161 * ..
162 *
163 * =====================================================================
164 *
165 * .. Parameters ..
166  DOUBLE PRECISION ONE, ZERO
167  parameter( one = 1.0d+0, zero = 0.0d+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  DOUBLE PRECISION AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
181 * ..
182 * .. Local Arrays ..
183  INTEGER ISEED( 4 ), ISEEDY( 4 )
184  DOUBLE PRECISION RESULT( NTESTS ), Z( 3 )
185 * ..
186 * .. External Functions ..
187  INTEGER IDAMAX
188  DOUBLE PRECISION DASUM, DGET06, DLANST
189  EXTERNAL idamax, dasum, dget06, dlanst
190 * ..
191 * .. External Subroutines ..
192  EXTERNAL alaerh, alahd, alasum, dcopy, derrgt, dget04,
195  $ dscal
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 ) = 'Double 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 derrgt( 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 DLATB4.
247 *
248  CALL dlatb4( 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 = 'DLATMS'
258  CALL dlatms( n, n, dist, iseed, TYPE, rwork, mode, cond,
259  $ anorm, kl, ku, 'B', a, 2, work, info )
260 *
261 * Check the error code from DLATMS.
262 *
263  IF( info.NE.0 ) THEN
264  CALL alaerh( path, 'DLATMS', 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 dlarnv( 2, iseed, n, d )
290  CALL dlarnv( 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 = idamax( n, d, 1 )
308  dmax = d( ix )
309  CALL dscal( n, anorm / dmax, d, 1 )
310  CALL dscal( 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 dcopy( n, d, 1, d( n+1 ), 1 )
365  IF( n.GT.1 )
366  $ CALL dcopy( 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 dpttrf( n, d( n+1 ), e( n+1 ), info )
373 *
374 * Check error code from DPTTRF.
375 *
376  IF( info.NE.izero ) THEN
377  CALL alaerh( path, 'DPTTRF', 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 dptt01( 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 = dlanst( '1', n, d, e )
405 *
406 * Use DPTTRS 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 dpttrs( n, 1, d( n+1 ), e( n+1 ), x, lda, info )
416  ainvnm = max( ainvnm, dasum( 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 dlarnv( 2, iseed, n, xact( ix ) )
428  ix = ix + lda
429  60 CONTINUE
430 *
431 * Set the right hand side.
432 *
433  CALL dlaptm( 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 dlacpy( 'Full', n, nrhs, b, lda, x, lda )
440  CALL dpttrs( n, nrhs, d( n+1 ), e( n+1 ), x, lda, info )
441 *
442 * Check error code from DPTTRS.
443 *
444  IF( info.NE.0 )
445  $ CALL alaerh( path, 'DPTTRS', info, 0, ' ', n, n, -1,
446  $ -1, nrhs, imat, nfail, nerrs, nout )
447 *
448  CALL dlacpy( 'Full', n, nrhs, b, lda, work, lda )
449  CALL dptt02( 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 dget04( 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 = 'DPTRFS'
462  CALL dptrfs( 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 DPTRFS.
466 *
467  IF( info.NE.0 )
468  $ CALL alaerh( path, 'DPTRFS', info, 0, ' ', n, n, -1,
469  $ -1, nrhs, imat, nfail, nerrs, nout )
470 *
471  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
472  $ result( 4 ) )
473  CALL dptt05( 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 = 'DPTCON'
497  CALL dptcon( n, d( n+1 ), e( n+1 ), anorm, rcond, rwork,
498  $ info )
499 *
500 * Check error code from DPTCON.
501 *
502  IF( info.NE.0 )
503  $ CALL alaerh( path, 'DPTCON', info, 0, ' ', n, n, -1, -1,
504  $ -1, imat, nfail, nerrs, nout )
505 *
506  result( 7 ) = dget06( 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 DCHKPT
531 *
532  END
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:97
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.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 dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dchkpt(DOTYPE, NN, NVAL, NNS, NSVAL, THRESH, TSTERR, A, D, E, B, X, XACT, WORK, RWORK, NOUT)
DCHKPT
Definition: dchkpt.f:146
subroutine dlaptm(N, NRHS, ALPHA, D, E, X, LDX, BETA, B, LDB)
DLAPTM
Definition: dlaptm.f:116
subroutine dptt02(N, NRHS, D, E, X, LDX, B, LDB, RESID)
DPTT02
Definition: dptt02.f:104
subroutine dget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
DGET04
Definition: dget04.f:102
subroutine derrgt(PATH, NUNIT)
DERRGT
Definition: derrgt.f:55
subroutine dlatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
DLATB4
Definition: dlatb4.f:120
subroutine dptt05(N, NRHS, D, E, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
DPTT05
Definition: dptt05.f:150
subroutine dptt01(N, D, E, DF, EF, WORK, RESID)
DPTT01
Definition: dptt01.f:91
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:321
subroutine dptrfs(N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR, BERR, WORK, INFO)
DPTRFS
Definition: dptrfs.f:163
subroutine dptcon(N, D, E, ANORM, RCOND, WORK, INFO)
DPTCON
Definition: dptcon.f:118
subroutine dpttrf(N, D, E, INFO)
DPTTRF
Definition: dpttrf.f:91
subroutine dpttrs(N, NRHS, D, E, B, LDB, INFO)
DPTTRS
Definition: dpttrs.f:109