LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, convertion 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_DLATMS
          D_WORK_DLATMS is DOUBLE PRECISION array, dimension ( 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.
Date
November 2013

Definition at line 245 of file ddrvrfp.f.

245 *
246 * -- LAPACK test routine (version 3.5.0) --
247 * -- LAPACK is a software package provided by Univ. of Tennessee, --
248 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249 * November 2013
250 *
251 * .. Scalar Arguments ..
252  INTEGER nn, nns, nnt, nout
253  DOUBLE PRECISION thresh
254 * ..
255 * .. Array Arguments ..
256  INTEGER nval( nn ), nsval( nns ), ntval( nnt )
257  DOUBLE PRECISION a( * )
258  DOUBLE PRECISION ainv( * )
259  DOUBLE PRECISION asav( * )
260  DOUBLE PRECISION b( * )
261  DOUBLE PRECISION bsav( * )
262  DOUBLE PRECISION afac( * )
263  DOUBLE PRECISION arf( * )
264  DOUBLE PRECISION arfinv( * )
265  DOUBLE PRECISION xact( * )
266  DOUBLE PRECISION x( * )
267  DOUBLE PRECISION d_work_dlatms( * )
268  DOUBLE PRECISION d_work_dpot01( * )
269  DOUBLE PRECISION d_temp_dpot02( * )
270  DOUBLE PRECISION d_temp_dpot03( * )
271  DOUBLE PRECISION d_work_dlansy( * )
272  DOUBLE PRECISION d_work_dpot02( * )
273  DOUBLE PRECISION d_work_dpot03( * )
274 * ..
275 *
276 * =====================================================================
277 *
278 * .. Parameters ..
279  DOUBLE PRECISION one, zero
280  parameter ( one = 1.0d+0, zero = 0.0d+0 )
281  INTEGER ntests
282  parameter ( ntests = 4 )
283 * ..
284 * .. Local Scalars ..
285  LOGICAL zerot
286  INTEGER i, info, iuplo, lda, ldb, imat, nerrs, nfail,
287  + nrhs, nrun, izero, ioff, k, nt, n, iform, iin,
288  + iit, iis
289  CHARACTER dist, ctype, uplo, cform
290  INTEGER kl, ku, mode
291  DOUBLE PRECISION anorm, ainvnm, cndnum, rcondc
292 * ..
293 * .. Local Arrays ..
294  CHARACTER uplos( 2 ), forms( 2 )
295  INTEGER iseed( 4 ), iseedy( 4 )
296  DOUBLE PRECISION result( ntests )
297 * ..
298 * .. External Functions ..
299  DOUBLE PRECISION dlansy
300  EXTERNAL dlansy
301 * ..
302 * .. External Subroutines ..
303  EXTERNAL aladhd, alaerh, alasvm, dget04, dtfttr, dlacpy,
304  + dlarhs, dlatb4, dlatms, dpftri, dpftrf, dpftrs,
305  + dpot01, dpot02, dpot03, dpotri, dpotrf, dtrttf
306 * ..
307 * .. Scalars in Common ..
308  CHARACTER*32 srnamt
309 * ..
310 * .. Common blocks ..
311  COMMON / srnamc / srnamt
312 * ..
313 * .. Data statements ..
314  DATA iseedy / 1988, 1989, 1990, 1991 /
315  DATA uplos / 'U', 'L' /
316  DATA forms / 'N', 'T' /
317 * ..
318 * .. Executable Statements ..
319 *
320 * Initialize constants and the random number seed.
321 *
322  nrun = 0
323  nfail = 0
324  nerrs = 0
325  DO 10 i = 1, 4
326  iseed( i ) = iseedy( i )
327  10 CONTINUE
328 *
329  DO 130 iin = 1, nn
330 *
331  n = nval( iin )
332  lda = max( n, 1 )
333  ldb = max( n, 1 )
334 *
335  DO 980 iis = 1, nns
336 *
337  nrhs = nsval( iis )
338 *
339  DO 120 iit = 1, nnt
340 *
341  imat = ntval( iit )
342 *
343 * If N.EQ.0, only consider the first type
344 *
345  IF( n.EQ.0 .AND. iit.GE.1 ) GO TO 120
346 *
347 * Skip types 3, 4, or 5 if the matrix size is too small.
348 *
349  IF( imat.EQ.4 .AND. n.LE.1 ) GO TO 120
350  IF( imat.EQ.5 .AND. n.LE.2 ) GO TO 120
351 *
352 * Do first for UPLO = 'U', then for UPLO = 'L'
353 *
354  DO 110 iuplo = 1, 2
355  uplo = uplos( iuplo )
356 *
357 * Do first for CFORM = 'N', then for CFORM = 'C'
358 *
359  DO 100 iform = 1, 2
360  cform = forms( iform )
361 *
362 * Set up parameters with DLATB4 and generate a test
363 * matrix with DLATMS.
364 *
365  CALL dlatb4( 'DPO', imat, n, n, ctype, kl, ku,
366  + anorm, mode, cndnum, dist )
367 *
368  srnamt = 'DLATMS'
369  CALL dlatms( n, n, dist, iseed, ctype,
370  + d_work_dlatms,
371  + mode, cndnum, anorm, kl, ku, uplo, a,
372  + lda, d_work_dlatms, info )
373 *
374 * Check error code from DLATMS.
375 *
376  IF( info.NE.0 ) THEN
377  CALL alaerh( 'DPF', 'DLATMS', info, 0, uplo, n,
378  + n, -1, -1, -1, iit, nfail, nerrs,
379  + nout )
380  GO TO 100
381  END IF
382 *
383 * For types 3-5, zero one row and column of the matrix to
384 * test that INFO is returned correctly.
385 *
386  zerot = imat.GE.3 .AND. imat.LE.5
387  IF( zerot ) THEN
388  IF( iit.EQ.3 ) THEN
389  izero = 1
390  ELSE IF( iit.EQ.4 ) THEN
391  izero = n
392  ELSE
393  izero = n / 2 + 1
394  END IF
395  ioff = ( izero-1 )*lda
396 *
397 * Set row and column IZERO of A to 0.
398 *
399  IF( iuplo.EQ.1 ) THEN
400  DO 20 i = 1, izero - 1
401  a( ioff+i ) = zero
402  20 CONTINUE
403  ioff = ioff + izero
404  DO 30 i = izero, n
405  a( ioff ) = zero
406  ioff = ioff + lda
407  30 CONTINUE
408  ELSE
409  ioff = izero
410  DO 40 i = 1, izero - 1
411  a( ioff ) = zero
412  ioff = ioff + lda
413  40 CONTINUE
414  ioff = ioff - izero
415  DO 50 i = izero, n
416  a( ioff+i ) = zero
417  50 CONTINUE
418  END IF
419  ELSE
420  izero = 0
421  END IF
422 *
423 * Save a copy of the matrix A in ASAV.
424 *
425  CALL dlacpy( uplo, n, n, a, lda, asav, lda )
426 *
427 * Compute the condition number of A (RCONDC).
428 *
429  IF( zerot ) THEN
430  rcondc = zero
431  ELSE
432 *
433 * Compute the 1-norm of A.
434 *
435  anorm = dlansy( '1', uplo, n, a, lda,
436  + d_work_dlansy )
437 *
438 * Factor the matrix A.
439 *
440  CALL dpotrf( uplo, n, a, lda, info )
441 *
442 * Form the inverse of A.
443 *
444  CALL dpotri( uplo, n, a, lda, info )
445 
446  IF ( n .NE. 0 ) THEN
447 
448 *
449 * Compute the 1-norm condition number of A.
450 *
451  ainvnm = dlansy( '1', uplo, n, a, lda,
452  + d_work_dlansy )
453  rcondc = ( one / anorm ) / ainvnm
454 *
455 * Restore the matrix A.
456 *
457  CALL dlacpy( uplo, n, n, asav, lda, a, lda )
458  END IF
459 *
460  END IF
461 *
462 * Form an exact solution and set the right hand side.
463 *
464  srnamt = 'DLARHS'
465  CALL dlarhs( 'DPO', 'N', uplo, ' ', n, n, kl, ku,
466  + nrhs, a, lda, xact, lda, b, lda,
467  + iseed, info )
468  CALL dlacpy( 'Full', n, nrhs, b, lda, bsav, lda )
469 *
470 * Compute the L*L' or U'*U factorization of the
471 * matrix and solve the system.
472 *
473  CALL dlacpy( uplo, n, n, a, lda, afac, lda )
474  CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldb )
475 *
476  srnamt = 'DTRTTF'
477  CALL dtrttf( cform, uplo, n, afac, lda, arf, info )
478  srnamt = 'DPFTRF'
479  CALL dpftrf( cform, uplo, n, arf, info )
480 *
481 * Check error code from DPFTRF.
482 *
483  IF( info.NE.izero ) THEN
484 *
485 * LANGOU: there is a small hick here: IZERO should
486 * always be INFO however if INFO is ZERO, ALAERH does not
487 * complain.
488 *
489  CALL alaerh( 'DPF', 'DPFSV ', info, izero,
490  + uplo, n, n, -1, -1, nrhs, iit,
491  + nfail, nerrs, nout )
492  GO TO 100
493  END IF
494 *
495 * Skip the tests if INFO is not 0.
496 *
497  IF( info.NE.0 ) THEN
498  GO TO 100
499  END IF
500 *
501  srnamt = 'DPFTRS'
502  CALL dpftrs( cform, uplo, n, nrhs, arf, x, ldb,
503  + info )
504 *
505  srnamt = 'DTFTTR'
506  CALL dtfttr( cform, uplo, n, arf, afac, lda, info )
507 *
508 * Reconstruct matrix from factors and compute
509 * residual.
510 *
511  CALL dlacpy( uplo, n, n, afac, lda, asav, lda )
512  CALL dpot01( uplo, n, a, lda, afac, lda,
513  + d_work_dpot01, result( 1 ) )
514  CALL dlacpy( uplo, n, n, asav, lda, afac, lda )
515 *
516 * Form the inverse and compute the residual.
517 *
518  IF(mod(n,2).EQ.0)THEN
519  CALL dlacpy( 'A', n+1, n/2, arf, n+1, arfinv,
520  + n+1 )
521  ELSE
522  CALL dlacpy( 'A', n, (n+1)/2, arf, n, arfinv,
523  + n )
524  END IF
525 *
526  srnamt = 'DPFTRI'
527  CALL dpftri( cform, uplo, n, arfinv , info )
528 *
529  srnamt = 'DTFTTR'
530  CALL dtfttr( cform, uplo, n, arfinv, ainv, lda,
531  + info )
532 *
533 * Check error code from DPFTRI.
534 *
535  IF( info.NE.0 )
536  + CALL alaerh( 'DPO', 'DPFTRI', info, 0, uplo, n,
537  + n, -1, -1, -1, imat, nfail, nerrs,
538  + nout )
539 *
540  CALL dpot03( uplo, n, a, lda, ainv, lda,
541  + d_temp_dpot03, lda, d_work_dpot03,
542  + rcondc, result( 2 ) )
543 *
544 * Compute residual of the computed solution.
545 *
546  CALL dlacpy( 'Full', n, nrhs, b, lda,
547  + d_temp_dpot02, lda )
548  CALL dpot02( uplo, n, nrhs, a, lda, x, lda,
549  + d_temp_dpot02, lda, d_work_dpot02,
550  + result( 3 ) )
551 *
552 * Check solution from generated exact solution.
553 
554  CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
555  + result( 4 ) )
556  nt = 4
557 *
558 * Print information about the tests that did not
559 * pass the threshold.
560 *
561  DO 60 k = 1, nt
562  IF( result( k ).GE.thresh ) THEN
563  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
564  + CALL aladhd( nout, 'DPF' )
565  WRITE( nout, fmt = 9999 )'DPFSV ', uplo,
566  + n, iit, k, result( k )
567  nfail = nfail + 1
568  END IF
569  60 CONTINUE
570  nrun = nrun + nt
571  100 CONTINUE
572  110 CONTINUE
573  120 CONTINUE
574  980 CONTINUE
575  130 CONTINUE
576 *
577 * Print a summary of the results.
578 *
579  CALL alasvm( 'DPF', nout, nfail, nrun, nerrs )
580 *
581  9999 FORMAT( 1x, a6, ', UPLO=''', a1, ''', N =', i5, ', type ', i1,
582  + ', test(', i1, ')=', g12.5 )
583 *
584  RETURN
585 *
586 * End of DDRVRFP
587 *
alasvm
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
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, or the element of largest absolute value of a real symmetric matrix.
Definition: dlansy.f:124
alaerh
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:149
dlarhs
subroutine dlarhs(PATH, XTYPE, UPLO, TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, LDB, ISEED, INFO)
DLARHS
Definition: dlarhs.f:206
dpotrf
subroutine dpotrf(UPLO, N, A, LDA, INFO)
DPOTRF
Definition: dpotrf.f:109
dpftri
subroutine dpftri(TRANSR, UPLO, N, A, INFO)
DPFTRI
Definition: dpftri.f:193
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:198
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:105
dpot03
subroutine dpot03(UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK, RCOND, RESID)
DPOT03
Definition: dpot03.f:127
dlatb4
subroutine dlatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
DLATB4
Definition: dlatb4.f:122
dget04
subroutine dget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
DGET04
Definition: dget04.f:104
aladhd
subroutine aladhd(IOUNIT, PATH)
ALADHD
Definition: aladhd.f:80
dpftrf
subroutine dpftrf(TRANSR, UPLO, N, A, INFO)
DPFTRF
Definition: dpftrf.f:200
dpot01
subroutine dpot01(UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID)
DPOT01
Definition: dpot01.f:106
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:196
dlatms
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:323
dpftrs
subroutine dpftrs(TRANSR, UPLO, N, NRHS, A, B, LDB, INFO)
DPFTRS
Definition: dpftrs.f:201
dpot02
subroutine dpot02(UPLO, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
DPOT02
Definition: dpot02.f:129
dpotri
subroutine dpotri(UPLO, N, A, LDA, INFO)
DPOTRI
Definition: dpotri.f:97

Here is the call graph for this function:

Here is the caller graph for this function: