LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ddrvrfp.f
Go to the documentation of this file.
1*> \brief \b DDRVRFP
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 DDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL,
12* + THRESH, A, ASAV, AFAC, AINV, B,
13* + BSAV, XACT, X, ARF, ARFINV,
14* + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02,
15* + D_TEMP_DPOT03, D_WORK_DLANSY,
16* + D_WORK_DPOT02, D_WORK_DPOT03 )
17*
18* .. Scalar Arguments ..
19* INTEGER NN, NNS, NNT, NOUT
20* DOUBLE PRECISION THRESH
21* ..
22* .. Array Arguments ..
23* INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
24* DOUBLE PRECISION A( * )
25* DOUBLE PRECISION AINV( * )
26* DOUBLE PRECISION ASAV( * )
27* DOUBLE PRECISION B( * )
28* DOUBLE PRECISION BSAV( * )
29* DOUBLE PRECISION AFAC( * )
30* DOUBLE PRECISION ARF( * )
31* DOUBLE PRECISION ARFINV( * )
32* DOUBLE PRECISION XACT( * )
33* DOUBLE PRECISION X( * )
34* DOUBLE PRECISION D_WORK_DLATMS( * )
35* DOUBLE PRECISION D_WORK_DPOT01( * )
36* DOUBLE PRECISION D_TEMP_DPOT02( * )
37* DOUBLE PRECISION D_TEMP_DPOT03( * )
38* DOUBLE PRECISION D_WORK_DLANSY( * )
39* DOUBLE PRECISION D_WORK_DPOT02( * )
40* DOUBLE PRECISION D_WORK_DPOT03( * )
41* ..
42*
43*
44*> \par Purpose:
45* =============
46*>
47*> \verbatim
48*>
49*> DDRVRFP tests the LAPACK RFP routines:
50*> DPFTRF, DPFTRS, and DPFTRI.
51*>
52*> This testing routine follow the same tests as DDRVPO (test for the full
53*> format Symmetric Positive Definite solver).
54*>
55*> The tests are performed in Full Format, conversion back and forth from
56*> full format to RFP format are performed using the routines DTRTTF and
57*> DTFTTR.
58*>
59*> First, a specific matrix A of size N is created. There is nine types of
60*> different matrixes possible.
61*> 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS)
62*> 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS
63*> *3. First row and column zero 8. Scaled near underflow
64*> *4. Last row and column zero 9. Scaled near overflow
65*> *5. Middle row and column zero
66*> (* - tests error exits from DPFTRF, no test ratios are computed)
67*> A solution XACT of size N-by-NRHS is created and the associated right
68*> hand side B as well. Then DPFTRF is called to compute L (or U), the
69*> Cholesky factor of A. Then L (or U) is used to solve the linear system
70*> of equations AX = B. This gives X. Then L (or U) is used to compute the
71*> inverse of A, AINV. The following four tests are then performed:
72*> (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or
73*> norm( U'*U - A ) / ( N * norm(A) * EPS ),
74*> (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
75*> (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
76*> (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
77*> where EPS is the machine precision, RCOND the condition number of A, and
78*> norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4).
79*> Errors occur when INFO parameter is not as expected. Failures occur when
80*> a test ratios is greater than THRES.
81*> \endverbatim
82*
83* Arguments:
84* ==========
85*
86*> \param[in] NOUT
87*> \verbatim
88*> NOUT is INTEGER
89*> The unit number for output.
90*> \endverbatim
91*>
92*> \param[in] NN
93*> \verbatim
94*> NN is INTEGER
95*> The number of values of N contained in the vector NVAL.
96*> \endverbatim
97*>
98*> \param[in] NVAL
99*> \verbatim
100*> NVAL is INTEGER array, dimension (NN)
101*> The values of the matrix dimension N.
102*> \endverbatim
103*>
104*> \param[in] NNS
105*> \verbatim
106*> NNS is INTEGER
107*> The number of values of NRHS contained in the vector NSVAL.
108*> \endverbatim
109*>
110*> \param[in] NSVAL
111*> \verbatim
112*> NSVAL is INTEGER array, dimension (NNS)
113*> The values of the number of right-hand sides NRHS.
114*> \endverbatim
115*>
116*> \param[in] NNT
117*> \verbatim
118*> NNT is INTEGER
119*> The number of values of MATRIX TYPE contained in the vector NTVAL.
120*> \endverbatim
121*>
122*> \param[in] NTVAL
123*> \verbatim
124*> NTVAL is INTEGER array, dimension (NNT)
125*> The values of matrix type (between 0 and 9 for PO/PP/PF matrices).
126*> \endverbatim
127*>
128*> \param[in] THRESH
129*> \verbatim
130*> THRESH is DOUBLE PRECISION
131*> The threshold value for the test ratios. A result is
132*> included in the output file if RESULT >= THRESH. To have
133*> every test ratio printed, use THRESH = 0.
134*> \endverbatim
135*>
136*> \param[out] A
137*> \verbatim
138*> A is DOUBLE PRECISION array, dimension (NMAX*NMAX)
139*> \endverbatim
140*>
141*> \param[out] ASAV
142*> \verbatim
143*> ASAV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
144*> \endverbatim
145*>
146*> \param[out] AFAC
147*> \verbatim
148*> AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX)
149*> \endverbatim
150*>
151*> \param[out] AINV
152*> \verbatim
153*> AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
154*> \endverbatim
155*>
156*> \param[out] B
157*> \verbatim
158*> B is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
159*> \endverbatim
160*>
161*> \param[out] BSAV
162*> \verbatim
163*> BSAV is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
164*> \endverbatim
165*>
166*> \param[out] XACT
167*> \verbatim
168*> XACT is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
169*> \endverbatim
170*>
171*> \param[out] X
172*> \verbatim
173*> X is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
174*> \endverbatim
175*>
176*> \param[out] ARF
177*> \verbatim
178*> ARF is DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
179*> \endverbatim
180*>
181*> \param[out] ARFINV
182*> \verbatim
183*> ARFINV is DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
184*> \endverbatim
185*>
186*> \param[out] D_WORK_DLATMS
187*> \verbatim
188*> D_WORK_DLATMS is DOUBLE PRECISION array, dimension ( 3*NMAX )
189*> \endverbatim
190*>
191*> \param[out] D_WORK_DPOT01
192*> \verbatim
193*> D_WORK_DPOT01 is DOUBLE PRECISION array, dimension ( NMAX )
194*> \endverbatim
195*>
196*> \param[out] D_TEMP_DPOT02
197*> \verbatim
198*> D_TEMP_DPOT02 is DOUBLE PRECISION array, dimension ( NMAX*MAXRHS )
199*> \endverbatim
200*>
201*> \param[out] D_TEMP_DPOT03
202*> \verbatim
203*> D_TEMP_DPOT03 is DOUBLE PRECISION array, dimension ( NMAX*NMAX )
204*> \endverbatim
205*>
206*> \param[out] D_WORK_DLANSY
207*> \verbatim
208*> D_WORK_DLANSY is DOUBLE PRECISION array, dimension ( NMAX )
209*> \endverbatim
210*>
211*> \param[out] D_WORK_DPOT02
212*> \verbatim
213*> D_WORK_DPOT02 is DOUBLE PRECISION array, dimension ( NMAX )
214*> \endverbatim
215*>
216*> \param[out] D_WORK_DPOT03
217*> \verbatim
218*> D_WORK_DPOT03 is DOUBLE PRECISION array, dimension ( NMAX )
219*> \endverbatim
220*
221* Authors:
222* ========
223*
224*> \author Univ. of Tennessee
225*> \author Univ. of California Berkeley
226*> \author Univ. of Colorado Denver
227*> \author NAG Ltd.
228*
229*> \ingroup double_lin
230*
231* =====================================================================
232 SUBROUTINE ddrvrfp( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL,
233 + THRESH, A, ASAV, AFAC, AINV, B,
234 + BSAV, XACT, X, ARF, ARFINV,
235 + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02,
236 + D_TEMP_DPOT03, D_WORK_DLANSY,
237 + D_WORK_DPOT02, D_WORK_DPOT03 )
238*
239* -- LAPACK test routine --
240* -- LAPACK is a software package provided by Univ. of Tennessee, --
241* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242*
243* .. Scalar Arguments ..
244 INTEGER NN, NNS, NNT, NOUT
245 DOUBLE PRECISION THRESH
246* ..
247* .. Array Arguments ..
248 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
249 DOUBLE PRECISION A( * )
250 DOUBLE PRECISION AINV( * )
251 DOUBLE PRECISION ASAV( * )
252 DOUBLE PRECISION B( * )
253 DOUBLE PRECISION BSAV( * )
254 DOUBLE PRECISION AFAC( * )
255 DOUBLE PRECISION ARF( * )
256 DOUBLE PRECISION ARFINV( * )
257 DOUBLE PRECISION XACT( * )
258 DOUBLE PRECISION X( * )
259 DOUBLE PRECISION D_WORK_DLATMS( * )
260 DOUBLE PRECISION D_WORK_DPOT01( * )
261 DOUBLE PRECISION D_TEMP_DPOT02( * )
262 DOUBLE PRECISION D_TEMP_DPOT03( * )
263 DOUBLE PRECISION D_WORK_DLANSY( * )
264 DOUBLE PRECISION D_WORK_DPOT02( * )
265 DOUBLE PRECISION D_WORK_DPOT03( * )
266* ..
267*
268* =====================================================================
269*
270* .. Parameters ..
271 DOUBLE PRECISION ONE, ZERO
272 PARAMETER ( ONE = 1.0d+0, zero = 0.0d+0 )
273 INTEGER NTESTS
274 PARAMETER ( NTESTS = 4 )
275* ..
276* .. Local Scalars ..
277 LOGICAL ZEROT
278 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL,
279 + nrhs, nrun, izero, ioff, k, nt, n, iform, iin,
280 + iit, iis
281 CHARACTER DIST, CTYPE, UPLO, CFORM
282 INTEGER KL, KU, MODE
283 DOUBLE PRECISION ANORM, AINVNM, CNDNUM, RCONDC
284* ..
285* .. Local Arrays ..
286 CHARACTER UPLOS( 2 ), FORMS( 2 )
287 INTEGER ISEED( 4 ), ISEEDY( 4 )
288 DOUBLE PRECISION RESULT( NTESTS )
289* ..
290* .. External Functions ..
291 DOUBLE PRECISION DLANSY
292 EXTERNAL dlansy
293* ..
294* .. External Subroutines ..
295 EXTERNAL aladhd, alaerh, alasvm, dget04, dtfttr, dlacpy,
298* ..
299* .. Scalars in Common ..
300 CHARACTER*32 SRNAMT
301* ..
302* .. Common blocks ..
303 COMMON / SRNAMC / SRNAMT
304* ..
305* .. Data statements ..
306 DATA iseedy / 1988, 1989, 1990, 1991 /
307 DATA uplos / 'U', 'L' /
308 DATA forms / 'N', 'T' /
309* ..
310* .. Executable Statements ..
311*
312* Initialize constants and the random number seed.
313*
314 nrun = 0
315 nfail = 0
316 nerrs = 0
317 DO 10 i = 1, 4
318 iseed( i ) = iseedy( i )
319 10 CONTINUE
320*
321 DO 130 iin = 1, nn
322*
323 n = nval( iin )
324 lda = max( n, 1 )
325 ldb = max( n, 1 )
326*
327 DO 980 iis = 1, nns
328*
329 nrhs = nsval( iis )
330*
331 DO 120 iit = 1, nnt
332*
333 imat = ntval( iit )
334*
335* If N.EQ.0, only consider the first type
336*
337 IF( n.EQ.0 .AND. iit.GE.1 ) GO TO 120
338*
339* Skip types 3, 4, or 5 if the matrix size is too small.
340*
341 IF( imat.EQ.4 .AND. n.LE.1 ) GO TO 120
342 IF( imat.EQ.5 .AND. n.LE.2 ) GO TO 120
343*
344* Do first for UPLO = 'U', then for UPLO = 'L'
345*
346 DO 110 iuplo = 1, 2
347 uplo = uplos( iuplo )
348*
349* Do first for CFORM = 'N', then for CFORM = 'C'
350*
351 DO 100 iform = 1, 2
352 cform = forms( iform )
353*
354* Set up parameters with DLATB4 and generate a test
355* matrix with DLATMS.
356*
357 CALL dlatb4( 'DPO', imat, n, n, ctype, kl, ku,
358 + anorm, mode, cndnum, dist )
359*
360 srnamt = 'DLATMS'
361 CALL dlatms( n, n, dist, iseed, ctype,
362 + d_work_dlatms,
363 + mode, cndnum, anorm, kl, ku, uplo, a,
364 + lda, d_work_dlatms, info )
365*
366* Check error code from DLATMS.
367*
368 IF( info.NE.0 ) THEN
369 CALL alaerh( 'DPF', 'DLATMS', info, 0, uplo, n,
370 + n, -1, -1, -1, iit, nfail, nerrs,
371 + nout )
372 GO TO 100
373 END IF
374*
375* For types 3-5, zero one row and column of the matrix to
376* test that INFO is returned correctly.
377*
378 zerot = imat.GE.3 .AND. imat.LE.5
379 IF( zerot ) THEN
380 IF( iit.EQ.3 ) THEN
381 izero = 1
382 ELSE IF( iit.EQ.4 ) THEN
383 izero = n
384 ELSE
385 izero = n / 2 + 1
386 END IF
387 ioff = ( izero-1 )*lda
388*
389* Set row and column IZERO of A to 0.
390*
391 IF( iuplo.EQ.1 ) THEN
392 DO 20 i = 1, izero - 1
393 a( ioff+i ) = zero
394 20 CONTINUE
395 ioff = ioff + izero
396 DO 30 i = izero, n
397 a( ioff ) = zero
398 ioff = ioff + lda
399 30 CONTINUE
400 ELSE
401 ioff = izero
402 DO 40 i = 1, izero - 1
403 a( ioff ) = zero
404 ioff = ioff + lda
405 40 CONTINUE
406 ioff = ioff - izero
407 DO 50 i = izero, n
408 a( ioff+i ) = zero
409 50 CONTINUE
410 END IF
411 ELSE
412 izero = 0
413 END IF
414*
415* Save a copy of the matrix A in ASAV.
416*
417 CALL dlacpy( uplo, n, n, a, lda, asav, lda )
418*
419* Compute the condition number of A (RCONDC).
420*
421 IF( zerot ) THEN
422 rcondc = zero
423 ELSE
424*
425* Compute the 1-norm of A.
426*
427 anorm = dlansy( '1', uplo, n, a, lda,
428 + d_work_dlansy )
429*
430* Factor the matrix A.
431*
432 CALL dpotrf( uplo, n, a, lda, info )
433*
434* Form the inverse of A.
435*
436 CALL dpotri( uplo, n, a, lda, info )
437
438 IF ( n .NE. 0 ) THEN
439*
440* Compute the 1-norm condition number of A.
441*
442 ainvnm = dlansy( '1', uplo, n, a, lda,
443 + d_work_dlansy )
444 rcondc = ( one / anorm ) / ainvnm
445*
446* Restore the matrix A.
447*
448 CALL dlacpy( uplo, n, n, asav, lda, a, lda )
449 END IF
450*
451 END IF
452*
453* Form an exact solution and set the right hand side.
454*
455 srnamt = 'DLARHS'
456 CALL dlarhs( 'DPO', 'N', uplo, ' ', n, n, kl, ku,
457 + nrhs, a, lda, xact, lda, b, lda,
458 + iseed, info )
459 CALL dlacpy( 'Full', n, nrhs, b, lda, bsav, lda )
460*
461* Compute the L*L' or U'*U factorization of the
462* matrix and solve the system.
463*
464 CALL dlacpy( uplo, n, n, a, lda, afac, lda )
465 CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldb )
466*
467 srnamt = 'DTRTTF'
468 CALL dtrttf( cform, uplo, n, afac, lda, arf, info )
469 srnamt = 'DPFTRF'
470 CALL dpftrf( cform, uplo, n, arf, info )
471*
472* Check error code from DPFTRF.
473*
474 IF( info.NE.izero ) THEN
475*
476* LANGOU: there is a small hick here: IZERO should
477* always be INFO however if INFO is ZERO, ALAERH does not
478* complain.
479*
480 CALL alaerh( 'DPF', 'DPFSV ', info, izero,
481 + uplo, n, n, -1, -1, nrhs, iit,
482 + nfail, nerrs, nout )
483 GO TO 100
484 END IF
485*
486* Skip the tests if INFO is not 0.
487*
488 IF( info.NE.0 ) THEN
489 GO TO 100
490 END IF
491*
492 srnamt = 'DPFTRS'
493 CALL dpftrs( cform, uplo, n, nrhs, arf, x, ldb,
494 + info )
495*
496 srnamt = 'DTFTTR'
497 CALL dtfttr( cform, uplo, n, arf, afac, lda, info )
498*
499* Reconstruct matrix from factors and compute
500* residual.
501*
502 CALL dlacpy( uplo, n, n, afac, lda, asav, lda )
503 CALL dpot01( uplo, n, a, lda, afac, lda,
504 + d_work_dpot01, result( 1 ) )
505 CALL dlacpy( uplo, n, n, asav, lda, afac, lda )
506*
507* Form the inverse and compute the residual.
508*
509 IF(mod(n,2).EQ.0)THEN
510 CALL dlacpy( 'A', n+1, n/2, arf, n+1, arfinv,
511 + n+1 )
512 ELSE
513 CALL dlacpy( 'A', n, (n+1)/2, arf, n, arfinv,
514 + n )
515 END IF
516*
517 srnamt = 'DPFTRI'
518 CALL dpftri( cform, uplo, n, arfinv , info )
519*
520 srnamt = 'DTFTTR'
521 CALL dtfttr( cform, uplo, n, arfinv, ainv, lda,
522 + info )
523*
524* Check error code from DPFTRI.
525*
526 IF( info.NE.0 )
527 + CALL alaerh( 'DPO', 'DPFTRI', info, 0, uplo, n,
528 + n, -1, -1, -1, imat, nfail, nerrs,
529 + nout )
530*
531 CALL dpot03( uplo, n, a, lda, ainv, lda,
532 + d_temp_dpot03, lda, d_work_dpot03,
533 + rcondc, result( 2 ) )
534*
535* Compute residual of the computed solution.
536*
537 CALL dlacpy( 'Full', n, nrhs, b, lda,
538 + d_temp_dpot02, lda )
539 CALL dpot02( uplo, n, nrhs, a, lda, x, lda,
540 + d_temp_dpot02, lda, d_work_dpot02,
541 + result( 3 ) )
542*
543* Check solution from generated exact solution.
544
545 CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
546 + result( 4 ) )
547 nt = 4
548*
549* Print information about the tests that did not
550* pass the threshold.
551*
552 DO 60 k = 1, nt
553 IF( result( k ).GE.thresh ) THEN
554 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
555 + CALL aladhd( nout, 'DPF' )
556 WRITE( nout, fmt = 9999 )'DPFSV ', uplo,
557 + n, iit, k, result( k )
558 nfail = nfail + 1
559 END IF
560 60 CONTINUE
561 nrun = nrun + nt
562 100 CONTINUE
563 110 CONTINUE
564 120 CONTINUE
565 980 CONTINUE
566 130 CONTINUE
567*
568* Print a summary of the results.
569*
570 CALL alasvm( 'DPF', nout, nfail, nrun, nerrs )
571*
572 9999 FORMAT( 1x, a6, ', UPLO=''', a1, ''', N =', i5, ', type ', i1,
573 + ', test(', i1, ')=', g12.5 )
574*
575 RETURN
576*
577* End of DDRVRFP
578*
579 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine dlarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
DLARHS
Definition dlarhs.f:205
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 ddrvrfp(nout, nn, nval, nns, nsval, nnt, ntval, thresh, a, asav, afac, ainv, b, bsav, xact, x, arf, arfinv, d_work_dlatms, d_work_dpot01, d_temp_dpot02, d_temp_dpot03, d_work_dlansy, d_work_dpot02, d_work_dpot03)
DDRVRFP
Definition ddrvrfp.f:238
subroutine dget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
DGET04
Definition dget04.f:102
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 dpot01(uplo, n, a, lda, afac, ldafac, rwork, resid)
DPOT01
Definition dpot01.f:104
subroutine dpot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DPOT02
Definition dpot02.f:127
subroutine dpot03(uplo, n, a, lda, ainv, ldainv, work, ldwork, rwork, rcond, resid)
DPOT03
Definition dpot03.f:125
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 dpftrf(transr, uplo, n, a, info)
DPFTRF
Definition dpftrf.f:198
subroutine dpftri(transr, uplo, n, a, info)
DPFTRI
Definition dpftri.f:191
subroutine dpftrs(transr, uplo, n, nrhs, a, b, ldb, info)
DPFTRS
Definition dpftrs.f:199
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
subroutine dpotri(uplo, n, a, lda, info)
DPOTRI
Definition dpotri.f:95
subroutine dtfttr(transr, uplo, n, arf, a, lda, info)
DTFTTR copies a triangular matrix from the rectangular full packed format (TF) to the standard full f...
Definition dtfttr.f:196
subroutine dtrttf(transr, uplo, n, a, lda, arf, info)
DTRTTF copies a triangular matrix from the standard full format (TR) to the rectangular full packed f...
Definition dtrttf.f:194