LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cdrvpt.f
Go to the documentation of this file.
1*> \brief \b CDRVPT
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 CDRVPT( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D,
12* E, B, X, XACT, WORK, RWORK, NOUT )
13*
14* .. Scalar Arguments ..
15* LOGICAL TSTERR
16* INTEGER NN, NOUT, NRHS
17* REAL THRESH
18* ..
19* .. Array Arguments ..
20* LOGICAL DOTYPE( * )
21* INTEGER NVAL( * )
22* REAL D( * ), RWORK( * )
23* COMPLEX A( * ), B( * ), E( * ), WORK( * ), X( * ),
24* $ XACT( * )
25* ..
26*
27*
28*> \par Purpose:
29* =============
30*>
31*> \verbatim
32*>
33*> CDRVPT tests CPTSV and -SVX.
34*> \endverbatim
35*
36* Arguments:
37* ==========
38*
39*> \param[in] DOTYPE
40*> \verbatim
41*> DOTYPE is LOGICAL array, dimension (NTYPES)
42*> The matrix types to be used for testing. Matrices of type j
43*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
44*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
45*> \endverbatim
46*>
47*> \param[in] NN
48*> \verbatim
49*> NN is INTEGER
50*> The number of values of N contained in the vector NVAL.
51*> \endverbatim
52*>
53*> \param[in] NVAL
54*> \verbatim
55*> NVAL is INTEGER array, dimension (NN)
56*> The values of the matrix dimension N.
57*> \endverbatim
58*>
59*> \param[in] NRHS
60*> \verbatim
61*> NRHS is INTEGER
62*> The number of right hand side vectors to be generated for
63*> each linear system.
64*> \endverbatim
65*>
66*> \param[in] THRESH
67*> \verbatim
68*> THRESH is REAL
69*> The threshold value for the test ratios. A result is
70*> included in the output file if RESULT >= THRESH. To have
71*> every test ratio printed, use THRESH = 0.
72*> \endverbatim
73*>
74*> \param[in] TSTERR
75*> \verbatim
76*> TSTERR is LOGICAL
77*> Flag that indicates whether error exits are to be tested.
78*> \endverbatim
79*>
80*> \param[out] A
81*> \verbatim
82*> A is COMPLEX array, dimension (NMAX*2)
83*> \endverbatim
84*>
85*> \param[out] D
86*> \verbatim
87*> D is REAL array, dimension (NMAX*2)
88*> \endverbatim
89*>
90*> \param[out] E
91*> \verbatim
92*> E is COMPLEX array, dimension (NMAX*2)
93*> \endverbatim
94*>
95*> \param[out] B
96*> \verbatim
97*> B is COMPLEX array, dimension (NMAX*NRHS)
98*> \endverbatim
99*>
100*> \param[out] X
101*> \verbatim
102*> X is COMPLEX array, dimension (NMAX*NRHS)
103*> \endverbatim
104*>
105*> \param[out] XACT
106*> \verbatim
107*> XACT is COMPLEX array, dimension (NMAX*NRHS)
108*> \endverbatim
109*>
110*> \param[out] WORK
111*> \verbatim
112*> WORK is COMPLEX array, dimension
113*> (NMAX*max(3,NRHS))
114*> \endverbatim
115*>
116*> \param[out] RWORK
117*> \verbatim
118*> RWORK is REAL array, dimension (NMAX+2*NRHS)
119*> \endverbatim
120*>
121*> \param[in] NOUT
122*> \verbatim
123*> NOUT is INTEGER
124*> The unit number for output.
125*> \endverbatim
126*
127* Authors:
128* ========
129*
130*> \author Univ. of Tennessee
131*> \author Univ. of California Berkeley
132*> \author Univ. of Colorado Denver
133*> \author NAG Ltd.
134*
135*> \ingroup complex_lin
136*
137* =====================================================================
138 SUBROUTINE cdrvpt( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, D,
139 $ E, B, X, XACT, WORK, RWORK, NOUT )
140*
141* -- LAPACK test routine --
142* -- LAPACK is a software package provided by Univ. of Tennessee, --
143* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144*
145* .. Scalar Arguments ..
146 LOGICAL TSTERR
147 INTEGER NN, NOUT, NRHS
148 REAL THRESH
149* ..
150* .. Array Arguments ..
151 LOGICAL DOTYPE( * )
152 INTEGER NVAL( * )
153 REAL D( * ), RWORK( * )
154 COMPLEX A( * ), B( * ), E( * ), WORK( * ), X( * ),
155 $ xact( * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameters ..
161 REAL ONE, ZERO
162 parameter( one = 1.0e+0, zero = 0.0e+0 )
163 INTEGER NTYPES
164 parameter( ntypes = 12 )
165 INTEGER NTESTS
166 parameter( ntests = 6 )
167* ..
168* .. Local Scalars ..
169 LOGICAL ZEROT
170 CHARACTER DIST, FACT, TYPE
171 CHARACTER*3 PATH
172 INTEGER I, IA, IFACT, IMAT, IN, INFO, IX, IZERO, J, K,
173 $ k1, kl, ku, lda, mode, n, nerrs, nfail, nimat,
174 $ nrun, nt
175 REAL AINVNM, ANORM, COND, DMAX, RCOND, RCONDC
176* ..
177* .. Local Arrays ..
178 INTEGER ISEED( 4 ), ISEEDY( 4 )
179 REAL RESULT( NTESTS ), Z( 3 )
180* ..
181* .. External Functions ..
182 INTEGER ISAMAX
183 REAL CLANHT, SCASUM, SGET06
184 EXTERNAL isamax, clanht, scasum, sget06
185* ..
186* .. External Subroutines ..
187 EXTERNAL aladhd, alaerh, alasvm, ccopy, cerrvx, cget04,
191* ..
192* .. Intrinsic Functions ..
193 INTRINSIC abs, cmplx, max
194* ..
195* .. Scalars in Common ..
196 LOGICAL LERR, OK
197 CHARACTER*32 SRNAMT
198 INTEGER INFOT, NUNIT
199* ..
200* .. Common blocks ..
201 COMMON / infoc / infot, nunit, ok, lerr
202 COMMON / srnamc / srnamt
203* ..
204* .. Data statements ..
205 DATA iseedy / 0, 0, 0, 1 /
206* ..
207* .. Executable Statements ..
208*
209 path( 1: 1 ) = 'Complex precision'
210 path( 2: 3 ) = 'PT'
211 nrun = 0
212 nfail = 0
213 nerrs = 0
214 DO 10 i = 1, 4
215 iseed( i ) = iseedy( i )
216 10 CONTINUE
217*
218* Test the error exits
219*
220 IF( tsterr )
221 $ CALL cerrvx( path, nout )
222 infot = 0
223*
224 DO 120 in = 1, nn
225*
226* Do for each value of N in NVAL.
227*
228 n = nval( in )
229 lda = max( 1, n )
230 nimat = ntypes
231 IF( n.LE.0 )
232 $ nimat = 1
233*
234 DO 110 imat = 1, nimat
235*
236* Do the tests only if DOTYPE( IMAT ) is true.
237*
238 IF( n.GT.0 .AND. .NOT.dotype( imat ) )
239 $ GO TO 110
240*
241* Set up parameters with CLATB4.
242*
243 CALL clatb4( path, imat, n, n, TYPE, kl, ku, anorm, mode,
244 $ cond, dist )
245*
246 zerot = imat.GE.8 .AND. imat.LE.10
247 IF( imat.LE.6 ) THEN
248*
249* Type 1-6: generate a symmetric tridiagonal matrix of
250* known condition number in lower triangular band storage.
251*
252 srnamt = 'CLATMS'
253 CALL clatms( n, n, dist, iseed, TYPE, rwork, mode, cond,
254 $ anorm, kl, ku, 'B', a, 2, work, info )
255*
256* Check the error code from CLATMS.
257*
258 IF( info.NE.0 ) THEN
259 CALL alaerh( path, 'CLATMS', info, 0, ' ', n, n, kl,
260 $ ku, -1, imat, nfail, nerrs, nout )
261 GO TO 110
262 END IF
263 izero = 0
264*
265* Copy the matrix to D and E.
266*
267 ia = 1
268 DO 20 i = 1, n - 1
269 d( i ) = real( a( ia ) )
270 e( i ) = a( ia+1 )
271 ia = ia + 2
272 20 CONTINUE
273 IF( n.GT.0 )
274 $ d( n ) = real( a( ia ) )
275 ELSE
276*
277* Type 7-12: generate a diagonally dominant matrix with
278* unknown condition number in the vectors D and E.
279*
280 IF( .NOT.zerot .OR. .NOT.dotype( 7 ) ) THEN
281*
282* Let D and E have values from [-1,1].
283*
284 CALL slarnv( 2, iseed, n, d )
285 CALL clarnv( 2, iseed, n-1, e )
286*
287* Make the tridiagonal matrix diagonally dominant.
288*
289 IF( n.EQ.1 ) THEN
290 d( 1 ) = abs( d( 1 ) )
291 ELSE
292 d( 1 ) = abs( d( 1 ) ) + abs( e( 1 ) )
293 d( n ) = abs( d( n ) ) + abs( e( n-1 ) )
294 DO 30 i = 2, n - 1
295 d( i ) = abs( d( i ) ) + abs( e( i ) ) +
296 $ abs( e( i-1 ) )
297 30 CONTINUE
298 END IF
299*
300* Scale D and E so the maximum element is ANORM.
301*
302 ix = isamax( n, d, 1 )
303 dmax = d( ix )
304 CALL sscal( n, anorm / dmax, d, 1 )
305 IF( n.GT.1 )
306 $ CALL csscal( n-1, anorm / dmax, e, 1 )
307*
308 ELSE IF( izero.GT.0 ) THEN
309*
310* Reuse the last matrix by copying back the zeroed out
311* elements.
312*
313 IF( izero.EQ.1 ) THEN
314 d( 1 ) = z( 2 )
315 IF( n.GT.1 )
316 $ e( 1 ) = z( 3 )
317 ELSE IF( izero.EQ.n ) THEN
318 e( n-1 ) = z( 1 )
319 d( n ) = z( 2 )
320 ELSE
321 e( izero-1 ) = z( 1 )
322 d( izero ) = z( 2 )
323 e( izero ) = z( 3 )
324 END IF
325 END IF
326*
327* For types 8-10, set one row and column of the matrix to
328* zero.
329*
330 izero = 0
331 IF( imat.EQ.8 ) THEN
332 izero = 1
333 z( 2 ) = d( 1 )
334 d( 1 ) = zero
335 IF( n.GT.1 ) THEN
336 z( 3 ) = real( e( 1 ) )
337 e( 1 ) = zero
338 END IF
339 ELSE IF( imat.EQ.9 ) THEN
340 izero = n
341 IF( n.GT.1 ) THEN
342 z( 1 ) = real( e( n-1 ) )
343 e( n-1 ) = zero
344 END IF
345 z( 2 ) = d( n )
346 d( n ) = zero
347 ELSE IF( imat.EQ.10 ) THEN
348 izero = ( n+1 ) / 2
349 IF( izero.GT.1 ) THEN
350 z( 1 ) = real( e( izero-1 ) )
351 e( izero-1 ) = zero
352 z( 3 ) = real( e( izero ) )
353 e( izero ) = zero
354 END IF
355 z( 2 ) = d( izero )
356 d( izero ) = zero
357 END IF
358 END IF
359*
360* Generate NRHS random solution vectors.
361*
362 ix = 1
363 DO 40 j = 1, nrhs
364 CALL clarnv( 2, iseed, n, xact( ix ) )
365 ix = ix + lda
366 40 CONTINUE
367*
368* Set the right hand side.
369*
370 CALL claptm( 'Lower', n, nrhs, one, d, e, xact, lda, zero,
371 $ b, lda )
372*
373 DO 100 ifact = 1, 2
374 IF( ifact.EQ.1 ) THEN
375 fact = 'F'
376 ELSE
377 fact = 'N'
378 END IF
379*
380* Compute the condition number for comparison with
381* the value returned by CPTSVX.
382*
383 IF( zerot ) THEN
384 IF( ifact.EQ.1 )
385 $ GO TO 100
386 rcondc = zero
387*
388 ELSE IF( ifact.EQ.1 ) THEN
389*
390* Compute the 1-norm of A.
391*
392 anorm = clanht( '1', n, d, e )
393*
394 CALL scopy( n, d, 1, d( n+1 ), 1 )
395 IF( n.GT.1 )
396 $ CALL ccopy( n-1, e, 1, e( n+1 ), 1 )
397*
398* Factor the matrix A.
399*
400 CALL cpttrf( n, d( n+1 ), e( n+1 ), info )
401*
402* Use CPTTRS to solve for one column at a time of
403* inv(A), computing the maximum column sum as we go.
404*
405 ainvnm = zero
406 DO 60 i = 1, n
407 DO 50 j = 1, n
408 x( j ) = zero
409 50 CONTINUE
410 x( i ) = one
411 CALL cpttrs( 'Lower', n, 1, d( n+1 ), e( n+1 ), x,
412 $ lda, info )
413 ainvnm = max( ainvnm, scasum( n, x, 1 ) )
414 60 CONTINUE
415*
416* Compute the 1-norm condition number of A.
417*
418 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
419 rcondc = one
420 ELSE
421 rcondc = ( one / anorm ) / ainvnm
422 END IF
423 END IF
424*
425 IF( ifact.EQ.2 ) THEN
426*
427* --- Test CPTSV --
428*
429 CALL scopy( n, d, 1, d( n+1 ), 1 )
430 IF( n.GT.1 )
431 $ CALL ccopy( n-1, e, 1, e( n+1 ), 1 )
432 CALL clacpy( 'Full', n, nrhs, b, lda, x, lda )
433*
434* Factor A as L*D*L' and solve the system A*X = B.
435*
436 srnamt = 'CPTSV '
437 CALL cptsv( n, nrhs, d( n+1 ), e( n+1 ), x, lda,
438 $ info )
439*
440* Check error code from CPTSV .
441*
442 IF( info.NE.izero )
443 $ CALL alaerh( path, 'CPTSV ', info, izero, ' ', n,
444 $ n, 1, 1, nrhs, imat, nfail, nerrs,
445 $ nout )
446 nt = 0
447 IF( izero.EQ.0 ) THEN
448*
449* Check the factorization by computing the ratio
450* norm(L*D*L' - A) / (n * norm(A) * EPS )
451*
452 CALL cptt01( n, d, e, d( n+1 ), e( n+1 ), work,
453 $ result( 1 ) )
454*
455* Compute the residual in the solution.
456*
457 CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
458 CALL cptt02( 'Lower', n, nrhs, d, e, x, lda, work,
459 $ lda, result( 2 ) )
460*
461* Check solution from generated exact solution.
462*
463 CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
464 $ result( 3 ) )
465 nt = 3
466 END IF
467*
468* Print information about the tests that did not pass
469* the threshold.
470*
471 DO 70 k = 1, nt
472 IF( result( k ).GE.thresh ) THEN
473 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
474 $ CALL aladhd( nout, path )
475 WRITE( nout, fmt = 9999 )'CPTSV ', n, imat, k,
476 $ result( k )
477 nfail = nfail + 1
478 END IF
479 70 CONTINUE
480 nrun = nrun + nt
481 END IF
482*
483* --- Test CPTSVX ---
484*
485 IF( ifact.GT.1 ) THEN
486*
487* Initialize D( N+1:2*N ) and E( N+1:2*N ) to zero.
488*
489 DO 80 i = 1, n - 1
490 d( n+i ) = zero
491 e( n+i ) = zero
492 80 CONTINUE
493 IF( n.GT.0 )
494 $ d( n+n ) = zero
495 END IF
496*
497 CALL claset( 'Full', n, nrhs, cmplx( zero ),
498 $ cmplx( zero ), x, lda )
499*
500* Solve the system and compute the condition number and
501* error bounds using CPTSVX.
502*
503 srnamt = 'CPTSVX'
504 CALL cptsvx( fact, n, nrhs, d, e, d( n+1 ), e( n+1 ), b,
505 $ lda, x, lda, rcond, rwork, rwork( nrhs+1 ),
506 $ work, rwork( 2*nrhs+1 ), info )
507*
508* Check the error code from CPTSVX.
509*
510 IF( info.NE.izero )
511 $ CALL alaerh( path, 'CPTSVX', info, izero, fact, n, n,
512 $ 1, 1, nrhs, imat, nfail, nerrs, nout )
513 IF( izero.EQ.0 ) THEN
514 IF( ifact.EQ.2 ) THEN
515*
516* Check the factorization by computing the ratio
517* norm(L*D*L' - A) / (n * norm(A) * EPS )
518*
519 k1 = 1
520 CALL cptt01( n, d, e, d( n+1 ), e( n+1 ), work,
521 $ result( 1 ) )
522 ELSE
523 k1 = 2
524 END IF
525*
526* Compute the residual in the solution.
527*
528 CALL clacpy( 'Full', n, nrhs, b, lda, work, lda )
529 CALL cptt02( 'Lower', n, nrhs, d, e, x, lda, work,
530 $ lda, result( 2 ) )
531*
532* Check solution from generated exact solution.
533*
534 CALL cget04( n, nrhs, x, lda, xact, lda, rcondc,
535 $ result( 3 ) )
536*
537* Check error bounds from iterative refinement.
538*
539 CALL cptt05( n, nrhs, d, e, b, lda, x, lda, xact, lda,
540 $ rwork, rwork( nrhs+1 ), result( 4 ) )
541 ELSE
542 k1 = 6
543 END IF
544*
545* Check the reciprocal of the condition number.
546*
547 result( 6 ) = sget06( rcond, rcondc )
548*
549* Print information about the tests that did not pass
550* the threshold.
551*
552 DO 90 k = k1, 6
553 IF( result( k ).GE.thresh ) THEN
554 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
555 $ CALL aladhd( nout, path )
556 WRITE( nout, fmt = 9998 )'CPTSVX', fact, n, imat,
557 $ k, result( k )
558 nfail = nfail + 1
559 END IF
560 90 CONTINUE
561 nrun = nrun + 7 - k1
562 100 CONTINUE
563 110 CONTINUE
564 120 CONTINUE
565*
566* Print a summary of the results.
567*
568 CALL alasvm( path, nout, nfail, nrun, nerrs )
569*
570 9999 FORMAT( 1x, a, ', N =', i5, ', type ', i2, ', test ', i2,
571 $ ', ratio = ', g12.5 )
572 9998 FORMAT( 1x, a, ', FACT=''', a1, ''', N =', i5, ', type ', i2,
573 $ ', test ', i2, ', ratio = ', g12.5 )
574 RETURN
575*
576* End of CDRVPT
577*
578 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
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 cdrvpt(dotype, nn, nval, nrhs, thresh, tsterr, a, d, e, b, x, xact, work, rwork, nout)
CDRVPT
Definition cdrvpt.f:140
subroutine cerrvx(path, nunit)
CERRVX
Definition cerrvx.f:55
subroutine cget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
CGET04
Definition cget04.f:102
subroutine claptm(uplo, n, nrhs, alpha, d, e, x, ldx, beta, b, ldb)
CLAPTM
Definition claptm.f:129
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 cptt01(n, d, e, df, ef, work, resid)
CPTT01
Definition cptt01.f:92
subroutine cptt02(uplo, n, nrhs, d, e, x, ldx, b, ldb, resid)
CPTT02
Definition cptt02.f:115
subroutine cptt05(n, nrhs, d, e, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
CPTT05
Definition cptt05.f:150
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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 slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:99
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
subroutine cptsv(n, nrhs, d, e, b, ldb, info)
CPTSV computes the solution to system of linear equations A * X = B for PT matrices
Definition cptsv.f:115
subroutine cptsvx(fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
CPTSVX computes the solution to system of linear equations A * X = B for PT matrices
Definition cptsvx.f:234
subroutine cpttrf(n, d, e, info)
CPTTRF
Definition cpttrf.f:92
subroutine cpttrs(uplo, n, nrhs, d, e, b, ldb, info)
CPTTRS
Definition cpttrs.f:121
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79