LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ddrvrfp()

subroutine ddrvrfp ( integer nout,
integer nn,
integer, dimension( nn ) nval,
integer nns,
integer, dimension( nns ) nsval,
integer nnt,
integer, dimension( nnt ) ntval,
double precision thresh,
double precision, dimension( * ) a,
double precision, dimension( * ) asav,
double precision, dimension( * ) afac,
double precision, dimension( * ) ainv,
double precision, dimension( * ) b,
double precision, dimension( * ) bsav,
double precision, dimension( * ) xact,
double precision, dimension( * ) x,
double precision, dimension( * ) arf,
double precision, dimension( * ) arfinv,
double precision, dimension( * ) d_work_dlatms,
double precision, dimension( * ) d_work_dpot01,
double precision, dimension( * ) d_temp_dpot02,
double precision, dimension( * ) d_temp_dpot03,
double precision, dimension( * ) d_work_dlansy,
double precision, dimension( * ) d_work_dpot02,
double precision, dimension( * ) d_work_dpot03 )

DDRVRFP

Purpose:
!>
!> DDRVRFP tests the LAPACK RFP routines:
!>     DPFTRF, DPFTRS, and DPFTRI.
!>
!> This testing routine follow the same tests as DDRVPO (test for the full
!> format Symmetric Positive Definite solver).
!>
!> The tests are performed in Full Format, conversion back and forth from
!> full format to RFP format are performed using the routines DTRTTF and
!> DTFTTR.
!>
!> First, a specific matrix A of size N is created. There is nine types of
!> different matrixes possible.
!>  1. Diagonal                        6. Random, CNDNUM = sqrt(0.1/EPS)
!>  2. Random, CNDNUM = 2              7. Random, CNDNUM = 0.1/EPS
!> *3. First row and column zero       8. Scaled near underflow
!> *4. Last row and column zero        9. Scaled near overflow
!> *5. Middle row and column zero
!> (* - tests error exits from DPFTRF, no test ratios are computed)
!> A solution XACT of size N-by-NRHS is created and the associated right
!> hand side B as well. Then DPFTRF is called to compute L (or U), the
!> Cholesky factor of A. Then L (or U) is used to solve the linear system
!> of equations AX = B. This gives X. Then L (or U) is used to compute the
!> inverse of A, AINV. The following four tests are then performed:
!> (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or
!>     norm( U'*U - A ) / ( N * norm(A) * EPS ),
!> (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
!> (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
!> (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
!> where EPS is the machine precision, RCOND the condition number of A, and
!> norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4).
!> Errors occur when INFO parameter is not as expected. Failures occur when
!> a test ratios is greater than THRES.
!>
Parameters
[in]NOUT
!>          NOUT is INTEGER
!>                The unit number for output.
!>
[in]NN
!>          NN is INTEGER
!>                The number of values of N contained in the vector NVAL.
!>
[in]NVAL
!>          NVAL is INTEGER array, dimension (NN)
!>                The values of the matrix dimension N.
!>
[in]NNS
!>          NNS is INTEGER
!>                The number of values of NRHS contained in the vector NSVAL.
!>
[in]NSVAL
!>          NSVAL is INTEGER array, dimension (NNS)
!>                The values of the number of right-hand sides NRHS.
!>
[in]NNT
!>          NNT is INTEGER
!>                The number of values of MATRIX TYPE contained in the vector NTVAL.
!>
[in]NTVAL
!>          NTVAL is INTEGER array, dimension (NNT)
!>                The values of matrix type (between 0 and 9 for PO/PP/PF matrices).
!>
[in]THRESH
!>          THRESH is DOUBLE PRECISION
!>                The threshold value for the test ratios.  A result is
!>                included in the output file if RESULT >= THRESH.  To have
!>                every test ratio printed, use THRESH = 0.
!>
[out]A
!>          A is DOUBLE PRECISION array, dimension (NMAX*NMAX)
!>
[out]ASAV
!>          ASAV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
!>
[out]AFAC
!>          AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX)
!>
[out]AINV
!>          AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
!>
[out]B
!>          B is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
!>
[out]BSAV
!>          BSAV is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
!>
[out]XACT
!>          XACT is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
!>
[out]X
!>          X is DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
!>
[out]ARF
!>          ARF is DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
!>
[out]ARFINV
!>          ARFINV is DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
!>
[out]D_WORK_DLATMS
!>          D_WORK_DLATMS is DOUBLE PRECISION array, dimension ( 3*NMAX )
!>
[out]D_WORK_DPOT01
!>          D_WORK_DPOT01 is DOUBLE PRECISION array, dimension ( NMAX )
!>
[out]D_TEMP_DPOT02
!>          D_TEMP_DPOT02 is DOUBLE PRECISION array, dimension ( NMAX*MAXRHS )
!>
[out]D_TEMP_DPOT03
!>          D_TEMP_DPOT03 is DOUBLE PRECISION array, dimension ( NMAX*NMAX )
!>
[out]D_WORK_DLANSY
!>          D_WORK_DLANSY is DOUBLE PRECISION array, dimension ( NMAX )
!>
[out]D_WORK_DPOT02
!>          D_WORK_DPOT02 is DOUBLE PRECISION array, dimension ( NMAX )
!>
[out]D_WORK_DPOT03
!>          D_WORK_DPOT03 is DOUBLE PRECISION array, dimension ( NMAX )
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 232 of file ddrvrfp.f.

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,
296 + dlarhs, dlatb4, dlatms, dpftri, dpftrf, dpftrs,
297 + dpot01, dpot02, dpot03, dpotri, dpotrf, dtrttf
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*
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
dlarhs
subroutine dlarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
DLARHS
Definition dlarhs.f:205
aladhd
subroutine aladhd(iounit, path)
ALADHD
Definition aladhd.f:90
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
dget04
subroutine dget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
DGET04
Definition dget04.f:102
dlatb4
subroutine dlatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
DLATB4
Definition dlatb4.f:120
dlatms
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
dpot01
subroutine dpot01(uplo, n, a, lda, afac, ldafac, rwork, resid)
DPOT01
Definition dpot01.f:104
dpot02
subroutine dpot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DPOT02
Definition dpot02.f:127
dpot03
subroutine dpot03(uplo, n, a, lda, ainv, ldainv, work, ldwork, rwork, rcond, resid)
DPOT03
Definition dpot03.f:125
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlansy
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansy.f:120
dpftrf
subroutine dpftrf(transr, uplo, n, a, info)
DPFTRF
Definition dpftrf.f:196
dpftri
subroutine dpftri(transr, uplo, n, a, info)
DPFTRI
Definition dpftri.f:189
dpftrs
subroutine dpftrs(transr, uplo, n, nrhs, a, b, ldb, info)
DPFTRS
Definition dpftrs.f:197
dpotrf
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:105
dpotri
subroutine dpotri(uplo, n, a, lda, info)
DPOTRI
Definition dpotri.f:93
dtfttr
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:194
dtrttf
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:192
Here is the call graph for this function:
Here is the caller graph for this function: