LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ddrvpt.f
Go to the documentation of this file.
1*> \brief \b DDRVPT
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 DDRVPT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D,
12* E, B, X, XACT, WORK, RWORK, NOUT )
13*
14* .. Scalar Arguments ..
15* LOGICAL TSTERR
16* INTEGER NN, NOUT, NRHS
17* DOUBLE PRECISION THRESH
18* ..
19* .. Array Arguments ..
20* LOGICAL DOTYPE( * )
21* INTEGER 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*> DDRVPT tests DPTSV and -SVX.
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] NRHS
59*> \verbatim
60*> NRHS is INTEGER
61*> The number of right hand side vectors to be generated for
62*> each linear system.
63*> \endverbatim
64*>
65*> \param[in] THRESH
66*> \verbatim
67*> THRESH is DOUBLE PRECISION
68*> The threshold value for the test ratios. A result is
69*> included in the output file if RESULT >= THRESH. To have
70*> every test ratio printed, use THRESH = 0.
71*> \endverbatim
72*>
73*> \param[in] TSTERR
74*> \verbatim
75*> TSTERR is LOGICAL
76*> Flag that indicates whether error exits are to be tested.
77*> \endverbatim
78*>
79*> \param[out] A
80*> \verbatim
81*> A is DOUBLE PRECISION array, dimension (NMAX*2)
82*> \endverbatim
83*>
84*> \param[out] D
85*> \verbatim
86*> D is DOUBLE PRECISION array, dimension (NMAX*2)
87*> \endverbatim
88*>
89*> \param[out] E
90*> \verbatim
91*> E is DOUBLE PRECISION array, dimension (NMAX*2)
92*> \endverbatim
93*>
94*> \param[out] B
95*> \verbatim
96*> B is DOUBLE PRECISION array, dimension (NMAX*NRHS)
97*> \endverbatim
98*>
99*> \param[out] X
100*> \verbatim
101*> X is DOUBLE PRECISION array, dimension (NMAX*NRHS)
102*> \endverbatim
103*>
104*> \param[out] XACT
105*> \verbatim
106*> XACT is DOUBLE PRECISION array, dimension (NMAX*NRHS)
107*> \endverbatim
108*>
109*> \param[out] WORK
110*> \verbatim
111*> WORK is DOUBLE PRECISION array, dimension
112*> (NMAX*max(3,NRHS))
113*> \endverbatim
114*>
115*> \param[out] RWORK
116*> \verbatim
117*> RWORK is DOUBLE PRECISION array, dimension
118*> (max(NMAX,2*NRHS))
119*> \endverbatim
120*>
121*> \param[in] NOUT
122*> \verbatim
123*> NOUT is INTEGER
124*> The unit number for output.
125*> \endverbatim
126*
127* Authors:
128* ========
129*
130*> \author Univ. of Tennessee
131*> \author Univ. of California Berkeley
132*> \author Univ. of Colorado Denver
133*> \author NAG Ltd.
134*
135*> \ingroup double_lin
136*
137* =====================================================================
138 SUBROUTINE ddrvpt( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D,
139 $ E, B, X, XACT, WORK, RWORK, NOUT )
140*
141* -- LAPACK test routine --
142* -- LAPACK is a software package provided by Univ. of Tennessee, --
143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144*
145* .. Scalar Arguments ..
146 LOGICAL TSTERR
147 INTEGER NN, NOUT, NRHS
148 DOUBLE PRECISION THRESH
149* ..
150* .. Array Arguments ..
151 LOGICAL DOTYPE( * )
152 INTEGER NVAL( * )
153 DOUBLE PRECISION A( * ), B( * ), D( * ), E( * ), RWORK( * ),
154 $ work( * ), x( * ), xact( * )
155* ..
156*
157* =====================================================================
158*
159* .. Parameters ..
160 DOUBLE PRECISION ONE, ZERO
161 parameter( one = 1.0d+0, zero = 0.0d+0 )
162 INTEGER NTYPES
163 parameter( ntypes = 12 )
164 INTEGER NTESTS
165 parameter( ntests = 6 )
166* ..
167* .. Local Scalars ..
168 LOGICAL ZEROT
169 CHARACTER DIST, FACT, TYPE
170 CHARACTER*3 PATH
171 INTEGER I, IA, IFACT, IMAT, IN, INFO, IX, IZERO, J, K,
172 $ k1, kl, ku, lda, mode, n, nerrs, nfail, nimat,
173 $ nrun, nt
174 DOUBLE PRECISION AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
175* ..
176* .. Local Arrays ..
177 INTEGER ISEED( 4 ), ISEEDY( 4 )
178 DOUBLE PRECISION RESULT( NTESTS ), Z( 3 )
179* ..
180* .. External Functions ..
181 INTEGER IDAMAX
182 DOUBLE PRECISION DASUM, DGET06, DLANST
183 EXTERNAL idamax, dasum, dget06, dlanst
184* ..
185* .. External Subroutines ..
186 EXTERNAL aladhd, alaerh, alasvm, dcopy, derrvx, dget04,
189 $ dpttrs, dscal
190* ..
191* .. Intrinsic Functions ..
192 INTRINSIC abs, max
193* ..
194* .. Scalars in Common ..
195 LOGICAL LERR, OK
196 CHARACTER*32 SRNAMT
197 INTEGER INFOT, NUNIT
198* ..
199* .. Common blocks ..
200 COMMON / infoc / infot, nunit, ok, lerr
201 COMMON / srnamc / srnamt
202* ..
203* .. Data statements ..
204 DATA iseedy / 0, 0, 0, 1 /
205* ..
206* .. Executable Statements ..
207*
208 path( 1: 1 ) = 'Double precision'
209 path( 2: 3 ) = 'PT'
210 nrun = 0
211 nfail = 0
212 nerrs = 0
213 DO 10 i = 1, 4
214 iseed( i ) = iseedy( i )
215 10 CONTINUE
216*
217* Test the error exits
218*
219 IF( tsterr )
220 $ CALL derrvx( path, nout )
221 infot = 0
222*
223 DO 120 in = 1, nn
224*
225* Do for each value of N in NVAL.
226*
227 n = nval( in )
228 lda = max( 1, n )
229 nimat = ntypes
230 IF( n.LE.0 )
231 $ nimat = 1
232*
233 DO 110 imat = 1, nimat
234*
235* Do the tests only if DOTYPE( IMAT ) is true.
236*
237 IF( n.GT.0 .AND. .NOT.dotype( imat ) )
238 $ GO TO 110
239*
240* Set up parameters with DLATB4.
241*
242 CALL dlatb4( path, imat, n, n, TYPE, kl, ku, anorm, mode,
243 $ cond, dist )
244*
245 zerot = imat.GE.8 .AND. imat.LE.10
246 IF( imat.LE.6 ) THEN
247*
248* Type 1-6: generate a symmetric tridiagonal matrix of
249* known condition number in lower triangular band storage.
250*
251 srnamt = 'DLATMS'
252 CALL dlatms( n, n, dist, iseed, TYPE, rwork, mode, cond,
253 $ anorm, kl, ku, 'B', a, 2, work, info )
254*
255* Check the error code from DLATMS.
256*
257 IF( info.NE.0 ) THEN
258 CALL alaerh( path, 'DLATMS', info, 0, ' ', n, n, kl,
259 $ ku, -1, imat, nfail, nerrs, nout )
260 GO TO 110
261 END IF
262 izero = 0
263*
264* Copy the matrix to D and E.
265*
266 ia = 1
267 DO 20 i = 1, n - 1
268 d( i ) = a( ia )
269 e( i ) = a( ia+1 )
270 ia = ia + 2
271 20 CONTINUE
272 IF( n.GT.0 )
273 $ d( n ) = a( ia )
274 ELSE
275*
276* Type 7-12: generate a diagonally dominant matrix with
277* unknown condition number in the vectors D and E.
278*
279 IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
280*
281* Let D and E have values from [-1,1].
282*
283 CALL dlarnv( 2, iseed, n, d )
284 CALL dlarnv( 2, iseed, n-1, e )
285*
286* Make the tridiagonal matrix diagonally dominant.
287*
288 IF( n.EQ.1 ) THEN
289 d( 1 ) = abs( d( 1 ) )
290 ELSE
291 d( 1 ) = abs( d( 1 ) ) + abs( e( 1 ) )
292 d( n ) = abs( d( n ) ) + abs( e( n-1 ) )
293 DO 30 i = 2, n - 1
294 d( i ) = abs( d( i ) ) + abs( e( i ) ) +
295 $ abs( e( i-1 ) )
296 30 CONTINUE
297 END IF
298*
299* Scale D and E so the maximum element is ANORM.
300*
301 ix = idamax( n, d, 1 )
302 dmax = d( ix )
303 CALL dscal( n, anorm / dmax, d, 1 )
304 IF( n.GT.1 )
305 $ CALL dscal( n-1, anorm / dmax, e, 1 )
306*
307 ELSE IF( izero.GT.0 ) THEN
308*
309* Reuse the last matrix by copying back the zeroed out
310* elements.
311*
312 IF( izero.EQ.1 ) THEN
313 d( 1 ) = z( 2 )
314 IF( n.GT.1 )
315 $ e( 1 ) = z( 3 )
316 ELSE IF( izero.EQ.n ) THEN
317 e( n-1 ) = z( 1 )
318 d( n ) = z( 2 )
319 ELSE
320 e( izero-1 ) = z( 1 )
321 d( izero ) = z( 2 )
322 e( izero ) = z( 3 )
323 END IF
324 END IF
325*
326* For types 8-10, set one row and column of the matrix to
327* zero.
328*
329 izero = 0
330 IF( imat.EQ.8 ) THEN
331 izero = 1
332 z( 2 ) = d( 1 )
333 d( 1 ) = zero
334 IF( n.GT.1 ) THEN
335 z( 3 ) = e( 1 )
336 e( 1 ) = zero
337 END IF
338 ELSE IF( imat.EQ.9 ) THEN
339 izero = n
340 IF( n.GT.1 ) THEN
341 z( 1 ) = e( n-1 )
342 e( n-1 ) = zero
343 END IF
344 z( 2 ) = d( n )
345 d( n ) = zero
346 ELSE IF( imat.EQ.10 ) THEN
347 izero = ( n+1 ) / 2
348 IF( izero.GT.1 ) THEN
349 z( 1 ) = e( izero-1 )
350 z( 3 ) = e( izero )
351 e( izero-1 ) = zero
352 e( izero ) = zero
353 END IF
354 z( 2 ) = d( izero )
355 d( izero ) = zero
356 END IF
357 END IF
358*
359* Generate NRHS random solution vectors.
360*
361 ix = 1
362 DO 40 j = 1, nrhs
363 CALL dlarnv( 2, iseed, n, xact( ix ) )
364 ix = ix + lda
365 40 CONTINUE
366*
367* Set the right hand side.
368*
369 CALL dlaptm( n, nrhs, one, d, e, xact, lda, zero, b, lda )
370*
371 DO 100 ifact = 1, 2
372 IF( ifact.EQ.1 ) THEN
373 fact = 'F'
374 ELSE
375 fact = 'N'
376 END IF
377*
378* Compute the condition number for comparison with
379* the value returned by DPTSVX.
380*
381 IF( zerot ) THEN
382 IF( ifact.EQ.1 )
383 $ GO TO 100
384 rcondc = zero
385*
386 ELSE IF( ifact.EQ.1 ) THEN
387*
388* Compute the 1-norm of A.
389*
390 anorm = dlanst( '1', n, d, e )
391*
392 CALL dcopy( n, d, 1, d( n+1 ), 1 )
393 IF( n.GT.1 )
394 $ CALL dcopy( n-1, e, 1, e( n+1 ), 1 )
395*
396* Factor the matrix A.
397*
398 CALL dpttrf( n, d( n+1 ), e( n+1 ), info )
399*
400* Use DPTTRS to solve for one column at a time of
401* inv(A), computing the maximum column sum as we go.
402*
403 ainvnm = zero
404 DO 60 i = 1, n
405 DO 50 j = 1, n
406 x( j ) = zero
407 50 CONTINUE
408 x( i ) = one
409 CALL dpttrs( n, 1, d( n+1 ), e( n+1 ), x, lda,
410 $ info )
411 ainvnm = max( ainvnm, dasum( n, x, 1 ) )
412 60 CONTINUE
413*
414* Compute the 1-norm condition number of A.
415*
416 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
417 rcondc = one
418 ELSE
419 rcondc = ( one / anorm ) / ainvnm
420 END IF
421 END IF
422*
423 IF( ifact.EQ.2 ) THEN
424*
425* --- Test DPTSV --
426*
427 CALL dcopy( n, d, 1, d( n+1 ), 1 )
428 IF( n.GT.1 )
429 $ CALL dcopy( n-1, e, 1, e( n+1 ), 1 )
430 CALL dlacpy( 'Full', n, nrhs, b, lda, x, lda )
431*
432* Factor A as L*D*L' and solve the system A*X = B.
433*
434 srnamt = 'DPTSV '
435 CALL dptsv( n, nrhs, d( n+1 ), e( n+1 ), x, lda,
436 $ info )
437*
438* Check error code from DPTSV .
439*
440 IF( info.NE.izero )
441 $ CALL alaerh( path, 'DPTSV ', info, izero, ' ', n,
442 $ n, 1, 1, nrhs, imat, nfail, nerrs,
443 $ nout )
444 nt = 0
445 IF( izero.EQ.0 ) THEN
446*
447* Check the factorization by computing the ratio
448* norm(L*D*L' - A) / (n * norm(A) * EPS )
449*
450 CALL dptt01( n, d, e, d( n+1 ), e( n+1 ), work,
451 $ result( 1 ) )
452*
453* Compute the residual in the solution.
454*
455 CALL dlacpy( 'Full', n, nrhs, b, lda, work, lda )
456 CALL dptt02( n, nrhs, d, e, x, lda, work, lda,
457 $ result( 2 ) )
458*
459* Check solution from generated exact solution.
460*
461 CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
462 $ result( 3 ) )
463 nt = 3
464 END IF
465*
466* Print information about the tests that did not pass
467* the threshold.
468*
469 DO 70 k = 1, nt
470 IF( result( k ).GE.thresh ) THEN
471 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
472 $ CALL aladhd( nout, path )
473 WRITE( nout, fmt = 9999 )'DPTSV ', n, imat, k,
474 $ result( k )
475 nfail = nfail + 1
476 END IF
477 70 CONTINUE
478 nrun = nrun + nt
479 END IF
480*
481* --- Test DPTSVX ---
482*
483 IF( ifact.GT.1 ) THEN
484*
485* Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
486*
487 DO 80 i = 1, n - 1
488 d( n+i ) = zero
489 e( n+i ) = zero
490 80 CONTINUE
491 IF( n.GT.0 )
492 $ d( n+n ) = zero
493 END IF
494*
495 CALL dlaset( 'Full', n, nrhs, zero, zero, x, lda )
496*
497* Solve the system and compute the condition number and
498* error bounds using DPTSVX.
499*
500 srnamt = 'DPTSVX'
501 CALL dptsvx( fact, n, nrhs, d, e, d( n+1 ), e( n+1 ), b,
502 $ lda, x, lda, rcond, rwork, rwork( nrhs+1 ),
503 $ work, info )
504*
505* Check the error code from DPTSVX.
506*
507 IF( info.NE.izero )
508 $ CALL alaerh( path, 'DPTSVX', info, izero, fact, n, n,
509 $ 1, 1, nrhs, imat, nfail, nerrs, nout )
510 IF( izero.EQ.0 ) THEN
511 IF( ifact.EQ.2 ) THEN
512*
513* Check the factorization by computing the ratio
514* norm(L*D*L' - A) / (n * norm(A) * EPS )
515*
516 k1 = 1
517 CALL dptt01( n, d, e, d( n+1 ), e( n+1 ), work,
518 $ result( 1 ) )
519 ELSE
520 k1 = 2
521 END IF
522*
523* Compute the residual in the solution.
524*
525 CALL dlacpy( 'Full', n, nrhs, b, lda, work, lda )
526 CALL dptt02( n, nrhs, d, e, x, lda, work, lda,
527 $ result( 2 ) )
528*
529* Check solution from generated exact solution.
530*
531 CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
532 $ result( 3 ) )
533*
534* Check error bounds from iterative refinement.
535*
536 CALL dptt05( n, nrhs, d, e, b, lda, x, lda, xact, lda,
537 $ rwork, rwork( nrhs+1 ), result( 4 ) )
538 ELSE
539 k1 = 6
540 END IF
541*
542* Check the reciprocal of the condition number.
543*
544 result( 6 ) = dget06( rcond, rcondc )
545*
546* Print information about the tests that did not pass
547* the threshold.
548*
549 DO 90 k = k1, 6
550 IF( result( k ).GE.thresh ) THEN
551 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
552 $ CALL aladhd( nout, path )
553 WRITE( nout, fmt = 9998 )'DPTSVX', fact, n, imat,
554 $ k, result( k )
555 nfail = nfail + 1
556 END IF
557 90 CONTINUE
558 nrun = nrun + 7 - k1
559 100 CONTINUE
560 110 CONTINUE
561 120 CONTINUE
562*
563* Print a summary of the results.
564*
565 CALL alasvm( path, nout, nfail, nrun, nerrs )
566*
567 9999 FORMAT( 1x, a, ', N =', i5, ', type ', i2, ', test ', i2,
568 $ ', ratio = ', g12.5 )
569 9998 FORMAT( 1x, a, ', FACT=''', a1, ''', N =', i5, ', type ', i2,
570 $ ', test ', i2, ', ratio = ', g12.5 )
571 RETURN
572*
573* End of DDRVPT
574*
575 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine aladhd(iounit, path)
ALADHD
Definition aladhd.f:90
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
subroutine ddrvpt(dotype, nn, nval, nrhs, thresh, tsterr, a, d, e, b, x, xact, work, rwork, nout)
DDRVPT
Definition ddrvpt.f:140
subroutine derrvx(path, nunit)
DERRVX
Definition derrvx.f:55
subroutine dget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
DGET04
Definition dget04.f:102
subroutine dlaptm(n, nrhs, alpha, d, e, x, ldx, beta, b, ldb)
DLAPTM
Definition dlaptm.f:116
subroutine dlatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
DLATB4
Definition dlatb4.f:120
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
subroutine dptt01(n, d, e, df, ef, work, resid)
DPTT01
Definition dptt01.f:91
subroutine dptt02(n, nrhs, d, e, x, ldx, b, ldb, resid)
DPTT02
Definition dptt02.f:104
subroutine dptt05(n, nrhs, d, e, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
DPTT05
Definition dptt05.f:150
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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 dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dptsv(n, nrhs, d, e, b, ldb, info)
DPTSV computes the solution to system of linear equations A * X = B for PT matrices
Definition dptsv.f:114
subroutine dptsvx(fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx, rcond, ferr, berr, work, info)
DPTSVX computes the solution to system of linear equations A * X = B for PT matrices
Definition dptsvx.f:228
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
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79