171
172
173
174
175
176
177 LOGICAL TSTERR
178 INTEGER NMAX, NN, NNB, NNS, NOUT
179 DOUBLE PRECISION THRESH
180
181
182 LOGICAL DOTYPE( * )
183 INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
184 DOUBLE PRECISION A( * ), AFAC( * ), AINV( * ), B( * ),
185 $ RWORK( * ), WORK( * ), X( * ), XACT( * )
186
187
188
189
190
191 DOUBLE PRECISION ZERO, ONE
192 parameter( zero = 0.0d+0, one = 1.0d+0 )
193 DOUBLE PRECISION EIGHT, SEVTEN
194 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
195 INTEGER NTYPES
196 parameter( ntypes = 10 )
197 INTEGER NTESTS
198 parameter( ntests = 7 )
199
200
201 LOGICAL TRFCON, ZEROT
202 CHARACTER DIST, TYPE, UPLO, XTYPE
203 CHARACTER*3 PATH, MATPATH
204 INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
205 $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
206 $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
207 DOUBLE PRECISION ALPHA, ANORM, CNDNUM, CONST, DTEMP, SING_MAX,
208 $ SING_MIN, RCOND, RCONDC
209
210
211 CHARACTER UPLOS( 2 )
212 INTEGER ISEED( 4 ), ISEEDY( 4 )
213 DOUBLE PRECISION BLOCK( 2, 2 ), DDUMMY( 1 ), RESULT( NTESTS )
214
215
216 DOUBLE PRECISION DGET06, DLANGE, DLANSY
218
219
224
225
226 INTRINSIC max, min, sqrt
227
228
229 LOGICAL LERR, OK
230 CHARACTER*32 SRNAMT
231 INTEGER INFOT, NUNIT
232
233
234 COMMON / infoc / infot, nunit, ok, lerr
235 COMMON / srnamc / srnamt
236
237
238 DATA iseedy / 1988, 1989, 1990, 1991 /
239 DATA uplos / 'U', 'L' /
240
241
242
243
244
245 alpha = ( one+sqrt( sevten ) ) / eight
246
247
248
249 path( 1: 1 ) = 'Double precision'
250 path( 2: 3 ) = 'SR'
251
252
253
254 matpath( 1: 1 ) = 'Double precision'
255 matpath( 2: 3 ) = 'SY'
256
257 nrun = 0
258 nfail = 0
259 nerrs = 0
260 DO 10 i = 1, 4
261 iseed( i ) = iseedy( i )
262 10 CONTINUE
263
264
265
266 IF( tsterr )
267 $
CALL derrsy( path, nout )
268 infot = 0
269
270
271
272
274
275
276
277 DO 270 in = 1, nn
278 n = nval( in )
279 lda = max( n, 1 )
280 xtype = 'N'
281 nimat = ntypes
282 IF( n.LE.0 )
283 $ nimat = 1
284
285 izero = 0
286
287
288
289 DO 260 imat = 1, nimat
290
291
292
293 IF( .NOT.dotype( imat ) )
294 $ GO TO 260
295
296
297
298 zerot = imat.GE.3 .AND. imat.LE.6
299 IF( zerot .AND. n.LT.imat-2 )
300 $ GO TO 260
301
302
303
304 DO 250 iuplo = 1, 2
305 uplo = uplos( iuplo )
306
307
308
309
310
311
312 CALL dlatb4( matpath, imat, n, n,
TYPE, KL, KU, ANORM,
313 $ MODE, CNDNUM, DIST )
314
315
316
317 srnamt = 'DLATMS'
318 CALL dlatms( n, n, dist, iseed,
TYPE, RWORK, MODE,
319 $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
320 $ INFO )
321
322
323
324 IF( info.NE.0 ) THEN
325 CALL alaerh( path,
'DLATMS', info, 0, uplo, n, n, -1,
326 $ -1, -1, imat, nfail, nerrs, nout )
327
328
329
330 GO TO 250
331 END IF
332
333
334
335
336
337 IF( zerot ) THEN
338 IF( imat.EQ.3 ) THEN
339 izero = 1
340 ELSE IF( imat.EQ.4 ) THEN
341 izero = n
342 ELSE
343 izero = n / 2 + 1
344 END IF
345
346 IF( imat.LT.6 ) THEN
347
348
349
350 IF( iuplo.EQ.1 ) THEN
351 ioff = ( izero-1 )*lda
352 DO 20 i = 1, izero - 1
353 a( ioff+i ) = zero
354 20 CONTINUE
355 ioff = ioff + izero
356 DO 30 i = izero, n
357 a( ioff ) = zero
358 ioff = ioff + lda
359 30 CONTINUE
360 ELSE
361 ioff = izero
362 DO 40 i = 1, izero - 1
363 a( ioff ) = zero
364 ioff = ioff + lda
365 40 CONTINUE
366 ioff = ioff - izero
367 DO 50 i = izero, n
368 a( ioff+i ) = zero
369 50 CONTINUE
370 END IF
371 ELSE
372 IF( iuplo.EQ.1 ) THEN
373
374
375
376 ioff = 0
377 DO 70 j = 1, n
378 i2 = min( j, izero )
379 DO 60 i = 1, i2
380 a( ioff+i ) = zero
381 60 CONTINUE
382 ioff = ioff + lda
383 70 CONTINUE
384 ELSE
385
386
387
388 ioff = 0
389 DO 90 j = 1, n
390 i1 = max( j, izero )
391 DO 80 i = i1, n
392 a( ioff+i ) = zero
393 80 CONTINUE
394 ioff = ioff + lda
395 90 CONTINUE
396 END IF
397 END IF
398 ELSE
399 izero = 0
400 END IF
401
402
403
404
405
406
407 DO 240 inb = 1, nnb
408
409
410
411
412 nb = nbval( inb )
414
415
416
417
418
419 CALL dlacpy( uplo, n, n, a, lda, afac, lda )
420
421
422
423
424
425
426 lwork = max( 2, nb )*lda
427 srnamt = 'DSYTRF_ROOK'
429 $ lwork, info )
430
431
432
433
434 k = izero
435 IF( k.GT.0 ) THEN
436 100 CONTINUE
437 IF( iwork( k ).LT.0 ) THEN
438 IF( iwork( k ).NE.-k ) THEN
439 k = -iwork( k )
440 GO TO 100
441 END IF
442 ELSE IF( iwork( k ).NE.k ) THEN
443 k = iwork( k )
444 GO TO 100
445 END IF
446 END IF
447
448
449
450 IF( info.NE.k)
451 $
CALL alaerh( path,
'DSYTRF_ROOK', info, k,
452 $ uplo, n, n, -1, -1, nb, imat,
453 $ nfail, nerrs, nout )
454
455
456
457 IF( info.NE.0 ) THEN
458 trfcon = .true.
459 ELSE
460 trfcon = .false.
461 END IF
462
463
464
465
466 CALL dsyt01_rook( uplo, n, a, lda, afac, lda, iwork,
467 $ ainv, lda, rwork, result( 1 ) )
468 nt = 1
469
470
471
472
473
474
475
476 IF( inb.EQ.1 .AND. .NOT.trfcon ) THEN
477 CALL dlacpy( uplo, n, n, afac, lda, ainv, lda )
478 srnamt = 'DSYTRI_ROOK'
480 $ info )
481
482
483
484 IF( info.NE.0 )
485 $
CALL alaerh( path,
'DSYTRI_ROOK', info, -1,
486 $ uplo, n, n, -1, -1, -1, imat,
487 $ nfail, nerrs, nout )
488
489
490
491
492 CALL dpot03( uplo, n, a, lda, ainv, lda, work, lda,
493 $ rwork, rcondc, result( 2 ) )
494 nt = 2
495 END IF
496
497
498
499
500 DO 110 k = 1, nt
501 IF( result( k ).GE.thresh ) THEN
502 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
503 $
CALL alahd( nout, path )
504 WRITE( nout, fmt = 9999 )uplo, n, nb, imat, k,
505 $ result( k )
506 nfail = nfail + 1
507 END IF
508 110 CONTINUE
509 nrun = nrun + nt
510
511
512
513
514 result( 3 ) = zero
515 dtemp = zero
516
517 const = one / ( one-alpha )
518
519 IF( iuplo.EQ.1 ) THEN
520
521
522
523 k = n
524 120 CONTINUE
525 IF( k.LE.1 )
526 $ GO TO 130
527
528 IF( iwork( k ).GT.zero ) THEN
529
530
531
532
533 dtemp =
dlange(
'M', k-1, 1,
534 $ afac( ( k-1 )*lda+1 ), lda, rwork )
535 ELSE
536
537
538
539
540 dtemp =
dlange(
'M', k-2, 2,
541 $ afac( ( k-2 )*lda+1 ), lda, rwork )
542 k = k - 1
543
544 END IF
545
546
547
548 dtemp = dtemp - const + thresh
549 IF( dtemp.GT.result( 3 ) )
550 $ result( 3 ) = dtemp
551
552 k = k - 1
553
554 GO TO 120
555 130 CONTINUE
556
557 ELSE
558
559
560
561 k = 1
562 140 CONTINUE
563 IF( k.GE.n )
564 $ GO TO 150
565
566 IF( iwork( k ).GT.zero ) THEN
567
568
569
570
571 dtemp =
dlange(
'M', n-k, 1,
572 $ afac( ( k-1 )*lda+k+1 ), lda, rwork )
573 ELSE
574
575
576
577
578 dtemp =
dlange(
'M', n-k-1, 2,
579 $ afac( ( k-1 )*lda+k+2 ), lda, rwork )
580 k = k + 1
581
582 END IF
583
584
585
586 dtemp = dtemp - const + thresh
587 IF( dtemp.GT.result( 3 ) )
588 $ result( 3 ) = dtemp
589
590 k = k + 1
591
592 GO TO 140
593 150 CONTINUE
594 END IF
595
596
597
598
599
600
601 result( 4 ) = zero
602 dtemp = zero
603
604 const = ( one+alpha ) / ( one-alpha )
605 CALL dlacpy( uplo, n, n, afac, lda, ainv, lda )
606
607 IF( iuplo.EQ.1 ) THEN
608
609
610
611 k = n
612 160 CONTINUE
613 IF( k.LE.1 )
614 $ GO TO 170
615
616 IF( iwork( k ).LT.zero ) THEN
617
618
619
620
621
622 block( 1, 1 ) = afac( ( k-2 )*lda+k-1 )
623 block( 1, 2 ) = afac( (k-1)*lda+k-1 )
624 block( 2, 1 ) = block( 1, 2 )
625 block( 2, 2 ) = afac( (k-1)*lda+k )
626
627 CALL dgesvd(
'N',
'N', 2, 2, block, 2, rwork,
628 $ ddummy, 1, ddummy, 1,
629 $ work, 10, info )
630
631 sing_max = rwork( 1 )
632 sing_min = rwork( 2 )
633
634 dtemp = sing_max / sing_min
635
636
637
638 dtemp = dtemp - const + thresh
639 IF( dtemp.GT.result( 4 ) )
640 $ result( 4 ) = dtemp
641 k = k - 1
642
643 END IF
644
645 k = k - 1
646
647 GO TO 160
648 170 CONTINUE
649
650 ELSE
651
652
653
654 k = 1
655 180 CONTINUE
656 IF( k.GE.n )
657 $ GO TO 190
658
659 IF( iwork( k ).LT.zero ) THEN
660
661
662
663
664
665 block( 1, 1 ) = afac( ( k-1 )*lda+k )
666 block( 2, 1 ) = afac( ( k-1 )*lda+k+1 )
667 block( 1, 2 ) = block( 2, 1 )
668 block( 2, 2 ) = afac( k*lda+k+1 )
669
670 CALL dgesvd(
'N',
'N', 2, 2, block, 2, rwork,
671 $ ddummy, 1, ddummy, 1,
672 $ work, 10, info )
673
674
675 sing_max = rwork( 1 )
676 sing_min = rwork( 2 )
677
678 dtemp = sing_max / sing_min
679
680
681
682 dtemp = dtemp - const + thresh
683 IF( dtemp.GT.result( 4 ) )
684 $ result( 4 ) = dtemp
685 k = k + 1
686
687 END IF
688
689 k = k + 1
690
691 GO TO 180
692 190 CONTINUE
693 END IF
694
695
696
697
698 DO 200 k = 3, 4
699 IF( result( k ).GE.thresh ) THEN
700 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
701 $
CALL alahd( nout, path )
702 WRITE( nout, fmt = 9999 )uplo, n, nb, imat, k,
703 $ result( k )
704 nfail = nfail + 1
705 END IF
706 200 CONTINUE
707 nrun = nrun + 2
708
709
710
711
712 IF( inb.GT.1 )
713 $ GO TO 240
714
715
716
717 IF( trfcon ) THEN
718 rcondc = zero
719 GO TO 230
720 END IF
721
722
723
724 DO 220 irhs = 1, nns
725 nrhs = nsval( irhs )
726
727
728
729
730
731
732
733 srnamt = 'DLARHS'
734 CALL dlarhs( matpath, xtype, uplo,
' ', n, n,
735 $ kl, ku, nrhs, a, lda, xact, lda,
736 $ b, lda, iseed, info )
737 CALL dlacpy(
'Full', n, nrhs, b, lda, x, lda )
738
739 srnamt = 'DSYTRS_ROOK'
741 $ x, lda, info )
742
743
744
745 IF( info.NE.0 )
746 $
CALL alaerh( path,
'DSYTRS_ROOK', info, 0,
747 $ uplo, n, n, -1, -1, nrhs, imat,
748 $ nfail, nerrs, nout )
749
750 CALL dlacpy(
'Full', n, nrhs, b, lda, work, lda )
751
752
753
754 CALL dpot02( uplo, n, nrhs, a, lda, x, lda, work,
755 $ lda, rwork, result( 5 ) )
756
757
758
759
760 CALL dget04( n, nrhs, x, lda, xact, lda, rcondc,
761 $ result( 6 ) )
762
763
764
765
766 DO 210 k = 5, 6
767 IF( result( k ).GE.thresh ) THEN
768 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
769 $
CALL alahd( nout, path )
770 WRITE( nout, fmt = 9998 )uplo, n, nrhs,
771 $ imat, k, result( k )
772 nfail = nfail + 1
773 END IF
774 210 CONTINUE
775 nrun = nrun + 2
776
777
778
779 220 CONTINUE
780
781
782
783
784 230 CONTINUE
785 anorm =
dlansy(
'1', uplo, n, a, lda, rwork )
786 srnamt = 'DSYCON_ROOK'
787 CALL dsycon_rook( uplo, n, afac, lda, iwork, anorm,
788 $ rcond, work, iwork( n+1 ), info )
789
790
791
792 IF( info.NE.0 )
793 $
CALL alaerh( path,
'DSYCON_ROOK', info, 0,
794 $ uplo, n, n, -1, -1, -1, imat,
795 $ nfail, nerrs, nout )
796
797
798
799 result( 7 ) =
dget06( rcond, rcondc )
800
801
802
803
804 IF( result( 7 ).GE.thresh ) THEN
805 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
806 $
CALL alahd( nout, path )
807 WRITE( nout, fmt = 9997 )uplo, n, imat, 7,
808 $ result( 7 )
809 nfail = nfail + 1
810 END IF
811 nrun = nrun + 1
812 240 CONTINUE
813
814 250 CONTINUE
815 260 CONTINUE
816 270 CONTINUE
817
818
819
820 CALL alasum( path, nout, nfail, nrun, nerrs )
821
822 9999 FORMAT( ' UPLO = ''', a1, ''', N =', i5, ', NB =', i4, ', type ',
823 $ i2, ', test ', i2, ', ratio =', g12.5 )
824 9998 FORMAT( ' UPLO = ''', a1, ''', N =', i5, ', NRHS=', i3, ', type ',
825 $ i2, ', test(', i2, ') =', g12.5 )
826 9997 FORMAT( ' UPLO = ''', a1, ''', N =', i5, ',', 10x, ' type ', i2,
827 $ ', test(', i2, ') =', g12.5 )
828 RETURN
829
830
831
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
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 alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
subroutine alahd(iounit, path)
ALAHD
subroutine derrsy(path, nunit)
DERRSY
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 dpot02(uplo, n, nrhs, a, lda, x, ldx, b, ldb, rwork, resid)
DPOT02
subroutine dpot03(uplo, n, a, lda, ainv, ldainv, work, ldwork, rwork, rcond, resid)
DPOT03
subroutine dsyt01_rook(uplo, n, a, lda, afac, ldafac, ipiv, c, ldc, rwork, resid)
DSYT01_ROOK
subroutine dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
DGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine dsycon_rook(uplo, n, a, lda, ipiv, anorm, rcond, work, iwork, info)
DSYCON_ROOK
subroutine dsytrf_rook(uplo, n, a, lda, ipiv, work, lwork, info)
DSYTRF_ROOK
subroutine dsytri_rook(uplo, n, a, lda, ipiv, work, info)
DSYTRI_ROOK
subroutine dsytrs_rook(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
DSYTRS_ROOK
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
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 dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...