LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, convertion 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_SLATMS
          S_WORK_SLATMS is REAL array, dimension ( 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.
Date
November 2013

Definition at line 245 of file sdrvrfp.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  REAL thresh
254 * ..
255 * .. Array Arguments ..
256  INTEGER nval( nn ), nsval( nns ), ntval( nnt )
257  REAL a( * )
258  REAL ainv( * )
259  REAL asav( * )
260  REAL b( * )
261  REAL bsav( * )
262  REAL afac( * )
263  REAL arf( * )
264  REAL arfinv( * )
265  REAL xact( * )
266  REAL x( * )
267  REAL s_work_slatms( * )
268  REAL s_work_spot01( * )
269  REAL s_temp_spot02( * )
270  REAL s_temp_spot03( * )
271  REAL s_work_slansy( * )
272  REAL s_work_spot02( * )
273  REAL s_work_spot03( * )
274 * ..
275 *
276 * =====================================================================
277 *
278 * .. Parameters ..
279  REAL one, zero
280  parameter ( one = 1.0e+0, zero = 0.0e+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  REAL anorm, ainvnm, cndnum, rcondc
292 * ..
293 * .. Local Arrays ..
294  CHARACTER uplos( 2 ), forms( 2 )
295  INTEGER iseed( 4 ), iseedy( 4 )
296  REAL result( ntests )
297 * ..
298 * .. External Functions ..
299  REAL slansy
300  EXTERNAL slansy
301 * ..
302 * .. External Subroutines ..
303  EXTERNAL aladhd, alaerh, alasvm, sget04, stfttr, slacpy,
304  + slarhs, slatb4, slatms, spftri, spftrf, spftrs,
305  + spot01, spot02, spot03, spotri, spotrf, strttf
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 SLATB4 and generate a test
363 * matrix with SLATMS.
364 *
365  CALL slatb4( 'SPO', imat, n, n, ctype, kl, ku,
366  + anorm, mode, cndnum, dist )
367 *
368  srnamt = 'SLATMS'
369  CALL slatms( n, n, dist, iseed, ctype,
370  + s_work_slatms,
371  + mode, cndnum, anorm, kl, ku, uplo, a,
372  + lda, s_work_slatms, info )
373 *
374 * Check error code from SLATMS.
375 *
376  IF( info.NE.0 ) THEN
377  CALL alaerh( 'SPF', 'SLATMS', 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 slacpy( 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 = slansy( '1', uplo, n, a, lda,
436  + s_work_slansy )
437 *
438 * Factor the matrix A.
439 *
440  CALL spotrf( uplo, n, a, lda, info )
441 *
442 * Form the inverse of A.
443 *
444  CALL spotri( uplo, n, a, lda, info )
445 
446  IF ( n .NE. 0 ) THEN
447 *
448 * Compute the 1-norm condition number of A.
449 *
450  ainvnm = slansy( '1', uplo, n, a, lda,
451  + s_work_slansy )
452  rcondc = ( one / anorm ) / ainvnm
453 *
454 * Restore the matrix A.
455 *
456  CALL slacpy( uplo, n, n, asav, lda, a, lda )
457  END IF
458 *
459  END IF
460 *
461 * Form an exact solution and set the right hand side.
462 *
463  srnamt = 'SLARHS'
464  CALL slarhs( 'SPO', 'N', uplo, ' ', n, n, kl, ku,
465  + nrhs, a, lda, xact, lda, b, lda,
466  + iseed, info )
467  CALL slacpy( 'Full', n, nrhs, b, lda, bsav, lda )
468 *
469 * Compute the L*L' or U'*U factorization of the
470 * matrix and solve the system.
471 *
472  CALL slacpy( uplo, n, n, a, lda, afac, lda )
473  CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldb )
474 *
475  srnamt = 'STRTTF'
476  CALL strttf( cform, uplo, n, afac, lda, arf, info )
477  srnamt = 'SPFTRF'
478  CALL spftrf( cform, uplo, n, arf, info )
479 *
480 * Check error code from SPFTRF.
481 *
482  IF( info.NE.izero ) THEN
483 *
484 * LANGOU: there is a small hick here: IZERO should
485 * always be INFO however if INFO is ZERO, ALAERH does not
486 * complain.
487 *
488  CALL alaerh( 'SPF', 'SPFSV ', info, izero,
489  + uplo, n, n, -1, -1, nrhs, iit,
490  + nfail, nerrs, nout )
491  GO TO 100
492  END IF
493 *
494 * Skip the tests if INFO is not 0.
495 *
496  IF( info.NE.0 ) THEN
497  GO TO 100
498  END IF
499 *
500  srnamt = 'SPFTRS'
501  CALL spftrs( cform, uplo, n, nrhs, arf, x, ldb,
502  + info )
503 *
504  srnamt = 'STFTTR'
505  CALL stfttr( cform, uplo, n, arf, afac, lda, info )
506 *
507 * Reconstruct matrix from factors and compute
508 * residual.
509 *
510  CALL slacpy( uplo, n, n, afac, lda, asav, lda )
511  CALL spot01( uplo, n, a, lda, afac, lda,
512  + s_work_spot01, result( 1 ) )
513  CALL slacpy( uplo, n, n, asav, lda, afac, lda )
514 *
515 * Form the inverse and compute the residual.
516 *
517  IF(mod(n,2).EQ.0)THEN
518  CALL slacpy( 'A', n+1, n/2, arf, n+1, arfinv,
519  + n+1 )
520  ELSE
521  CALL slacpy( 'A', n, (n+1)/2, arf, n, arfinv,
522  + n )
523  END IF
524 *
525  srnamt = 'SPFTRI'
526  CALL spftri( cform, uplo, n, arfinv , info )
527 *
528  srnamt = 'STFTTR'
529  CALL stfttr( cform, uplo, n, arfinv, ainv, lda,
530  + info )
531 *
532 * Check error code from SPFTRI.
533 *
534  IF( info.NE.0 )
535  + CALL alaerh( 'SPO', 'SPFTRI', info, 0, uplo, n,
536  + n, -1, -1, -1, imat, nfail, nerrs,
537  + nout )
538 *
539  CALL spot03( uplo, n, a, lda, ainv, lda,
540  + s_temp_spot03, lda, s_work_spot03,
541  + rcondc, result( 2 ) )
542 *
543 * Compute residual of the computed solution.
544 *
545  CALL slacpy( 'Full', n, nrhs, b, lda,
546  + s_temp_spot02, lda )
547  CALL spot02( uplo, n, nrhs, a, lda, x, lda,
548  + s_temp_spot02, lda, s_work_spot02,
549  + result( 3 ) )
550 *
551 * Check solution from generated exact solution.
552 
553  CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
554  + result( 4 ) )
555  nt = 4
556 *
557 * Print information about the tests that did not
558 * pass the threshold.
559 *
560  DO 60 k = 1, nt
561  IF( result( k ).GE.thresh ) THEN
562  IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
563  + CALL aladhd( nout, 'SPF' )
564  WRITE( nout, fmt = 9999 )'SPFSV ', uplo,
565  + n, iit, k, result( k )
566  nfail = nfail + 1
567  END IF
568  60 CONTINUE
569  nrun = nrun + nt
570  100 CONTINUE
571  110 CONTINUE
572  120 CONTINUE
573  980 CONTINUE
574  130 CONTINUE
575 *
576 * Print a summary of the results.
577 *
578  CALL alasvm( 'SPF', nout, nfail, nrun, nerrs )
579 *
580  9999 FORMAT( 1x, a6, ', UPLO=''', a1, ''', N =', i5, ', type ', i1,
581  + ', test(', i1, ')=', g12.5 )
582 *
583  RETURN
584 *
585 * End of SDRVRFP
586 *
alasvm
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
alaerh
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:149
slatb4
subroutine slatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
SLATB4
Definition: slatb4.f:122
slarhs
subroutine slarhs(PATH, XTYPE, UPLO, TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, LDB, ISEED, INFO)
SLARHS
Definition: slarhs.f:206
spot01
subroutine spot01(UPLO, N, A, LDA, AFAC, LDAFAC, RWORK, RESID)
SPOT01
Definition: spot01.f:106
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:105
slatms
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:323
spotrf
subroutine spotrf(UPLO, N, A, LDA, INFO)
SPOTRF
Definition: spotrf.f:109
spftrf
subroutine spftrf(TRANSR, UPLO, N, A, INFO)
SPFTRF
Definition: spftrf.f:200
sget04
subroutine sget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
SGET04
Definition: sget04.f:104
spot03
subroutine spot03(UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK, RCOND, RESID)
SPOT03
Definition: spot03.f:127
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:196
aladhd
subroutine aladhd(IOUNIT, PATH)
ALADHD
Definition: aladhd.f:80
spftri
subroutine spftri(TRANSR, UPLO, N, A, INFO)
SPFTRI
Definition: spftri.f:193
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:198
spot02
subroutine spot02(UPLO, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
SPOT02
Definition: spot02.f:129
spftrs
subroutine spftrs(TRANSR, UPLO, N, NRHS, A, B, LDB, INFO)
SPFTRS
Definition: spftrs.f:201
spotri
subroutine spotri(UPLO, N, A, LDA, INFO)
SPOTRI
Definition: spotri.f:97
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, or the element of largest absolute value of a real symmetric matrix.
Definition: slansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function: