LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zdrvrfp ( integer  NOUT,
integer  NN,
integer, dimension( nn )  NVAL,
integer  NNS,
integer, dimension( nns )  NSVAL,
integer  NNT,
integer, dimension( nnt )  NTVAL,
double precision  THRESH,
complex*16, dimension( * )  A,
complex*16, dimension( * )  ASAV,
complex*16, dimension( * )  AFAC,
complex*16, dimension( * )  AINV,
complex*16, dimension( * )  B,
complex*16, dimension( * )  BSAV,
complex*16, dimension( * )  XACT,
complex*16, dimension( * )  X,
complex*16, dimension( * )  ARF,
complex*16, dimension( * )  ARFINV,
complex*16, dimension( * )  Z_WORK_ZLATMS,
complex*16, dimension( * )  Z_WORK_ZPOT02,
complex*16, dimension( * )  Z_WORK_ZPOT03,
double precision, dimension( * )  D_WORK_ZLATMS,
double precision, dimension( * )  D_WORK_ZLANHE,
double precision, dimension( * )  D_WORK_ZPOT01,
double precision, dimension( * )  D_WORK_ZPOT02,
double precision, dimension( * )  D_WORK_ZPOT03 
)

ZDRVRFP

Purpose:
 ZDRVRFP tests the LAPACK RFP routines:
     ZPFTRF, ZPFTRS, and ZPFTRI.

 This testing routine follow the same tests as ZDRVPO (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 ZTRTTF and
 ZTFTTR.

 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 ZPFTRF, 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 ZPFTRF 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 COMPLEX*16 array, dimension (NMAX*NMAX)
[out]ASAV
          ASAV is COMPLEX*16 array, dimension (NMAX*NMAX)
[out]AFAC
          AFAC is COMPLEX*16 array, dimension (NMAX*NMAX)
[out]AINV
          AINV is COMPLEX*16 array, dimension (NMAX*NMAX)
[out]B
          B is COMPLEX*16 array, dimension (NMAX*MAXRHS)
[out]BSAV
          BSAV is COMPLEX*16 array, dimension (NMAX*MAXRHS)
[out]XACT
          XACT is COMPLEX*16 array, dimension (NMAX*MAXRHS)
[out]X
          X is COMPLEX*16 array, dimension (NMAX*MAXRHS)
[out]ARF
          ARF is COMPLEX*16 array, dimension ((NMAX*(NMAX+1))/2)
[out]ARFINV
          ARFINV is COMPLEX*16 array, dimension ((NMAX*(NMAX+1))/2)
[out]Z_WORK_ZLATMS
          Z_WORK_ZLATMS is COMPLEX*16 array, dimension ( 3*NMAX )
[out]Z_WORK_ZPOT02
          Z_WORK_ZPOT02 is COMPLEX*16 array, dimension ( NMAX*MAXRHS )
[out]Z_WORK_ZPOT03
          Z_WORK_ZPOT03 is COMPLEX*16 array, dimension ( NMAX*NMAX )
[out]D_WORK_ZLATMS
          D_WORK_ZLATMS is DOUBLE PRECISION array, dimension ( NMAX )
[out]D_WORK_ZLANHE
          D_WORK_ZLANHE is DOUBLE PRECISION array, dimension ( NMAX )
[out]D_WORK_ZPOT01
          D_WORK_ZPOT01 is DOUBLE PRECISION array, dimension ( NMAX )
[out]D_WORK_ZPOT02
          D_WORK_ZPOT02 is DOUBLE PRECISION array, dimension ( NMAX )
[out]D_WORK_ZPOT03
          D_WORK_ZPOT03 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 246 of file zdrvrfp.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: