LAPACK 3.12.0
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,
303 + claipd, clarhs, clatb4, clatms, cpftri, cpftrf,
304 + cpftrs, cpot01, cpot02, cpot03, cpotri, cpotrf,
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*
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
clarhs
subroutine clarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
CLARHS
Definition clarhs.f:208
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
cget04
subroutine cget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
CGET04
Definition cget04.f:102
claipd
subroutine claipd(n, a, inda, vinda)
CLAIPD
Definition claipd.f:83
clatb4
subroutine clatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
CLATB4
Definition clatb4.f:121
clatms
subroutine clatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
CLATMS
Definition clatms.f:332
cpot01
subroutine cpot01(uplo, n, a, lda, afac, ldafac, rwork, resid)
CPOT01
Definition cpot01.f:106
cpot02
subroutine cpot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
CPOT02
Definition cpot02.f:127
cpot03
subroutine cpot03(uplo, n, a, lda, ainv, ldainv, work, ldwork, rwork, rcond, resid)
CPOT03
Definition cpot03.f:126
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
clanhe
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:124
cpftrf
subroutine cpftrf(transr, uplo, n, a, info)
CPFTRF
Definition cpftrf.f:211
cpftri
subroutine cpftri(transr, uplo, n, a, info)
CPFTRI
Definition cpftri.f:212
cpftrs
subroutine cpftrs(transr, uplo, n, nrhs, a, b, ldb, info)
CPFTRS
Definition cpftrs.f:220
cpotrf
subroutine cpotrf(uplo, n, a, lda, info)
CPOTRF
Definition cpotrf.f:107
cpotri
subroutine cpotri(uplo, n, a, lda, info)
CPOTRI
Definition cpotri.f:95
ctfttr
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:216
ctrttf
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:216
Here is the call graph for this function:
Here is the caller graph for this function: