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