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

◆ sdrvrfp()

subroutine sdrvrfp ( integer  nout,
integer  nn,
integer, dimension( nn )  nval,
integer  nns,
integer, dimension( nns )  nsval,
integer  nnt,
integer, dimension( nnt )  ntval,
real  thresh,
real, dimension( * )  a,
real, dimension( * )  asav,
real, dimension( * )  afac,
real, dimension( * )  ainv,
real, dimension( * )  b,
real, dimension( * )  bsav,
real, dimension( * )  xact,
real, dimension( * )  x,
real, dimension( * )  arf,
real, dimension( * )  arfinv,
real, dimension( * )  s_work_slatms,
real, dimension( * )  s_work_spot01,
real, dimension( * )  s_temp_spot02,
real, dimension( * )  s_temp_spot03,
real, dimension( * )  s_work_slansy,
real, dimension( * )  s_work_spot02,
real, dimension( * )  s_work_spot03 
)

SDRVRFP

Purpose:
 SDRVRFP tests the LAPACK RFP routines:
     SPFTRF, SPFTRS, and SPFTRI.

 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 STRTTF and
 STFTTR.

 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 SPFTRF, 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 SPFTRF 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 REAL array, dimension (NMAX*NMAX)
[out]ASAV
          ASAV is REAL array, dimension (NMAX*NMAX)
[out]AFAC
          AFAC is REAL array, dimension (NMAX*NMAX)
[out]AINV
          AINV is REAL array, dimension (NMAX*NMAX)
[out]B
          B is REAL array, dimension (NMAX*MAXRHS)
[out]BSAV
          BSAV is REAL array, dimension (NMAX*MAXRHS)
[out]XACT
          XACT is REAL array, dimension (NMAX*MAXRHS)
[out]X
          X is REAL array, dimension (NMAX*MAXRHS)
[out]ARF
          ARF is REAL array, dimension ((NMAX*(NMAX+1))/2)
[out]ARFINV
          ARFINV is REAL array, dimension ((NMAX*(NMAX+1))/2)
[out]S_WORK_SLATMS
          S_WORK_SLATMS is REAL array, dimension ( 3*NMAX )
[out]S_WORK_SPOT01
          S_WORK_SPOT01 is REAL array, dimension ( NMAX )
[out]S_TEMP_SPOT02
          S_TEMP_SPOT02 is REAL array, dimension ( NMAX*MAXRHS )
[out]S_TEMP_SPOT03
          S_TEMP_SPOT03 is REAL array, dimension ( NMAX*NMAX )
[out]S_WORK_SLANSY
          S_WORK_SLANSY is REAL array, dimension ( NMAX )
[out]S_WORK_SPOT02
          S_WORK_SPOT02 is REAL array, dimension ( NMAX )
[out]S_WORK_SPOT03
          S_WORK_SPOT03 is REAL array, dimension ( NMAX )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 232 of file sdrvrfp.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 REAL THRESH
246* ..
247* .. Array Arguments ..
248 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
249 REAL A( * )
250 REAL AINV( * )
251 REAL ASAV( * )
252 REAL B( * )
253 REAL BSAV( * )
254 REAL AFAC( * )
255 REAL ARF( * )
256 REAL ARFINV( * )
257 REAL XACT( * )
258 REAL X( * )
259 REAL S_WORK_SLATMS( * )
260 REAL S_WORK_SPOT01( * )
261 REAL S_TEMP_SPOT02( * )
262 REAL S_TEMP_SPOT03( * )
263 REAL S_WORK_SLANSY( * )
264 REAL S_WORK_SPOT02( * )
265 REAL S_WORK_SPOT03( * )
266* ..
267*
268* =====================================================================
269*
270* .. Parameters ..
271 REAL ONE, ZERO
272 parameter( one = 1.0e+0, zero = 0.0e+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 REAL ANORM, AINVNM, CNDNUM, RCONDC
284* ..
285* .. Local Arrays ..
286 CHARACTER UPLOS( 2 ), FORMS( 2 )
287 INTEGER ISEED( 4 ), ISEEDY( 4 )
288 REAL RESULT( NTESTS )
289* ..
290* .. External Functions ..
291 REAL SLANSY
292 EXTERNAL slansy
293* ..
294* .. External Subroutines ..
295 EXTERNAL aladhd, alaerh, alasvm, sget04, stfttr, slacpy,
296 + slarhs, slatb4, slatms, spftri, spftrf, spftrs,
297 + spot01, spot02, spot03, spotri, spotrf, strttf
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 SLATB4 and generate a test
355* matrix with SLATMS.
356*
357 CALL slatb4( 'SPO', imat, n, n, ctype, kl, ku,
358 + anorm, mode, cndnum, dist )
359*
360 srnamt = 'SLATMS'
361 CALL slatms( n, n, dist, iseed, ctype,
362 + s_work_slatms,
363 + mode, cndnum, anorm, kl, ku, uplo, a,
364 + lda, s_work_slatms, info )
365*
366* Check error code from SLATMS.
367*
368 IF( info.NE.0 ) THEN
369 CALL alaerh( 'SPF', 'SLATMS', 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 slacpy( 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 = slansy( '1', uplo, n, a, lda,
428 + s_work_slansy )
429*
430* Factor the matrix A.
431*
432 CALL spotrf( uplo, n, a, lda, info )
433*
434* Form the inverse of A.
435*
436 CALL spotri( 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 = slansy( '1', uplo, n, a, lda,
443 + s_work_slansy )
444 rcondc = ( one / anorm ) / ainvnm
445*
446* Restore the matrix A.
447*
448 CALL slacpy( 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 = 'SLARHS'
456 CALL slarhs( 'SPO', 'N', uplo, ' ', n, n, kl, ku,
457 + nrhs, a, lda, xact, lda, b, lda,
458 + iseed, info )
459 CALL slacpy( '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 slacpy( uplo, n, n, a, lda, afac, lda )
465 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldb )
466*
467 srnamt = 'STRTTF'
468 CALL strttf( cform, uplo, n, afac, lda, arf, info )
469 srnamt = 'SPFTRF'
470 CALL spftrf( cform, uplo, n, arf, info )
471*
472* Check error code from SPFTRF.
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( 'SPF', 'SPFSV ', 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 = 'SPFTRS'
493 CALL spftrs( cform, uplo, n, nrhs, arf, x, ldb,
494 + info )
495*
496 srnamt = 'STFTTR'
497 CALL stfttr( cform, uplo, n, arf, afac, lda, info )
498*
499* Reconstruct matrix from factors and compute
500* residual.
501*
502 CALL slacpy( uplo, n, n, afac, lda, asav, lda )
503 CALL spot01( uplo, n, a, lda, afac, lda,
504 + s_work_spot01, result( 1 ) )
505 CALL slacpy( 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 slacpy( 'A', n+1, n/2, arf, n+1, arfinv,
511 + n+1 )
512 ELSE
513 CALL slacpy( 'A', n, (n+1)/2, arf, n, arfinv,
514 + n )
515 END IF
516*
517 srnamt = 'SPFTRI'
518 CALL spftri( cform, uplo, n, arfinv , info )
519*
520 srnamt = 'STFTTR'
521 CALL stfttr( cform, uplo, n, arfinv, ainv, lda,
522 + info )
523*
524* Check error code from SPFTRI.
525*
526 IF( info.NE.0 )
527 + CALL alaerh( 'SPO', 'SPFTRI', info, 0, uplo, n,
528 + n, -1, -1, -1, imat, nfail, nerrs,
529 + nout )
530*
531 CALL spot03( uplo, n, a, lda, ainv, lda,
532 + s_temp_spot03, lda, s_work_spot03,
533 + rcondc, result( 2 ) )
534*
535* Compute residual of the computed solution.
536*
537 CALL slacpy( 'Full', n, nrhs, b, lda,
538 + s_temp_spot02, lda )
539 CALL spot02( uplo, n, nrhs, a, lda, x, lda,
540 + s_temp_spot02, lda, s_work_spot02,
541 + result( 3 ) )
542*
543* Check solution from generated exact solution.
544
545 CALL sget04( 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, 'SPF' )
556 WRITE( nout, fmt = 9999 )'SPFSV ', 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( 'SPF', 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 SDRVRFP
578*
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
slarhs
subroutine slarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
SLARHS
Definition slarhs.f:205
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
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
slansy
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
spftrf
subroutine spftrf(transr, uplo, n, a, info)
SPFTRF
Definition spftrf.f:198
spftri
subroutine spftri(transr, uplo, n, a, info)
SPFTRI
Definition spftri.f:191
spftrs
subroutine spftrs(transr, uplo, n, nrhs, a, b, ldb, info)
SPFTRS
Definition spftrs.f:199
spotrf
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:107
spotri
subroutine spotri(uplo, n, a, lda, info)
SPOTRI
Definition spotri.f:95
stfttr
subroutine stfttr(transr, uplo, n, arf, a, lda, info)
STFTTR copies a triangular matrix from the rectangular full packed format (TF) to the standard full f...
Definition stfttr.f:196
strttf
subroutine strttf(transr, uplo, n, a, lda, arf, info)
STRTTF copies a triangular matrix from the standard full format (TR) to the rectangular full packed f...
Definition strttf.f:194
sget04
subroutine sget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
SGET04
Definition sget04.f:102
slatb4
subroutine slatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
SLATB4
Definition slatb4.f:120
slatms
subroutine slatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
SLATMS
Definition slatms.f:321
spot01
subroutine spot01(uplo, n, a, lda, afac, ldafac, rwork, resid)
SPOT01
Definition spot01.f:104
spot02
subroutine spot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
SPOT02
Definition spot02.f:127
spot03
subroutine spot03(uplo, n, a, lda, ainv, ldainv, work, ldwork, rwork, rcond, resid)
SPOT03
Definition spot03.f:125
Here is the call graph for this function:
Here is the caller graph for this function: