LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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,
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 *
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 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 dlarhs(PATH, XTYPE, UPLO, TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, LDB, ISEED, INFO)
DLARHS
Definition: dlarhs.f:205
subroutine dget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
DGET04
Definition: dget04.f:102
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 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 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
subroutine dpftrs(TRANSR, UPLO, N, NRHS, A, B, LDB, INFO)
DPFTRS
Definition: dpftrs.f:199
subroutine dpftri(TRANSR, UPLO, N, A, INFO)
DPFTRI
Definition: dpftri.f:191
subroutine dpftrf(TRANSR, UPLO, N, A, INFO)
DPFTRF
Definition: dpftrf.f:198
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 dpotrf(UPLO, N, A, LDA, INFO)
DPOTRF
Definition: dpotrf.f:107
subroutine dpotri(UPLO, N, A, LDA, INFO)
DPOTRI
Definition: dpotri.f:95
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
Here is the call graph for this function:
Here is the caller graph for this function: