172
173
174
175
176
177
178 LOGICAL TSTERR
179 INTEGER LA, LAFB, NN, NOUT, NRHS
180 DOUBLE PRECISION THRESH
181
182
183 LOGICAL DOTYPE( * )
184 INTEGER IWORK( * ), NVAL( * )
185 DOUBLE PRECISION A( * ), AFB( * ), ASAV( * ), B( * ), BSAV( * ),
186 $ RWORK( * ), S( * ), WORK( * ), X( * ),
187 $ XACT( * )
188
189
190
191
192
193 DOUBLE PRECISION ONE, ZERO
194 parameter( one = 1.0d+0, zero = 0.0d+0 )
195 INTEGER NTYPES
196 parameter( ntypes = 8 )
197 INTEGER NTESTS
198 parameter( ntests = 7 )
199 INTEGER NTRAN
200 parameter( ntran = 3 )
201
202
203 LOGICAL EQUIL, NOFACT, PREFAC, TRFCON, ZEROT
204 CHARACTER DIST, EQUED, FACT, TRANS, TYPE, XTYPE
205 CHARACTER*3 PATH
206 INTEGER I, I1, I2, IEQUED, IFACT, IKL, IKU, IMAT, IN,
207 $ INFO, IOFF, ITRAN, IZERO, J, K, K1, KL, KU,
208 $ LDA, LDAFB, LDB, MODE, N, NB, NBMIN, NERRS,
209 $ NFACT, NFAIL, NIMAT, NKL, NKU, NRUN, NT
210 DOUBLE PRECISION AINVNM, AMAX, ANORM, ANORMI, ANORMO, ANRMPV,
211 $ CNDNUM, COLCND, RCOND, RCONDC, RCONDI, RCONDO,
212 $ ROLDC, ROLDI, ROLDO, ROWCND, RPVGRW
213
214
215 CHARACTER EQUEDS( 4 ), FACTS( 3 ), TRANSS( NTRAN )
216 INTEGER ISEED( 4 ), ISEEDY( 4 )
217 DOUBLE PRECISION RESULT( NTESTS )
218
219
220 LOGICAL LSAME
221 DOUBLE PRECISION DGET06, DLAMCH, DLANGB, DLANGE, DLANTB
223
224
229
230
231 INTRINSIC abs, max, min
232
233
234 LOGICAL LERR, OK
235 CHARACTER*32 SRNAMT
236 INTEGER INFOT, NUNIT
237
238
239 COMMON / infoc / infot, nunit, ok, lerr
240 COMMON / srnamc / srnamt
241
242
243 DATA iseedy / 1988, 1989, 1990, 1991 /
244 DATA transs / 'N', 'T', 'C' /
245 DATA facts / 'F', 'N', 'E' /
246 DATA equeds / 'N', 'R', 'C', 'B' /
247
248
249
250
251
252 path( 1: 1 ) = 'Double precision'
253 path( 2: 3 ) = 'GB'
254 nrun = 0
255 nfail = 0
256 nerrs = 0
257 DO 10 i = 1, 4
258 iseed( i ) = iseedy( i )
259 10 CONTINUE
260
261
262
263 IF( tsterr )
264 $
CALL derrvx( path, nout )
265 infot = 0
266
267
268
269 nb = 1
270 nbmin = 2
273
274
275
276 DO 150 in = 1, nn
277 n = nval( in )
278 ldb = max( n, 1 )
279 xtype = 'N'
280
281
282
283 nkl = max( 1, min( n, 4 ) )
284 IF( n.EQ.0 )
285 $ nkl = 1
286 nku = nkl
287 nimat = ntypes
288 IF( n.LE.0 )
289 $ nimat = 1
290
291 DO 140 ikl = 1, nkl
292
293
294
295
296 IF( ikl.EQ.1 ) THEN
297 kl = 0
298 ELSE IF( ikl.EQ.2 ) THEN
299 kl = max( n-1, 0 )
300 ELSE IF( ikl.EQ.3 ) THEN
301 kl = ( 3*n-1 ) / 4
302 ELSE IF( ikl.EQ.4 ) THEN
303 kl = ( n+1 ) / 4
304 END IF
305 DO 130 iku = 1, nku
306
307
308
309
310
311 IF( iku.EQ.1 ) THEN
312 ku = 0
313 ELSE IF( iku.EQ.2 ) THEN
314 ku = max( n-1, 0 )
315 ELSE IF( iku.EQ.3 ) THEN
316 ku = ( 3*n-1 ) / 4
317 ELSE IF( iku.EQ.4 ) THEN
318 ku = ( n+1 ) / 4
319 END IF
320
321
322
323
324 lda = kl + ku + 1
325 ldafb = 2*kl + ku + 1
326 IF( lda*n.GT.la .OR. ldafb*n.GT.lafb ) THEN
327 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
328 $
CALL aladhd( nout, path )
329 IF( lda*n.GT.la ) THEN
330 WRITE( nout, fmt = 9999 )la, n, kl, ku,
331 $ n*( kl+ku+1 )
332 nerrs = nerrs + 1
333 END IF
334 IF( ldafb*n.GT.lafb ) THEN
335 WRITE( nout, fmt = 9998 )lafb, n, kl, ku,
336 $ n*( 2*kl+ku+1 )
337 nerrs = nerrs + 1
338 END IF
339 GO TO 130
340 END IF
341
342 DO 120 imat = 1, nimat
343
344
345
346 IF( .NOT.dotype( imat ) )
347 $ GO TO 120
348
349
350
351 zerot = imat.GE.2 .AND. imat.LE.4
352 IF( zerot .AND. n.LT.imat-1 )
353 $ GO TO 120
354
355
356
357
358 CALL dlatb4( path, imat, n, n,
TYPE, KL, KU, ANORM,
359 $ MODE, CNDNUM, DIST )
360 rcondc = one / cndnum
361
362 srnamt = 'DLATMS'
363 CALL dlatms( n, n, dist, iseed,
TYPE, RWORK, MODE,
364 $ CNDNUM, ANORM, KL, KU, 'Z', A, LDA, WORK,
365 $ INFO )
366
367
368
369 IF( info.NE.0 ) THEN
370 CALL alaerh( path,
'DLATMS', info, 0,
' ', n, n,
371 $ kl, ku, -1, imat, nfail, nerrs, nout )
372 GO TO 120
373 END IF
374
375
376
377
378 izero = 0
379 IF( zerot ) THEN
380 IF( imat.EQ.2 ) THEN
381 izero = 1
382 ELSE IF( imat.EQ.3 ) THEN
383 izero = n
384 ELSE
385 izero = n / 2 + 1
386 END IF
387 ioff = ( izero-1 )*lda
388 IF( imat.LT.4 ) THEN
389 i1 = max( 1, ku+2-izero )
390 i2 = min( kl+ku+1, ku+1+( n-izero ) )
391 DO 20 i = i1, i2
392 a( ioff+i ) = zero
393 20 CONTINUE
394 ELSE
395 DO 40 j = izero, n
396 DO 30 i = max( 1, ku+2-j ),
397 $ min( kl+ku+1, ku+1+( n-j ) )
398 a( ioff+i ) = zero
399 30 CONTINUE
400 ioff = ioff + lda
401 40 CONTINUE
402 END IF
403 END IF
404
405
406
407 CALL dlacpy(
'Full', kl+ku+1, n, a, lda, asav, lda )
408
409 DO 110 iequed = 1, 4
410 equed = equeds( iequed )
411 IF( iequed.EQ.1 ) THEN
412 nfact = 3
413 ELSE
414 nfact = 1
415 END IF
416
417 DO 100 ifact = 1, nfact
418 fact = facts( ifact )
419 prefac =
lsame( fact,
'F' )
420 nofact =
lsame( fact,
'N' )
421 equil =
lsame( fact,
'E' )
422
423 IF( zerot ) THEN
424 IF( prefac )
425 $ GO TO 100
426 rcondo = zero
427 rcondi = zero
428
429 ELSE IF( .NOT.nofact ) THEN
430
431
432
433
434
435
436 CALL dlacpy(
'Full', kl+ku+1, n, asav, lda,
437 $ afb( kl+1 ), ldafb )
438 IF( equil .OR. iequed.GT.1 ) THEN
439
440
441
442
443 CALL dgbequ( n, n, kl, ku, afb( kl+1 ),
444 $ ldafb, s, s( n+1 ), rowcnd,
445 $ colcnd, amax, info )
446 IF( info.EQ.0 .AND. n.GT.0 ) THEN
447 IF(
lsame( equed,
'R' ) )
THEN
448 rowcnd = zero
449 colcnd = one
450 ELSE IF(
lsame( equed,
'C' ) )
THEN
451 rowcnd = one
452 colcnd = zero
453 ELSE IF(
lsame( equed,
'B' ) )
THEN
454 rowcnd = zero
455 colcnd = zero
456 END IF
457
458
459
460 CALL dlaqgb( n, n, kl, ku, afb( kl+1 ),
461 $ ldafb, s, s( n+1 ),
462 $ rowcnd, colcnd, amax,
463 $ equed )
464 END IF
465 END IF
466
467
468
469
470 IF( equil ) THEN
471 roldo = rcondo
472 roldi = rcondi
473 END IF
474
475
476
477 anormo =
dlangb(
'1', n, kl, ku, afb( kl+1 ),
478 $ ldafb, rwork )
479 anormi =
dlangb(
'I', n, kl, ku, afb( kl+1 ),
480 $ ldafb, rwork )
481
482
483
484 CALL dgbtrf( n, n, kl, ku, afb, ldafb, iwork,
485 $ info )
486
487
488
489 CALL dlaset(
'Full', n, n, zero, one, work,
490 $ ldb )
491 srnamt = 'DGBTRS'
492 CALL dgbtrs(
'No transpose', n, kl, ku, n,
493 $ afb, ldafb, iwork, work, ldb,
494 $ info )
495
496
497
498 ainvnm =
dlange(
'1', n, n, work, ldb,
499 $ rwork )
500 IF( anormo.LE.zero .OR. ainvnm.LE.zero ) THEN
501 rcondo = one
502 ELSE
503 rcondo = ( one / anormo ) / ainvnm
504 END IF
505
506
507
508
509 ainvnm =
dlange(
'I', n, n, work, ldb,
510 $ rwork )
511 IF( anormi.LE.zero .OR. ainvnm.LE.zero ) THEN
512 rcondi = one
513 ELSE
514 rcondi = ( one / anormi ) / ainvnm
515 END IF
516 END IF
517
518 DO 90 itran = 1, ntran
519
520
521
522 trans = transs( itran )
523 IF( itran.EQ.1 ) THEN
524 rcondc = rcondo
525 ELSE
526 rcondc = rcondi
527 END IF
528
529
530
531 CALL dlacpy(
'Full', kl+ku+1, n, asav, lda,
532 $ a, lda )
533
534
535
536
537 srnamt = 'DLARHS'
538 CALL dlarhs( path, xtype,
'Full', trans, n,
539 $ n, kl, ku, nrhs, a, lda, xact,
540 $ ldb, b, ldb, iseed, info )
541 xtype = 'C'
542 CALL dlacpy(
'Full', n, nrhs, b, ldb, bsav,
543 $ ldb )
544
545 IF( nofact .AND. itran.EQ.1 ) THEN
546
547
548
549
550
551
552 CALL dlacpy(
'Full', kl+ku+1, n, a, lda,
553 $ afb( kl+1 ), ldafb )
554 CALL dlacpy(
'Full', n, nrhs, b, ldb, x,
555 $ ldb )
556
557 srnamt = 'DGBSV '
558 CALL dgbsv( n, kl, ku, nrhs, afb, ldafb,
559 $ iwork, x, ldb, info )
560
561
562
563 IF( info.NE.izero )
564 $
CALL alaerh( path,
'DGBSV ', info,
565 $ izero, ' ', n, n, kl, ku,
566 $ nrhs, imat, nfail, nerrs,
567 $ nout )
568
569
570
571
572 CALL dgbt01( n, n, kl, ku, a, lda, afb,
573 $ ldafb, iwork, work,
574 $ result( 1 ) )
575 nt = 1
576 IF( izero.EQ.0 ) THEN
577
578
579
580
581 CALL dlacpy(
'Full', n, nrhs, b, ldb,
582 $ work, ldb )
583 CALL dgbt02(
'No transpose', n, n, kl,
584 $ ku, nrhs, a, lda, x, ldb,
585 $ work, ldb, rwork,
586 $ result( 2 ) )
587
588
589
590
591 CALL dget04( n, nrhs, x, ldb, xact,
592 $ ldb, rcondc, result( 3 ) )
593 nt = 3
594 END IF
595
596
597
598
599 DO 50 k = 1, nt
600 IF( result( k ).GE.thresh ) THEN
601 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
602 $
CALL aladhd( nout, path )
603 WRITE( nout, fmt = 9997 )'DGBSV ',
604 $ n, kl, ku, imat, k, result( k )
605 nfail = nfail + 1
606 END IF
607 50 CONTINUE
608 nrun = nrun + nt
609 END IF
610
611
612
613 IF( .NOT.prefac )
614 $
CALL dlaset(
'Full', 2*kl+ku+1, n, zero,
615 $ zero, afb, ldafb )
616 CALL dlaset(
'Full', n, nrhs, zero, zero, x,
617 $ ldb )
618 IF( iequed.GT.1 .AND. n.GT.0 ) THEN
619
620
621
622
623 CALL dlaqgb( n, n, kl, ku, a, lda, s,
624 $ s( n+1 ), rowcnd, colcnd,
625 $ amax, equed )
626 END IF
627
628
629
630
631 srnamt = 'DGBSVX'
632 CALL dgbsvx( fact, trans, n, kl, ku, nrhs, a,
633 $ lda, afb, ldafb, iwork, equed,
634 $ s, s( n+1 ), b, ldb, x, ldb,
635 $ rcond, rwork, rwork( nrhs+1 ),
636 $ work, iwork( n+1 ), info )
637
638
639
640 IF( info.NE.izero )
641 $
CALL alaerh( path,
'DGBSVX', info, izero,
642 $ fact // trans, n, n, kl, ku,
643 $ nrhs, imat, nfail, nerrs,
644 $ nout )
645
646
647
648
649 IF( info.NE.0 .AND. info.LE.n) THEN
650 anrmpv = zero
651 DO 70 j = 1, info
652 DO 60 i = max( ku+2-j, 1 ),
653 $ min( n+ku+1-j, kl+ku+1 )
654 anrmpv = max( anrmpv,
655 $ abs( a( i+( j-1 )*lda ) ) )
656 60 CONTINUE
657 70 CONTINUE
658 rpvgrw =
dlantb(
'M',
'U',
'N', info,
659 $ min( info-1, kl+ku ),
660 $ afb( max( 1, kl+ku+2-info ) ),
661 $ ldafb, work )
662 IF( rpvgrw.EQ.zero ) THEN
663 rpvgrw = one
664 ELSE
665 rpvgrw = anrmpv / rpvgrw
666 END IF
667 ELSE
668 rpvgrw =
dlantb(
'M',
'U',
'N', n, kl+ku,
669 $ afb, ldafb, work )
670 IF( rpvgrw.EQ.zero ) THEN
671 rpvgrw = one
672 ELSE
673 rpvgrw =
dlangb(
'M', n, kl, ku, a,
674 $ lda, work ) / rpvgrw
675 END IF
676 END IF
677 result( 7 ) = abs( rpvgrw-work( 1 ) ) /
678 $ max( work( 1 ), rpvgrw ) /
680
681 IF( .NOT.prefac ) THEN
682
683
684
685
686 CALL dgbt01( n, n, kl, ku, a, lda, afb,
687 $ ldafb, iwork, work,
688 $ result( 1 ) )
689 k1 = 1
690 ELSE
691 k1 = 2
692 END IF
693
694 IF( info.EQ.0 ) THEN
695 trfcon = .false.
696
697
698
699 CALL dlacpy(
'Full', n, nrhs, bsav, ldb,
700 $ work, ldb )
701 CALL dgbt02( trans, n, n, kl, ku, nrhs,
702 $ asav, lda, x, ldb, work, ldb,
703 $ rwork( 2*nrhs+1 ),
704 $ result( 2 ) )
705
706
707
708
709 IF( nofact .OR. ( prefac .AND.
710 $
lsame( equed,
'N' ) ) )
THEN
711 CALL dget04( n, nrhs, x, ldb, xact,
712 $ ldb, rcondc, result( 3 ) )
713 ELSE
714 IF( itran.EQ.1 ) THEN
715 roldc = roldo
716 ELSE
717 roldc = roldi
718 END IF
719 CALL dget04( n, nrhs, x, ldb, xact,
720 $ ldb, roldc, result( 3 ) )
721 END IF
722
723
724
725
726 CALL dgbt05( trans, n, kl, ku, nrhs, asav,
727 $ lda, b, ldb, x, ldb, xact,
728 $ ldb, rwork, rwork( nrhs+1 ),
729 $ result( 4 ) )
730 ELSE
731 trfcon = .true.
732 END IF
733
734
735
736
737 result( 6 ) =
dget06( rcond, rcondc )
738
739
740
741
742 IF( .NOT.trfcon ) THEN
743 DO 80 k = k1, ntests
744 IF( result( k ).GE.thresh ) THEN
745 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
746 $
CALL aladhd( nout, path )
747 IF( prefac ) THEN
748 WRITE( nout, fmt = 9995 )
749 $ 'DGBSVX', fact, trans, n, kl,
750 $ ku, equed, imat, k,
751 $ result( k )
752 ELSE
753 WRITE( nout, fmt = 9996 )
754 $ 'DGBSVX', fact, trans, n, kl,
755 $ ku, imat, k, result( k )
756 END IF
757 nfail = nfail + 1
758 END IF
759 80 CONTINUE
760 nrun = nrun + ntests - k1 + 1
761 ELSE
762 IF( result( 1 ).GE.thresh .AND. .NOT.
763 $ prefac ) THEN
764 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
765 $
CALL aladhd( nout, path )
766 IF( prefac ) THEN
767 WRITE( nout, fmt = 9995 )'DGBSVX',
768 $ fact, trans, n, kl, ku, equed,
769 $ imat, 1, result( 1 )
770 ELSE
771 WRITE( nout, fmt = 9996 )'DGBSVX',
772 $ fact, trans, n, kl, ku, imat, 1,
773 $ result( 1 )
774 END IF
775 nfail = nfail + 1
776 nrun = nrun + 1
777 END IF
778 IF( result( 6 ).GE.thresh ) THEN
779 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
780 $
CALL aladhd( nout, path )
781 IF( prefac ) THEN
782 WRITE( nout, fmt = 9995 )'DGBSVX',
783 $ fact, trans, n, kl, ku, equed,
784 $ imat, 6, result( 6 )
785 ELSE
786 WRITE( nout, fmt = 9996 )'DGBSVX',
787 $ fact, trans, n, kl, ku, imat, 6,
788 $ result( 6 )
789 END IF
790 nfail = nfail + 1
791 nrun = nrun + 1
792 END IF
793 IF( result( 7 ).GE.thresh ) THEN
794 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
795 $
CALL aladhd( nout, path )
796 IF( prefac ) THEN
797 WRITE( nout, fmt = 9995 )'DGBSVX',
798 $ fact, trans, n, kl, ku, equed,
799 $ imat, 7, result( 7 )
800 ELSE
801 WRITE( nout, fmt = 9996 )'DGBSVX',
802 $ fact, trans, n, kl, ku, imat, 7,
803 $ result( 7 )
804 END IF
805 nfail = nfail + 1
806 nrun = nrun + 1
807 END IF
808
809 END IF
810 90 CONTINUE
811 100 CONTINUE
812 110 CONTINUE
813 120 CONTINUE
814 130 CONTINUE
815 140 CONTINUE
816 150 CONTINUE
817
818
819
820 CALL alasvm( path, nout, nfail, nrun, nerrs )
821
822 9999 FORMAT( ' *** In DDRVGB, LA=', i5, ' is too small for N=', i5,
823 $ ', KU=', i5, ', KL=', i5, / ' ==> Increase LA to at least ',
824 $ i5 )
825 9998 FORMAT( ' *** In DDRVGB, LAFB=', i5, ' is too small for N=', i5,
826 $ ', KU=', i5, ', KL=', i5, /
827 $ ' ==> Increase LAFB to at least ', i5 )
828 9997 FORMAT( 1x, a, ', N=', i5, ', KL=', i5, ', KU=', i5, ', type ',
829 $ i1, ', test(', i1, ')=', g12.5 )
830 9996 FORMAT( 1x, a, '( ''', a1, ''',''', a1, ''',', i5, ',', i5, ',',
831 $ i5, ',...), type ', i1, ', test(', i1, ')=', g12.5 )
832 9995 FORMAT( 1x, a, '( ''', a1, ''',''', a1, ''',', i5, ',', i5, ',',
833 $ i5, ',...), EQUED=''', a1, ''', type ', i1, ', test(', i1,
834 $ ')=', g12.5 )
835
836 RETURN
837
838
839
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine dlarhs(path, xtype, uplo, trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, iseed, info)
DLARHS
subroutine xlaenv(ispec, nvalue)
XLAENV
subroutine aladhd(iounit, path)
ALADHD
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
subroutine derrvx(path, nunit)
DERRVX
subroutine dgbt01(m, n, kl, ku, a, lda, afac, ldafac, ipiv, work, resid)
DGBT01
subroutine dgbt02(trans, m, n, kl, ku, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DGBT02
subroutine dgbt05(trans, n, kl, ku, nrhs, ab, ldab, b, ldb, x, ldx, xact, ldxact, ferr, berr, reslts)
DGBT05
subroutine dget04(n, nrhs, x, ldx, xact, ldxact, rcond, resid)
DGET04
double precision function dget06(rcond, rcondc)
DGET06
subroutine dlatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
DLATB4
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
subroutine dgbequ(m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
DGBEQU
subroutine dgbsv(n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
DGBSV computes the solution to system of linear equations A * X = B for GB matrices (simple driver)
subroutine dgbsvx(fact, trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, equed, r, c, b, ldb, x, ldx, rcond, ferr, berr, work, iwork, info)
DGBSVX computes the solution to system of linear equations A * X = B for GB matrices
subroutine dgbtrf(m, n, kl, ku, ab, ldab, ipiv, info)
DGBTRF
subroutine dgbtrs(trans, n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info)
DGBTRS
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function dlangb(norm, n, kl, ku, ab, ldab, work)
DLANGB returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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 ...
double precision function dlantb(norm, uplo, diag, n, k, ab, ldab, work)
DLANTB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine dlaqgb(m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, equed)
DLAQGB scales a general band matrix, using row and column scaling factors computed by sgbequ.
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.
logical function lsame(ca, cb)
LSAME