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

◆ cdrvrfp()

subroutine cdrvrfp ( integer nout,
integer nn,
integer, dimension( nn ) nval,
integer nns,
integer, dimension( nns ) nsval,
integer nnt,
integer, dimension( nnt ) ntval,
real thresh,
complex, dimension( * ) a,
complex, dimension( * ) asav,
complex, dimension( * ) afac,
complex, dimension( * ) ainv,
complex, dimension( * ) b,
complex, dimension( * ) bsav,
complex, dimension( * ) xact,
complex, dimension( * ) x,
complex, dimension( * ) arf,
complex, dimension( * ) arfinv,
complex, dimension( * ) c_work_clatms,
complex, dimension( * ) c_work_cpot02,
complex, dimension( * ) c_work_cpot03,
real, dimension( * ) s_work_clatms,
real, dimension( * ) s_work_clanhe,
real, dimension( * ) s_work_cpot01,
real, dimension( * ) s_work_cpot02,
real, dimension( * ) s_work_cpot03 )

CDRVRFP

Purpose:
!>
!> CDRVRFP tests the LAPACK RFP routines:
!>     CPFTRF, CPFTRS, and CPFTRI.
!>
!> This testing routine follow the same tests as CDRVPO (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 CTRTTF and
!> CTFTTR.
!>
!> 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 CPFTRF, 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 CPFTRF 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 REAL
!>                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 array, dimension (NMAX*NMAX)
!> 
[out]ASAV
!>          ASAV is COMPLEX array, dimension (NMAX*NMAX)
!> 
[out]AFAC
!>          AFAC is COMPLEX array, dimension (NMAX*NMAX)
!> 
[out]AINV
!>          AINV is COMPLEX array, dimension (NMAX*NMAX)
!> 
[out]B
!>          B is COMPLEX array, dimension (NMAX*MAXRHS)
!> 
[out]BSAV
!>          BSAV is COMPLEX array, dimension (NMAX*MAXRHS)
!> 
[out]XACT
!>          XACT is COMPLEX array, dimension (NMAX*MAXRHS)
!> 
[out]X
!>          X is COMPLEX array, dimension (NMAX*MAXRHS)
!> 
[out]ARF
!>          ARF is COMPLEX array, dimension ((NMAX*(NMAX+1))/2)
!> 
[out]ARFINV
!>          ARFINV is COMPLEX array, dimension ((NMAX*(NMAX+1))/2)
!> 
[out]C_WORK_CLATMS
!>          C_WORK_CLATMS is COMPLEX array, dimension ( 3*NMAX )
!> 
[out]C_WORK_CPOT02
!>          C_WORK_CPOT02 is COMPLEX array, dimension ( NMAX*MAXRHS )
!> 
[out]C_WORK_CPOT03
!>          C_WORK_CPOT03 is COMPLEX array, dimension ( NMAX*NMAX )
!> 
[out]S_WORK_CLATMS
!>          S_WORK_CLATMS is REAL array, dimension ( NMAX )
!> 
[out]S_WORK_CLANHE
!>          S_WORK_CLANHE is REAL array, dimension ( NMAX )
!> 
[out]S_WORK_CPOT01
!>          S_WORK_CPOT01 is REAL array, dimension ( NMAX )
!> 
[out]S_WORK_CPOT02
!>          S_WORK_CPOT02 is REAL array, dimension ( NMAX )
!> 
[out]S_WORK_CPOT03
!>          S_WORK_CPOT03 is REAL array, dimension ( NMAX )
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 238 of file cdrvrfp.f.

244*
245* -- LAPACK test routine --
246* -- LAPACK is a software package provided by Univ. of Tennessee, --
247* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
248*
249* .. Scalar Arguments ..
250 INTEGER NN, NNS, NNT, NOUT
251 REAL THRESH
252* ..
253* .. Array Arguments ..
254 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
255 COMPLEX A( * )
256 COMPLEX AINV( * )
257 COMPLEX ASAV( * )
258 COMPLEX B( * )
259 COMPLEX BSAV( * )
260 COMPLEX AFAC( * )
261 COMPLEX ARF( * )
262 COMPLEX ARFINV( * )
263 COMPLEX XACT( * )
264 COMPLEX X( * )
265 COMPLEX C_WORK_CLATMS( * )
266 COMPLEX C_WORK_CPOT02( * )
267 COMPLEX C_WORK_CPOT03( * )
268 REAL S_WORK_CLATMS( * )
269 REAL S_WORK_CLANHE( * )
270 REAL S_WORK_CPOT01( * )
271 REAL S_WORK_CPOT02( * )
272 REAL S_WORK_CPOT03( * )
273* ..
274*
275* =====================================================================
276*
277* .. Parameters ..
278 REAL ONE, ZERO
279 parameter( one = 1.0e+0, zero = 0.0e+0 )
280 INTEGER NTESTS
281 parameter( ntests = 4 )
282* ..
283* .. Local Scalars ..
284 LOGICAL ZEROT
285 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL,
286 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN,
287 + IIT, IIS
288 CHARACTER DIST, CTYPE, UPLO, CFORM
289 INTEGER KL, KU, MODE
290 REAL ANORM, AINVNM, CNDNUM, RCONDC
291* ..
292* .. Local Arrays ..
293 CHARACTER UPLOS( 2 ), FORMS( 2 )
294 INTEGER ISEED( 4 ), ISEEDY( 4 )
295 REAL RESULT( NTESTS )
296* ..
297* .. External Functions ..
298 REAL CLANHE
299 EXTERNAL clanhe
300* ..
301* .. External Subroutines ..
302 EXTERNAL aladhd, alaerh, alasvm, cget04, ctfttr, clacpy,
305 + ctrttf
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', 'C' /
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 CLATB4 and generate a test
363* matrix with CLATMS.
364*
365 CALL clatb4( 'CPO', imat, n, n, ctype, kl, ku,
366 + anorm, mode, cndnum, dist )
367*
368 srnamt = 'CLATMS'
369 CALL clatms( n, n, dist, iseed, ctype,
370 + s_work_clatms,
371 + mode, cndnum, anorm, kl, ku, uplo, a,
372 + lda, c_work_clatms, info )
373*
374* Check error code from CLATMS.
375*
376 IF( info.NE.0 ) THEN
377 CALL alaerh( 'CPF', 'CLATMS', 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* Set the imaginary part of the diagonals.
424*
425 CALL claipd( n, a, lda+1, 0 )
426*
427* Save a copy of the matrix A in ASAV.
428*
429 CALL clacpy( uplo, n, n, a, lda, asav, lda )
430*
431* Compute the condition number of A (RCONDC).
432*
433 IF( zerot ) THEN
434 rcondc = zero
435 ELSE
436*
437* Compute the 1-norm of A.
438*
439 anorm = clanhe( '1', uplo, n, a, lda,
440 + s_work_clanhe )
441*
442* Factor the matrix A.
443*
444 CALL cpotrf( uplo, n, a, lda, info )
445*
446* Form the inverse of A.
447*
448 CALL cpotri( uplo, n, a, lda, info )
449
450 IF ( n .NE. 0 ) THEN
451*
452* Compute the 1-norm condition number of A.
453*
454 ainvnm = clanhe( '1', uplo, n, a, lda,
455 + s_work_clanhe )
456 rcondc = ( one / anorm ) / ainvnm
457*
458* Restore the matrix A.
459*
460 CALL clacpy( uplo, n, n, asav, lda, a, lda )
461 END IF
462*
463 END IF
464*
465* Form an exact solution and set the right hand side.
466*
467 srnamt = 'CLARHS'
468 CALL clarhs( 'CPO', 'N', uplo, ' ', n, n, kl, ku,
469 + nrhs, a, lda, xact, lda, b, lda,
470 + iseed, info )
471 CALL clacpy( 'Full', n, nrhs, b, lda, bsav, lda )
472*
473* Compute the L*L' or U'*U factorization of the
474* matrix and solve the system.
475*
476 CALL clacpy( uplo, n, n, a, lda, afac, lda )
477 CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldb )
478*
479 srnamt = 'CTRTTF'
480 CALL ctrttf( cform, uplo, n, afac, lda, arf, info )
481 srnamt = 'CPFTRF'
482 CALL cpftrf( cform, uplo, n, arf, info )
483*
484* Check error code from CPFTRF.
485*
486 IF( info.NE.izero ) THEN
487*
488* LANGOU: there is a small hick here: IZERO should
489* always be INFO however if INFO is ZERO, ALAERH does not
490* complain.
491*
492 CALL alaerh( 'CPF', 'CPFSV ', info, izero,
493 + uplo, n, n, -1, -1, nrhs, iit,
494 + nfail, nerrs, nout )
495 GO TO 100
496 END IF
497*
498* Skip the tests if INFO is not 0.
499*
500 IF( info.NE.0 ) THEN
501 GO TO 100
502 END IF
503*
504 srnamt = 'CPFTRS'
505 CALL cpftrs( cform, uplo, n, nrhs, arf, x, ldb,
506 + info )
507*
508 srnamt = 'CTFTTR'
509 CALL ctfttr( cform, uplo, n, arf, afac, lda, info )
510*
511* Reconstruct matrix from factors and compute
512* residual.
513*
514 CALL clacpy( uplo, n, n, afac, lda, asav, lda )
515 CALL cpot01( uplo, n, a, lda, afac, lda,
516 + s_work_cpot01, result( 1 ) )
517 CALL clacpy( uplo, n, n, asav, lda, afac, lda )
518*
519* Form the inverse and compute the residual.
520*
521 IF(mod(n,2).EQ.0)THEN
522 CALL clacpy( 'A', n+1, n/2, arf, n+1, arfinv,
523 + n+1 )
524 ELSE
525 CALL clacpy( 'A', n, (n+1)/2, arf, n, arfinv,
526 + n )
527 END IF
528*
529 srnamt = 'CPFTRI'
530 CALL cpftri( cform, uplo, n, arfinv , info )
531*
532 srnamt = 'CTFTTR'
533 CALL ctfttr( cform, uplo, n, arfinv, ainv, lda,
534 + info )
535*
536* Check error code from CPFTRI.
537*
538 IF( info.NE.0 )
539 + CALL alaerh( 'CPO', 'CPFTRI', info, 0, uplo, n,
540 + n, -1, -1, -1, imat, nfail, nerrs,
541 + nout )
542*
543 CALL cpot03( uplo, n, a, lda, ainv, lda,
544 + c_work_cpot03, lda, s_work_cpot03,
545 + rcondc, result( 2 ) )
546*
547* Compute residual of the computed solution.
548*
549 CALL clacpy( 'Full', n, nrhs, b, lda,
550 + c_work_cpot02, lda )
551 CALL cpot02( uplo, n, nrhs, a, lda, x, lda,
552 + c_work_cpot02, lda, s_work_cpot02,
553 + result( 3 ) )
554*
555* Check solution from generated exact solution.
556*
557 CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
558 + result( 4 ) )
559 nt = 4
560*
561* Print information about the tests that did not
562* pass the threshold.
563*
564 DO 60 k = 1, nt
565 IF( result( k ).GE.thresh ) THEN
566 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
567 + CALL aladhd( nout, 'CPF' )
568 WRITE( nout, fmt = 9999 )'CPFSV ', uplo,
569 + n, iit, k, result( k )
570 nfail = nfail + 1
571 END IF
572 60 CONTINUE
573 nrun = nrun + nt
574 100 CONTINUE
575 110 CONTINUE
576 120 CONTINUE
577 980 CONTINUE
578 130 CONTINUE
579*
580* Print a summary of the results.
581*
582 CALL alasvm( 'CPF', nout, nfail, nrun, nerrs )
583*
584 9999 FORMAT( 1x, a6, ', UPLO=''', a1, ''', N =', i5, ', type ', i1,
585 + ', test(', i1, ')=', g12.5 )
586*
587 RETURN
588*
589* End of CDRVRFP
590*
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine clarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
CLARHS
Definition clarhs.f:208
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 cget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
CGET04
Definition cget04.f:102
subroutine claipd(n, a, inda, vinda)
CLAIPD
Definition claipd.f:83
subroutine clatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
CLATB4
Definition clatb4.f:121
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
Definition clatms.f:332
subroutine cpot01(uplo, n, a, lda, afac, ldafac, rwork, resid)
CPOT01
Definition cpot01.f:106
subroutine cpot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
CPOT02
Definition cpot02.f:127
subroutine cpot03(uplo, n, a, lda, ainv, ldainv, work, ldwork, rwork, rcond, resid)
CPOT03
Definition cpot03.f:126
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
real function clanhe(norm, uplo, n, a, lda, work)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhe.f:122
subroutine cpftrf(transr, uplo, n, a, info)
CPFTRF
Definition cpftrf.f:209
subroutine cpftri(transr, uplo, n, a, info)
CPFTRI
Definition cpftri.f:210
subroutine cpftrs(transr, uplo, n, nrhs, a, b, ldb, info)
CPFTRS
Definition cpftrs.f:218
subroutine cpotrf(uplo, n, a, lda, info)
CPOTRF
Definition cpotrf.f:105
subroutine cpotri(uplo, n, a, lda, info)
CPOTRI
Definition cpotri.f:93
subroutine ctfttr(transr, uplo, n, arf, a, lda, info)
CTFTTR copies a triangular matrix from the rectangular full packed format (TF) to the standard full f...
Definition ctfttr.f:214
subroutine ctrttf(transr, uplo, n, a, lda, arf, info)
CTRTTF copies a triangular matrix from the standard full format (TR) to the rectangular full packed f...
Definition ctrttf.f:214
Here is the call graph for this function:
Here is the caller graph for this function: