LAPACK 3.12.0
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:103
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:122
dpftrf
subroutine dpftrf(transr, uplo, n, a, info)
DPFTRF
Definition dpftrf.f:198
dpftri
subroutine dpftri(transr, uplo, n, a, info)
DPFTRI
Definition dpftri.f:191
dpftrs
subroutine dpftrs(transr, uplo, n, nrhs, a, b, ldb, info)
DPFTRS
Definition dpftrs.f:199
dpotrf
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
dpotri
subroutine dpotri(uplo, n, a, lda, info)
DPOTRI
Definition dpotri.f:95
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:196
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:194
Here is the call graph for this function:
Here is the caller graph for this function: