LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dchkpb()

subroutine dchkpb ( logical, dimension( * )  dotype,
integer  nn,
integer, dimension( * )  nval,
integer  nnb,
integer, dimension( * )  nbval,
integer  nns,
integer, dimension( * )  nsval,
double precision  thresh,
logical  tsterr,
integer  nmax,
double precision, dimension( * )  a,
double precision, dimension( * )  afac,
double precision, dimension( * )  ainv,
double precision, dimension( * )  b,
double precision, dimension( * )  x,
double precision, dimension( * )  xact,
double precision, dimension( * )  work,
double precision, dimension( * )  rwork,
integer, dimension( * )  iwork,
integer  nout 
)

DCHKPB

Purpose:
 DCHKPB tests DPBTRF, -TRS, -RFS, and -CON.
Parameters
[in]DOTYPE
          DOTYPE is LOGICAL array, dimension (NTYPES)
          The matrix types to be used for testing.  Matrices of type j
          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
[in]NN
          NN is INTEGER
          The number of values of N contained in the vector NVAL.
[in]NVAL
          NVAL is INTEGER array, dimension (NN)
          The values of the matrix dimension N.
[in]NNB
          NNB is INTEGER
          The number of values of NB contained in the vector NBVAL.
[in]NBVAL
          NBVAL is INTEGER array, dimension (NNB)
          The values of the blocksize NB.
[in]NNS
          NNS is INTEGER
          The number of values of NRHS contained in the vector NSVAL.
[in]NSVAL
          NSVAL is INTEGER array, dimension (NNS)
          The values of the number of right hand sides NRHS.
[in]THRESH
          THRESH is DOUBLE PRECISION
          The threshold value for the test ratios.  A result is
          included in the output file if RESULT >= THRESH.  To have
          every test ratio printed, use THRESH = 0.
[in]TSTERR
          TSTERR is LOGICAL
          Flag that indicates whether error exits are to be tested.
[in]NMAX
          NMAX is INTEGER
          The maximum value permitted for N, used in dimensioning the
          work arrays.
[out]A
          A is DOUBLE PRECISION array, dimension (NMAX*NMAX)
[out]AFAC
          AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX)
[out]AINV
          AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
[out]B
          B is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
          where NSMAX is the largest entry in NSVAL.
[out]X
          X is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
[out]XACT
          XACT is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
[out]WORK
          WORK is DOUBLE PRECISION array, dimension
                      (NMAX*max(3,NSMAX))
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension
                      (max(NMAX,2*NSMAX))
[out]IWORK
          IWORK is INTEGER array, dimension (NMAX)
[in]NOUT
          NOUT is INTEGER
          The unit number for output.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 169 of file dchkpb.f.

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 DOUBLE PRECISION THRESH
181* ..
182* .. Array Arguments ..
183 LOGICAL DOTYPE( * )
184 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
185 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ),
186 $ RWORK( * ), WORK( * ), X( * ), XACT( * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 DOUBLE PRECISION ONE, ZERO
193 parameter( one = 1.0d+0, zero = 0.0d+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 DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC
208* ..
209* .. Local Arrays ..
210 INTEGER ISEED( 4 ), ISEEDY( 4 ), KDVAL( NBW )
211 DOUBLE PRECISION RESULT( NTESTS )
212* ..
213* .. External Functions ..
214 DOUBLE PRECISION DGET06, DLANGE, DLANSB
215 EXTERNAL dget06, dlange, dlansb
216* ..
217* .. External Subroutines ..
218 EXTERNAL alaerh, alahd, alasum, dcopy, derrpo, dget04,
219 $ dlacpy, dlarhs, dlaset, dlatb4, dlatms, dpbcon,
220 $ dpbrfs, dpbt01, dpbt02, dpbt05, dpbtrf, dpbtrs,
221 $ dswap, 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 ) = 'Double 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 derrpo( 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 DLATB4 and generate a test
315* matrix with DLATMS.
316*
317 CALL dlatb4( path, imat, n, n, TYPE, KL, KU, ANORM,
318 $ MODE, CNDNUM, DIST )
319*
320 srnamt = 'DLATMS'
321 CALL dlatms( 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 DLATMS.
326*
327 IF( info.NE.0 ) THEN
328 CALL alaerh( path, 'DLATMS', 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 dcopy( izero-i1, work( iw ), 1,
342 $ a( ioff-izero+i1 ), 1 )
343 iw = iw + izero - i1
344 CALL dcopy( i2-izero+1, work( iw ), 1,
345 $ a( ioff ), max( ldab-1, 1 ) )
346 ELSE
347 ioff = ( i1-1 )*ldab + 1
348 CALL dcopy( 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 dcopy( 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 dswap( izero-i1, a( ioff-izero+i1 ), 1,
384 $ work( iw ), 1 )
385 iw = iw + izero - i1
386 CALL dswap( i2-izero+1, a( ioff ),
387 $ max( ldab-1, 1 ), work( iw ), 1 )
388 ELSE
389 ioff = ( i1-1 )*ldab + 1
390 CALL dswap( 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 dswap( 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 dlacpy( 'Full', kd+1, n, a, ldab, afac, ldab )
409 srnamt = 'DPBTRF'
410 CALL dpbtrf( uplo, n, kd, afac, ldab, info )
411*
412* Check error code from DPBTRF.
413*
414 IF( info.NE.izero ) THEN
415 CALL alaerh( path, 'DPBTRF', 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 dlacpy( 'Full', kd+1, n, afac, ldab, ainv,
431 $ ldab )
432 CALL dpbt01( 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 dlaset( 'Full', n, n, zero, one, ainv, lda )
455 srnamt = 'DPBTRS'
456 CALL dpbtrs( uplo, n, kd, n, afac, ldab, ainv, lda,
457 $ info )
458*
459* Compute RCONDC = 1/(norm(A) * norm(inv(A))).
460*
461 anorm = dlansb( '1', uplo, n, kd, a, ldab, rwork )
462 ainvnm = dlange( '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 = 'DLARHS'
476 CALL dlarhs( path, xtype, uplo, ' ', n, n, kd,
477 $ kd, nrhs, a, ldab, xact, lda, b,
478 $ lda, iseed, info )
479 CALL dlacpy( 'Full', n, nrhs, b, lda, x, lda )
480*
481 srnamt = 'DPBTRS'
482 CALL dpbtrs( uplo, n, kd, nrhs, afac, ldab, x,
483 $ lda, info )
484*
485* Check error code from DPBTRS.
486*
487 IF( info.NE.0 )
488 $ CALL alaerh( path, 'DPBTRS', info, 0, uplo,
489 $ n, n, kd, kd, nrhs, imat, nfail,
490 $ nerrs, nout )
491*
492 CALL dlacpy( 'Full', n, nrhs, b, lda, work,
493 $ lda )
494 CALL dpbt02( 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 dget04( 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 = 'DPBRFS'
507 CALL dpbrfs( 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 DPBRFS.
513*
514 IF( info.NE.0 )
515 $ CALL alaerh( path, 'DPBRFS', info, 0, uplo,
516 $ n, n, kd, kd, nrhs, imat, nfail,
517 $ nerrs, nout )
518*
519 CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
520 $ result( 4 ) )
521 CALL dpbt05( 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 = 'DPBCON'
544 CALL dpbcon( uplo, n, kd, afac, ldab, anorm, rcond,
545 $ work, iwork, info )
546*
547* Check error code from DPBCON.
548*
549 IF( info.NE.0 )
550 $ CALL alaerh( path, 'DPBCON', info, 0, uplo, n,
551 $ n, kd, kd, -1, imat, nfail, nerrs,
552 $ nout )
553*
554 result( 7 ) = dget06( 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 DCHKPB
585*
alasum
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
Definition alasum.f:73
dlarhs
subroutine dlarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
DLARHS
Definition dlarhs.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
derrpo
subroutine derrpo(path, nunit)
DERRPO
Definition derrpo.f:55
dget04
subroutine dget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
DGET04
Definition dget04.f:102
dget06
double precision function dget06(rcond, rcondc)
DGET06
Definition dget06.f:55
dlatb4
subroutine dlatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
DLATB4
Definition dlatb4.f:120
dlatms
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
dpbt01
subroutine dpbt01(uplo, n, kd, a, lda, afac, ldafac, rwork, resid)
DPBT01
Definition dpbt01.f:119
dpbt02
subroutine dpbt02(uplo, n, kd, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DPBT02
Definition dpbt02.f:136
dpbt05
subroutine dpbt05(uplo, n, kd, nrhs, ab, ldab, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
DPBT05
Definition dpbt05.f:171
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:114
dlansb
double precision function dlansb(norm, uplo, n, k, ab, ldab, work)
DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansb.f:129
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dpbcon
subroutine dpbcon(uplo, n, kd, ab, ldab, anorm, rcond, work, iwork, info)
DPBCON
Definition dpbcon.f:132
dpbrfs
subroutine dpbrfs(uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x, ldx, ferr, berr, work, iwork, info)
DPBRFS
Definition dpbrfs.f:189
dpbtrf
subroutine dpbtrf(uplo, n, kd, ab, ldab, info)
DPBTRF
Definition dpbtrf.f:142
dpbtrs
subroutine dpbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
DPBTRS
Definition dpbtrs.f:121
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: