LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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,
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
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:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:73
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine alahd(IOUNIT, PATH)
ALAHD
Definition: alahd.f:107
subroutine alaerh(PATH, SUBNAM, INFO, INFOE, OPTS, M, N, KL, KU, N5, IMAT, NFAIL, NERRS, NOUT)
ALAERH
Definition: alaerh.f:147
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:321
subroutine spbrfs(UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
SPBRFS
Definition: spbrfs.f:189
subroutine spbtrs(UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO)
SPBTRS
Definition: spbtrs.f:121
subroutine spbcon(UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK, IWORK, INFO)
SPBCON
Definition: spbcon.f:132
subroutine spbtrf(UPLO, N, KD, AB, LDAB, INFO)
SPBTRF
Definition: spbtrf.f:142
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine slarhs(PATH, XTYPE, UPLO, TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B, LDB, ISEED, INFO)
SLARHS
Definition: slarhs.f:205
subroutine spbt02(UPLO, N, KD, NRHS, A, LDA, X, LDX, B, LDB, RWORK, RESID)
SPBT02
Definition: spbt02.f:136
subroutine spbt01(UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK, RESID)
SPBT01
Definition: spbt01.f:119
subroutine serrpo(PATH, NUNIT)
SERRPO
Definition: serrpo.f:55
subroutine slatb4(PATH, IMAT, M, N, TYPE, KL, KU, ANORM, MODE, CNDNUM, DIST)
SLATB4
Definition: slatb4.f:120
subroutine sget04(N, NRHS, X, LDX, XACT, LDXACT, RCOND, RESID)
SGET04
Definition: sget04.f:102
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
subroutine spbt05(UPLO, N, KD, NRHS, AB, LDAB, B, LDB, X, LDX, XACT, LDXACT, FERR, BERR, RESLTS)
SPBT05
Definition: spbt05.f:171