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