LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sdrvrfp.f
Go to the documentation of this file.
1*> \brief \b SDRVRFP
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE SDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL,
12* + THRESH, A, ASAV, AFAC, AINV, B,
13* + BSAV, XACT, X, ARF, ARFINV,
14* + S_WORK_SLATMS, S_WORK_SPOT01, S_TEMP_SPOT02,
15* + S_TEMP_SPOT03, S_WORK_SLANSY,
16* + S_WORK_SPOT02, S_WORK_SPOT03 )
17*
18* .. Scalar Arguments ..
19* INTEGER NN, NNS, NNT, NOUT
20* REAL THRESH
21* ..
22* .. Array Arguments ..
23* INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
24* REAL A( * )
25* REAL AINV( * )
26* REAL ASAV( * )
27* REAL B( * )
28* REAL BSAV( * )
29* REAL AFAC( * )
30* REAL ARF( * )
31* REAL ARFINV( * )
32* REAL XACT( * )
33* REAL X( * )
34* REAL S_WORK_SLATMS( * )
35* REAL S_WORK_SPOT01( * )
36* REAL S_TEMP_SPOT02( * )
37* REAL S_TEMP_SPOT03( * )
38* REAL S_WORK_SLANSY( * )
39* REAL S_WORK_SPOT02( * )
40* REAL S_WORK_SPOT03( * )
41* ..
42*
43*
44*> \par Purpose:
45* =============
46*>
47*> \verbatim
48*>
49*> SDRVRFP tests the LAPACK RFP routines:
50*> SPFTRF, SPFTRS, and SPFTRI.
51*>
52*> This testing routine follow the same tests as DDRVPO (test for the full
53*> format Symmetric Positive Definite solver).
54*>
55*> The tests are performed in Full Format, conversion back and forth from
56*> full format to RFP format are performed using the routines STRTTF and
57*> STFTTR.
58*>
59*> First, a specific matrix A of size N is created. There is nine types of
60*> different matrixes possible.
61*> 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS)
62*> 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS
63*> *3. First row and column zero 8. Scaled near underflow
64*> *4. Last row and column zero 9. Scaled near overflow
65*> *5. Middle row and column zero
66*> (* - tests error exits from SPFTRF, no test ratios are computed)
67*> A solution XACT of size N-by-NRHS is created and the associated right
68*> hand side B as well. Then SPFTRF is called to compute L (or U), the
69*> Cholesky factor of A. Then L (or U) is used to solve the linear system
70*> of equations AX = B. This gives X. Then L (or U) is used to compute the
71*> inverse of A, AINV. The following four tests are then performed:
72*> (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or
73*> norm( U'*U - A ) / ( N * norm(A) * EPS ),
74*> (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
75*> (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
76*> (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
77*> where EPS is the machine precision, RCOND the condition number of A, and
78*> norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4).
79*> Errors occur when INFO parameter is not as expected. Failures occur when
80*> a test ratios is greater than THRES.
81*> \endverbatim
82*
83* Arguments:
84* ==========
85*
86*> \param[in] NOUT
87*> \verbatim
88*> NOUT is INTEGER
89*> The unit number for output.
90*> \endverbatim
91*>
92*> \param[in] NN
93*> \verbatim
94*> NN is INTEGER
95*> The number of values of N contained in the vector NVAL.
96*> \endverbatim
97*>
98*> \param[in] NVAL
99*> \verbatim
100*> NVAL is INTEGER array, dimension (NN)
101*> The values of the matrix dimension N.
102*> \endverbatim
103*>
104*> \param[in] NNS
105*> \verbatim
106*> NNS is INTEGER
107*> The number of values of NRHS contained in the vector NSVAL.
108*> \endverbatim
109*>
110*> \param[in] NSVAL
111*> \verbatim
112*> NSVAL is INTEGER array, dimension (NNS)
113*> The values of the number of right-hand sides NRHS.
114*> \endverbatim
115*>
116*> \param[in] NNT
117*> \verbatim
118*> NNT is INTEGER
119*> The number of values of MATRIX TYPE contained in the vector NTVAL.
120*> \endverbatim
121*>
122*> \param[in] NTVAL
123*> \verbatim
124*> NTVAL is INTEGER array, dimension (NNT)
125*> The values of matrix type (between 0 and 9 for PO/PP/PF matrices).
126*> \endverbatim
127*>
128*> \param[in] THRESH
129*> \verbatim
130*> THRESH is REAL
131*> The threshold value for the test ratios. A result is
132*> included in the output file if RESULT >= THRESH. To have
133*> every test ratio printed, use THRESH = 0.
134*> \endverbatim
135*>
136*> \param[out] A
137*> \verbatim
138*> A is REAL array, dimension (NMAX*NMAX)
139*> \endverbatim
140*>
141*> \param[out] ASAV
142*> \verbatim
143*> ASAV is REAL array, dimension (NMAX*NMAX)
144*> \endverbatim
145*>
146*> \param[out] AFAC
147*> \verbatim
148*> AFAC is REAL array, dimension (NMAX*NMAX)
149*> \endverbatim
150*>
151*> \param[out] AINV
152*> \verbatim
153*> AINV is REAL array, dimension (NMAX*NMAX)
154*> \endverbatim
155*>
156*> \param[out] B
157*> \verbatim
158*> B is REAL array, dimension (NMAX*MAXRHS)
159*> \endverbatim
160*>
161*> \param[out] BSAV
162*> \verbatim
163*> BSAV is REAL array, dimension (NMAX*MAXRHS)
164*> \endverbatim
165*>
166*> \param[out] XACT
167*> \verbatim
168*> XACT is REAL array, dimension (NMAX*MAXRHS)
169*> \endverbatim
170*>
171*> \param[out] X
172*> \verbatim
173*> X is REAL array, dimension (NMAX*MAXRHS)
174*> \endverbatim
175*>
176*> \param[out] ARF
177*> \verbatim
178*> ARF is REAL array, dimension ((NMAX*(NMAX+1))/2)
179*> \endverbatim
180*>
181*> \param[out] ARFINV
182*> \verbatim
183*> ARFINV is REAL array, dimension ((NMAX*(NMAX+1))/2)
184*> \endverbatim
185*>
186*> \param[out] S_WORK_SLATMS
187*> \verbatim
188*> S_WORK_SLATMS is REAL array, dimension ( 3*NMAX )
189*> \endverbatim
190*>
191*> \param[out] S_WORK_SPOT01
192*> \verbatim
193*> S_WORK_SPOT01 is REAL array, dimension ( NMAX )
194*> \endverbatim
195*>
196*> \param[out] S_TEMP_SPOT02
197*> \verbatim
198*> S_TEMP_SPOT02 is REAL array, dimension ( NMAX*MAXRHS )
199*> \endverbatim
200*>
201*> \param[out] S_TEMP_SPOT03
202*> \verbatim
203*> S_TEMP_SPOT03 is REAL array, dimension ( NMAX*NMAX )
204*> \endverbatim
205*>
206*> \param[out] S_WORK_SLANSY
207*> \verbatim
208*> S_WORK_SLANSY is REAL array, dimension ( NMAX )
209*> \endverbatim
210*>
211*> \param[out] S_WORK_SPOT02
212*> \verbatim
213*> S_WORK_SPOT02 is REAL array, dimension ( NMAX )
214*> \endverbatim
215*>
216*> \param[out] S_WORK_SPOT03
217*> \verbatim
218*> S_WORK_SPOT03 is REAL array, dimension ( NMAX )
219*> \endverbatim
220*
221* Authors:
222* ========
223*
224*> \author Univ. of Tennessee
225*> \author Univ. of California Berkeley
226*> \author Univ. of Colorado Denver
227*> \author NAG Ltd.
228*
229*> \ingroup single_lin
230*
231* =====================================================================
232 SUBROUTINE sdrvrfp( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL,
233 + THRESH, A, ASAV, AFAC, AINV, B,
234 + BSAV, XACT, X, ARF, ARFINV,
235 + S_WORK_SLATMS, S_WORK_SPOT01, S_TEMP_SPOT02,
236 + S_TEMP_SPOT03, S_WORK_SLANSY,
237 + S_WORK_SPOT02, S_WORK_SPOT03 )
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*
579 END
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:101
spftrf
subroutine spftrf(transr, uplo, n, a, info)
SPFTRF
Definition spftrf.f:196
spftri
subroutine spftri(transr, uplo, n, a, info)
SPFTRI
Definition spftri.f:189
spftrs
subroutine spftrs(transr, uplo, n, nrhs, a, b, ldb, info)
SPFTRS
Definition spftrs.f:197
spotrf
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:105
spotri
subroutine spotri(uplo, n, a, lda, info)
SPOTRI
Definition spotri.f:93
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:194
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:192
sdrvrfp
subroutine sdrvrfp(nout, nn, nval, nns, nsval, nnt, ntval, thresh, a, asav, afac, ainv, b, bsav, xact, x, arf, arfinv, s_work_slatms, s_work_spot01, s_temp_spot02, s_temp_spot03, s_work_slansy, s_work_spot02, s_work_spot03)
SDRVRFP
Definition sdrvrfp.f:238
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