LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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