LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
schkpb.f
Go to the documentation of this file.
1*> \brief \b SCHKPB
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 SCHKPB( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
12* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
13* XACT, WORK, RWORK, IWORK, NOUT )
14*
15* .. Scalar Arguments ..
16* LOGICAL TSTERR
17* INTEGER NMAX, NN, NNB, NNS, NOUT
18* REAL THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL DOTYPE( * )
22* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
23* REAL A( * ), AFAC( * ), AINV( * ), B( * ),
24* $ RWORK( * ), WORK( * ), X( * ), XACT( * )
25* ..
26*
27*
28*> \par Purpose:
29* =============
30*>
31*> \verbatim
32*>
33*> SCHKPB tests SPBTRF, -TRS, -RFS, and -CON.
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] NNB
60*> \verbatim
61*> NNB is INTEGER
62*> The number of values of NB contained in the vector NBVAL.
63*> \endverbatim
64*>
65*> \param[in] NBVAL
66*> \verbatim
67*> NBVAL is INTEGER array, dimension (NNB)
68*> The values of the blocksize NB.
69*> \endverbatim
70*>
71*> \param[in] NNS
72*> \verbatim
73*> NNS is INTEGER
74*> The number of values of NRHS contained in the vector NSVAL.
75*> \endverbatim
76*>
77*> \param[in] NSVAL
78*> \verbatim
79*> NSVAL is INTEGER array, dimension (NNS)
80*> The values of the number of right hand sides NRHS.
81*> \endverbatim
82*>
83*> \param[in] THRESH
84*> \verbatim
85*> THRESH is REAL
86*> The threshold value for the test ratios. A result is
87*> included in the output file if RESULT >= THRESH. To have
88*> every test ratio printed, use THRESH = 0.
89*> \endverbatim
90*>
91*> \param[in] TSTERR
92*> \verbatim
93*> TSTERR is LOGICAL
94*> Flag that indicates whether error exits are to be tested.
95*> \endverbatim
96*>
97*> \param[in] NMAX
98*> \verbatim
99*> NMAX is INTEGER
100*> The maximum value permitted for N, used in dimensioning the
101*> work arrays.
102*> \endverbatim
103*>
104*> \param[out] A
105*> \verbatim
106*> A is REAL array, dimension (NMAX*NMAX)
107*> \endverbatim
108*>
109*> \param[out] AFAC
110*> \verbatim
111*> AFAC is REAL array, dimension (NMAX*NMAX)
112*> \endverbatim
113*>
114*> \param[out] AINV
115*> \verbatim
116*> AINV is REAL array, dimension (NMAX*NMAX)
117*> \endverbatim
118*>
119*> \param[out] B
120*> \verbatim
121*> B is REAL array, dimension (NMAX*NSMAX)
122*> where NSMAX is the largest entry in NSVAL.
123*> \endverbatim
124*>
125*> \param[out] X
126*> \verbatim
127*> X is REAL array, dimension (NMAX*NSMAX)
128*> \endverbatim
129*>
130*> \param[out] XACT
131*> \verbatim
132*> XACT is REAL array, dimension (NMAX*NSMAX)
133*> \endverbatim
134*>
135*> \param[out] WORK
136*> \verbatim
137*> WORK is REAL array, dimension
138*> (NMAX*max(3,NSMAX))
139*> \endverbatim
140*>
141*> \param[out] RWORK
142*> \verbatim
143*> RWORK is REAL array, dimension
144*> (max(NMAX,2*NSMAX))
145*> \endverbatim
146*>
147*> \param[out] IWORK
148*> \verbatim
149*> IWORK is INTEGER array, dimension (NMAX)
150*> \endverbatim
151*>
152*> \param[in] NOUT
153*> \verbatim
154*> NOUT is INTEGER
155*> The unit number for output.
156*> \endverbatim
157*
158* Authors:
159* ========
160*
161*> \author Univ. of Tennessee
162*> \author Univ. of California Berkeley
163*> \author Univ. of Colorado Denver
164*> \author NAG Ltd.
165*
166*> \ingroup single_lin
167*
168* =====================================================================
169 SUBROUTINE schkpb( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
170 $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
171 $ XACT, WORK, RWORK, IWORK, NOUT )
172*
173* -- LAPACK test routine --
174* -- LAPACK is a software package provided by Univ. of Tennessee, --
175* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176*
177* .. Scalar Arguments ..
178 LOGICAL TSTERR
179 INTEGER NMAX, NN, NNB, NNS, NOUT
180 REAL THRESH
181* ..
182* .. Array Arguments ..
183 LOGICAL DOTYPE( * )
184 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
185 REAL A( * ), AFAC( * ), AINV( * ), B( * ),
186 $ rwork( * ), work( * ), x( * ), xact( * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 REAL ONE, ZERO
193 PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
194 INTEGER NTYPES, NTESTS
195 parameter( ntypes = 8, ntests = 7 )
196 INTEGER NBW
197 parameter( nbw = 4 )
198* ..
199* .. Local Scalars ..
200 LOGICAL ZEROT
201 CHARACTER DIST, PACKIT, TYPE, UPLO, XTYPE
202 CHARACTER*3 PATH
203 INTEGER I, I1, I2, IKD, IMAT, IN, INB, INFO, IOFF,
204 $ irhs, iuplo, iw, izero, k, kd, kl, koff, ku,
205 $ lda, ldab, mode, n, nb, nerrs, nfail, nimat,
206 $ nkd, nrhs, nrun
207 REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC
208* ..
209* .. Local Arrays ..
210 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
211 REAL RESULT( NTESTS )
212* ..
213* .. External Functions ..
214 REAL SGET06, SLANGE, SLANSB
215 EXTERNAL SGET06, SLANGE, SLANSB
216* ..
217* .. External Subroutines ..
218 EXTERNAL alaerh, alahd, alasum, scopy, serrpo, sget04,
219 $ slacpy, slarhs, slaset, slatb4, slatms, spbcon,
220 $ spbrfs, spbt01, spbt02, spbt05, spbtrf, spbtrs,
221 $ sswap, xlaenv
222* ..
223* .. Intrinsic Functions ..
224 INTRINSIC max, min
225* ..
226* .. Scalars in Common ..
227 LOGICAL LERR, OK
228 CHARACTER*32 SRNAMT
229 INTEGER INFOT, NUNIT
230* ..
231* .. Common blocks ..
232 COMMON / infoc / infot, nunit, ok, lerr
233 COMMON / srnamc / srnamt
234* ..
235* .. Data statements ..
236 DATA iseedy / 1988, 1989, 1990, 1991 /
237* ..
238* .. Executable Statements ..
239*
240* Initialize constants and the random number seed.
241*
242 path( 1: 1 ) = 'Single precision'
243 path( 2: 3 ) = 'PB'
244 nrun = 0
245 nfail = 0
246 nerrs = 0
247 DO 10 i = 1, 4
248 iseed( i ) = iseedy( i )
249 10 CONTINUE
250*
251* Test the error exits
252*
253 IF( tsterr )
254 $ CALL serrpo( path, nout )
255 infot = 0
256 CALL xlaenv( 2, 2 )
257 kdval( 1 ) = 0
258*
259* Do for each value of N in NVAL
260*
261 DO 90 in = 1, nn
262 n = nval( in )
263 lda = max( n, 1 )
264 xtype = 'N'
265*
266* Set limits on the number of loop iterations.
267*
268 nkd = max( 1, min( n, 4 ) )
269 nimat = ntypes
270 IF( n.EQ.0 )
271 $ nimat = 1
272*
273 kdval( 2 ) = n + ( n+1 ) / 4
274 kdval( 3 ) = ( 3*n-1 ) / 4
275 kdval( 4 ) = ( n+1 ) / 4
276*
277 DO 80 ikd = 1, nkd
278*
279* Do for KD = 0, (5*N+1)/4, (3N-1)/4, and (N+1)/4. This order
280* makes it easier to skip redundant values for small values
281* of N.
282*
283 kd = kdval( ikd )
284 ldab = kd + 1
285*
286* Do first for UPLO = 'U', then for UPLO = 'L'
287*
288 DO 70 iuplo = 1, 2
289 koff = 1
290 IF( iuplo.EQ.1 ) THEN
291 uplo = 'U'
292 koff = max( 1, kd+2-n )
293 packit = 'Q'
294 ELSE
295 uplo = 'L'
296 packit = 'B'
297 END IF
298*
299 DO 60 imat = 1, nimat
300*
301* Do the tests only if DOTYPE( IMAT ) is true.
302*
303 IF( .NOT.dotype( imat ) )
304 $ GO TO 60
305*
306* Skip types 2, 3, or 4 if the matrix size is too small.
307*
308 zerot = imat.GE.2 .AND. imat.LE.4
309 IF( zerot .AND. n.LT.imat-1 )
310 $ GO TO 60
311*
312 IF( .NOT.zerot .OR. .NOT.dotype( 1 ) ) THEN
313*
314* Set up parameters with SLATB4 and generate a test
315* matrix with SLATMS.
316*
317 CALL slatb4( path, imat, n, n, TYPE, kl, ku, anorm,
318 $ mode, cndnum, dist )
319*
320 srnamt = 'SLATMS'
321 CALL slatms( n, n, dist, iseed, TYPE, rwork, mode,
322 $ cndnum, anorm, kd, kd, packit,
323 $ a( koff ), ldab, work, info )
324*
325* Check error code from SLATMS.
326*
327 IF( info.NE.0 ) THEN
328 CALL alaerh( path, 'SLATMS', info, 0, uplo, n,
329 $ n, kd, kd, -1, imat, nfail, nerrs,
330 $ nout )
331 GO TO 60
332 END IF
333 ELSE IF( izero.GT.0 ) THEN
334*
335* Use the same matrix for types 3 and 4 as for type
336* 2 by copying back the zeroed out column,
337*
338 iw = 2*lda + 1
339 IF( iuplo.EQ.1 ) THEN
340 ioff = ( izero-1 )*ldab + kd + 1
341 CALL scopy( izero-i1, work( iw ), 1,
342 $ a( ioff-izero+i1 ), 1 )
343 iw = iw + izero - i1
344 CALL scopy( i2-izero+1, work( iw ), 1,
345 $ a( ioff ), max( ldab-1, 1 ) )
346 ELSE
347 ioff = ( i1-1 )*ldab + 1
348 CALL scopy( izero-i1, work( iw ), 1,
349 $ a( ioff+izero-i1 ),
350 $ max( ldab-1, 1 ) )
351 ioff = ( izero-1 )*ldab + 1
352 iw = iw + izero - i1
353 CALL scopy( i2-izero+1, work( iw ), 1,
354 $ a( ioff ), 1 )
355 END IF
356 END IF
357*
358* For types 2-4, zero one row and column of the matrix
359* to test that INFO is returned correctly.
360*
361 izero = 0
362 IF( zerot ) THEN
363 IF( imat.EQ.2 ) THEN
364 izero = 1
365 ELSE IF( imat.EQ.3 ) THEN
366 izero = n
367 ELSE
368 izero = n / 2 + 1
369 END IF
370*
371* Save the zeroed out row and column in WORK(*,3)
372*
373 iw = 2*lda
374 DO 20 i = 1, min( 2*kd+1, n )
375 work( iw+i ) = zero
376 20 CONTINUE
377 iw = iw + 1
378 i1 = max( izero-kd, 1 )
379 i2 = min( izero+kd, n )
380*
381 IF( iuplo.EQ.1 ) THEN
382 ioff = ( izero-1 )*ldab + kd + 1
383 CALL sswap( izero-i1, a( ioff-izero+i1 ), 1,
384 $ work( iw ), 1 )
385 iw = iw + izero - i1
386 CALL sswap( i2-izero+1, a( ioff ),
387 $ max( ldab-1, 1 ), work( iw ), 1 )
388 ELSE
389 ioff = ( i1-1 )*ldab + 1
390 CALL sswap( izero-i1, a( ioff+izero-i1 ),
391 $ max( ldab-1, 1 ), work( iw ), 1 )
392 ioff = ( izero-1 )*ldab + 1
393 iw = iw + izero - i1
394 CALL sswap( i2-izero+1, a( ioff ), 1,
395 $ work( iw ), 1 )
396 END IF
397 END IF
398*
399* Do for each value of NB in NBVAL
400*
401 DO 50 inb = 1, nnb
402 nb = nbval( inb )
403 CALL xlaenv( 1, nb )
404*
405* Compute the L*L' or U'*U factorization of the band
406* matrix.
407*
408 CALL slacpy( 'Full', kd+1, n, a, ldab, afac, ldab )
409 srnamt = 'SPBTRF'
410 CALL spbtrf( uplo, n, kd, afac, ldab, info )
411*
412* Check error code from SPBTRF.
413*
414 IF( info.NE.izero ) THEN
415 CALL alaerh( path, 'SPBTRF', info, izero, uplo,
416 $ n, n, kd, kd, nb, imat, nfail,
417 $ nerrs, nout )
418 GO TO 50
419 END IF
420*
421* Skip the tests if INFO is not 0.
422*
423 IF( info.NE.0 )
424 $ GO TO 50
425*
426*+ TEST 1
427* Reconstruct matrix from factors and compute
428* residual.
429*
430 CALL slacpy( 'Full', kd+1, n, afac, ldab, ainv,
431 $ ldab )
432 CALL spbt01( uplo, n, kd, a, ldab, ainv, ldab,
433 $ rwork, result( 1 ) )
434*
435* Print the test ratio if it is .GE. THRESH.
436*
437 IF( result( 1 ).GE.thresh ) THEN
438 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
439 $ CALL alahd( nout, path )
440 WRITE( nout, fmt = 9999 )uplo, n, kd, nb, imat,
441 $ 1, result( 1 )
442 nfail = nfail + 1
443 END IF
444 nrun = nrun + 1
445*
446* Only do other tests if this is the first blocksize.
447*
448 IF( inb.GT.1 )
449 $ GO TO 50
450*
451* Form the inverse of A so we can get a good estimate
452* of RCONDC = 1/(norm(A) * norm(inv(A))).
453*
454 CALL slaset( 'Full', n, n, zero, one, ainv, lda )
455 srnamt = 'SPBTRS'
456 CALL spbtrs( uplo, n, kd, n, afac, ldab, ainv, lda,
457 $ info )
458*
459* Compute RCONDC = 1/(norm(A) * norm(inv(A))).
460*
461 anorm = slansb( '1', uplo, n, kd, a, ldab, rwork )
462 ainvnm = slange( '1', n, n, ainv, lda, rwork )
463 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
464 rcondc = one
465 ELSE
466 rcondc = ( one / anorm ) / ainvnm
467 END IF
468*
469 DO 40 irhs = 1, nns
470 nrhs = nsval( irhs )
471*
472*+ TEST 2
473* Solve and compute residual for A * X = B.
474*
475 srnamt = 'SLARHS'
476 CALL slarhs( path, xtype, uplo, ' ', n, n, kd,
477 $ kd, nrhs, a, ldab, xact, lda, b,
478 $ lda, iseed, info )
479 CALL slacpy( 'Full', n, nrhs, b, lda, x, lda )
480*
481 srnamt = 'SPBTRS'
482 CALL spbtrs( uplo, n, kd, nrhs, afac, ldab, x,
483 $ lda, info )
484*
485* Check error code from SPBTRS.
486*
487 IF( info.NE.0 )
488 $ CALL alaerh( path, 'SPBTRS', info, 0, uplo,
489 $ n, n, kd, kd, nrhs, imat, nfail,
490 $ nerrs, nout )
491*
492 CALL slacpy( 'Full', n, nrhs, b, lda, work,
493 $ lda )
494 CALL spbt02( uplo, n, kd, nrhs, a, ldab, x, lda,
495 $ work, lda, rwork, result( 2 ) )
496*
497*+ TEST 3
498* Check solution from generated exact solution.
499*
500 CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
501 $ result( 3 ) )
502*
503*+ TESTS 4, 5, and 6
504* Use iterative refinement to improve the solution.
505*
506 srnamt = 'SPBRFS'
507 CALL spbrfs( uplo, n, kd, nrhs, a, ldab, afac,
508 $ ldab, b, lda, x, lda, rwork,
509 $ rwork( nrhs+1 ), work, iwork,
510 $ info )
511*
512* Check error code from SPBRFS.
513*
514 IF( info.NE.0 )
515 $ CALL alaerh( path, 'SPBRFS', info, 0, uplo,
516 $ n, n, kd, kd, nrhs, imat, nfail,
517 $ nerrs, nout )
518*
519 CALL sget04( n, nrhs, x, lda, xact, lda, rcondc,
520 $ result( 4 ) )
521 CALL spbt05( uplo, n, kd, nrhs, a, ldab, b, lda,
522 $ x, lda, xact, lda, rwork,
523 $ rwork( nrhs+1 ), result( 5 ) )
524*
525* Print information about the tests that did not
526* pass the threshold.
527*
528 DO 30 k = 2, 6
529 IF( result( k ).GE.thresh ) THEN
530 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
531 $ CALL alahd( nout, path )
532 WRITE( nout, fmt = 9998 )uplo, n, kd,
533 $ nrhs, imat, k, result( k )
534 nfail = nfail + 1
535 END IF
536 30 CONTINUE
537 nrun = nrun + 5
538 40 CONTINUE
539*
540*+ TEST 7
541* Get an estimate of RCOND = 1/CNDNUM.
542*
543 srnamt = 'SPBCON'
544 CALL spbcon( uplo, n, kd, afac, ldab, anorm, rcond,
545 $ work, iwork, info )
546*
547* Check error code from SPBCON.
548*
549 IF( info.NE.0 )
550 $ CALL alaerh( path, 'SPBCON', info, 0, uplo, n,
551 $ n, kd, kd, -1, imat, nfail, nerrs,
552 $ nout )
553*
554 result( 7 ) = sget06( rcond, rcondc )
555*
556* Print the test ratio if it is .GE. THRESH.
557*
558 IF( result( 7 ).GE.thresh ) THEN
559 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
560 $ CALL alahd( nout, path )
561 WRITE( nout, fmt = 9997 )uplo, n, kd, imat, 7,
562 $ result( 7 )
563 nfail = nfail + 1
564 END IF
565 nrun = nrun + 1
566 50 CONTINUE
567 60 CONTINUE
568 70 CONTINUE
569 80 CONTINUE
570 90 CONTINUE
571*
572* Print a summary of the results.
573*
574 CALL alasum( path, nout, nfail, nrun, nerrs )
575*
576 9999 FORMAT( ' UPLO=''', a1, ''', N=', i5, ', KD=', i5, ', NB=', i4,
577 $ ', type ', i2, ', test ', i2, ', ratio= ', g12.5 )
578 9998 FORMAT( ' UPLO=''', a1, ''', N=', i5, ', KD=', i5, ', NRHS=', i3,
579 $ ', type ', i2, ', test(', i2, ') = ', g12.5 )
580 9997 FORMAT( ' UPLO=''', a1, ''', N=', i5, ', KD=', i5, ',', 10x,
581 $ ' type ', i2, ', test(', i2, ') = ', g12.5 )
582 RETURN
583*
584* End of SCHKPB
585*
586 END
alasum
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
Definition alasum.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
xlaenv
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
alahd
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
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
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
spbcon
subroutine spbcon(uplo, n, kd, ab, ldab, anorm, rcond, work, iwork, info)
SPBCON
Definition spbcon.f:130
spbrfs
subroutine spbrfs(uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SPBRFS
Definition spbrfs.f:187
spbtrf
subroutine spbtrf(uplo, n, kd, ab, ldab, info)
SPBTRF
Definition spbtrf.f:140
spbtrs
subroutine spbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
SPBTRS
Definition spbtrs.f:119
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
schkpb
subroutine schkpb(dotype, nn, nval, nnb, nbval, nns, nsval, thresh, tsterr, nmax, a, afac, ainv, b, x, xact, work, rwork, iwork, nout)
SCHKPB
Definition schkpb.f:172
serrpo
subroutine serrpo(path, nunit)
SERRPO
Definition serrpo.f:55
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
spbt01
subroutine spbt01(uplo, n, kd, a, lda, afac, ldafac, rwork, resid)
SPBT01
Definition spbt01.f:119
spbt02
subroutine spbt02(uplo, n, kd, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
SPBT02
Definition spbt02.f:136
spbt05
subroutine spbt05(uplo, n, kd, nrhs, ab, ldab, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
SPBT05
Definition spbt05.f:171