301
302
303
304
305
306
307 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
308 $ NSIZE
309 REAL THRESH
310
311
312 LOGICAL BWORK( * )
313 INTEGER IWORK( * )
314 REAL A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
315 $ ALPHAR( * ), B( LDA, * ), BETA( * ),
316 $ BI( LDA, * ), DIF( * ), DIFTRU( * ),
317 $ LSCALE( * ), RESULT( 4 ), RSCALE( * ), S( * ),
318 $ STRU( * ), VL( LDA, * ), VR( LDA, * ),
319 $ WORK( * )
320
321
322
323
324
325 REAL ZERO, ONE, TEN, TNTH, HALF
326 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
327 $ tnth = 1.0e-1, half = 0.5d+0 )
328
329
330 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
331 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
332 REAL ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
333 $ ULP, ULPINV
334
335
336 REAL WEIGHT( 5 )
337
338
339 INTEGER ILAENV
340 REAL SLAMCH, SLANGE
342
343
345
346
347 INTRINSIC abs, max, sqrt
348
349
350
351
352
353 info = 0
354
355 nmax = 5
356
357 IF( nsize.LT.0 ) THEN
358 info = -1
359 ELSE IF( thresh.LT.zero ) THEN
360 info = -2
361 ELSE IF( nin.LE.0 ) THEN
362 info = -3
363 ELSE IF( nout.LE.0 ) THEN
364 info = -4
365 ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
366 info = -6
367 ELSE IF( liwork.LT.nmax+6 ) THEN
368 info = -26
369 END IF
370
371
372
373
374
375
376
377
378 minwrk = 1
379 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
380 minwrk = 2*nmax*nmax + 12*nmax + 16
381 maxwrk = 6*nmax + nmax*
ilaenv( 1,
'SGEQRF',
' ', nmax, 1, nmax,
382 $ 0 )
383 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
384 work( 1 ) = maxwrk
385 END IF
386
387 IF( lwork.LT.minwrk )
388 $ info = -24
389
390 IF( info.NE.0 ) THEN
391 CALL xerbla(
'SDRGVX', -info )
392 RETURN
393 END IF
394
395 n = 5
397 ulpinv = one / ulp
398 thrsh2 = ten*thresh
399 nerrs = 0
400 nptknt = 0
401 ntestt = 0
402
403 IF( nsize.EQ.0 )
404 $ GO TO 90
405
406
407
408 weight( 1 ) = tnth
409 weight( 2 ) = half
410 weight( 3 ) = one
411 weight( 4 ) = one / weight( 2 )
412 weight( 5 ) = one / weight( 1 )
413
414 DO 80 iptype = 1, 2
415 DO 70 iwa = 1, 5
416 DO 60 iwb = 1, 5
417 DO 50 iwx = 1, 5
418 DO 40 iwy = 1, 5
419
420
421
422 CALL slatm6( iptype, 5, a, lda, b, vr, lda, vl,
423 $ lda, weight( iwa ), weight( iwb ),
424 $ weight( iwx ), weight( iwy ), stru,
425 $ diftru )
426
427
428
429
430
431 CALL slacpy(
'F', n, n, a, lda, ai, lda )
432 CALL slacpy(
'F', n, n, b, lda, bi, lda )
433
434 CALL sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
435 $ lda, alphar, alphai, beta, vl, lda,
436 $ vr, lda, ilo, ihi, lscale, rscale,
437 $ anorm, bnorm, s, dif, work, lwork,
438 $ iwork, bwork, linfo )
439 IF( linfo.NE.0 ) THEN
440 result( 1 ) = ulpinv
441 WRITE( nout, fmt = 9999 )'SGGEVX', linfo, n,
442 $ iptype
443 GO TO 30
444 END IF
445
446
447
448 CALL slacpy(
'Full', n, n, ai, lda, work, n )
449 CALL slacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
450 $ n )
451 abnorm =
slange(
'Fro', n, 2*n, work, n, work )
452
453
454
455 result( 1 ) = zero
456 CALL sget52( .true., n, a, lda, b, lda, vl, lda,
457 $ alphar, alphai, beta, work,
458 $ result( 1 ) )
459 IF( result( 2 ).GT.thresh ) THEN
460 WRITE( nout, fmt = 9998 )'Left', 'SGGEVX',
461 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
462 END IF
463
464 result( 2 ) = zero
465 CALL sget52( .false., n, a, lda, b, lda, vr, lda,
466 $ alphar, alphai, beta, work,
467 $ result( 2 ) )
468 IF( result( 3 ).GT.thresh ) THEN
469 WRITE( nout, fmt = 9998 )'Right', 'SGGEVX',
470 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
471 END IF
472
473
474
475 result( 3 ) = zero
476 DO 10 i = 1, n
477 IF( s( i ).EQ.zero ) THEN
478 IF( stru( i ).GT.abnorm*ulp )
479 $ result( 3 ) = ulpinv
480 ELSE IF( stru( i ).EQ.zero ) THEN
481 IF( s( i ).GT.abnorm*ulp )
482 $ result( 3 ) = ulpinv
483 ELSE
484 work( i ) = max( abs( stru( i ) / s( i ) ),
485 $ abs( s( i ) / stru( i ) ) )
486 result( 3 ) = max( result( 3 ), work( i ) )
487 END IF
488 10 CONTINUE
489
490
491
492 result( 4 ) = zero
493 IF( dif( 1 ).EQ.zero ) THEN
494 IF( diftru( 1 ).GT.abnorm*ulp )
495 $ result( 4 ) = ulpinv
496 ELSE IF( diftru( 1 ).EQ.zero ) THEN
497 IF( dif( 1 ).GT.abnorm*ulp )
498 $ result( 4 ) = ulpinv
499 ELSE IF( dif( 5 ).EQ.zero ) THEN
500 IF( diftru( 5 ).GT.abnorm*ulp )
501 $ result( 4 ) = ulpinv
502 ELSE IF( diftru( 5 ).EQ.zero ) THEN
503 IF( dif( 5 ).GT.abnorm*ulp )
504 $ result( 4 ) = ulpinv
505 ELSE
506 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
507 $ abs( dif( 1 ) / diftru( 1 ) ) )
508 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
509 $ abs( dif( 5 ) / diftru( 5 ) ) )
510 result( 4 ) = max( ratio1, ratio2 )
511 END IF
512
513 ntestt = ntestt + 4
514
515
516
517 DO 20 j = 1, 4
518 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
519 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
520 $ THEN
521
522
523
524
525 IF( nerrs.EQ.0 ) THEN
526 WRITE( nout, fmt = 9997 )'SXV'
527
528
529
530
531
532 WRITE( nout, fmt = 9995 )
533 WRITE( nout, fmt = 9994 )
534 WRITE( nout, fmt = 9993 )
535
536
537
538 WRITE( nout, fmt = 9992 )'''',
539 $ 'transpose', ''''
540
541 END IF
542 nerrs = nerrs + 1
543 IF( result( j ).LT.10000.0 ) THEN
544 WRITE( nout, fmt = 9991 )iptype, iwa,
545 $ iwb, iwx, iwy, j, result( j )
546 ELSE
547 WRITE( nout, fmt = 9990 )iptype, iwa,
548 $ iwb, iwx, iwy, j, result( j )
549 END IF
550 END IF
551 20 CONTINUE
552
553 30 CONTINUE
554
555 40 CONTINUE
556 50 CONTINUE
557 60 CONTINUE
558 70 CONTINUE
559 80 CONTINUE
560
561 GO TO 150
562
563 90 CONTINUE
564
565
566
567
568 READ( nin, fmt = *, END = 150 )n
569 IF( n.EQ.0 )
570 $ GO TO 150
571 DO 100 i = 1, n
572 READ( nin, fmt = * )( a( i, j ), j = 1, n )
573 100 CONTINUE
574 DO 110 i = 1, n
575 READ( nin, fmt = * )( b( i, j ), j = 1, n )
576 110 CONTINUE
577 READ( nin, fmt = * )( stru( i ), i = 1, n )
578 READ( nin, fmt = * )( diftru( i ), i = 1, n )
579
580 nptknt = nptknt + 1
581
582
583
584
585
586 CALL slacpy(
'F', n, n, a, lda, ai, lda )
587 CALL slacpy(
'F', n, n, b, lda, bi, lda )
588
589 CALL sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alphar,
590 $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
591 $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
592 $ bwork, linfo )
593
594 IF( linfo.NE.0 ) THEN
595 result( 1 ) = ulpinv
596 WRITE( nout, fmt = 9987 )'SGGEVX', linfo, n, nptknt
597 GO TO 140
598 END IF
599
600
601
602 CALL slacpy(
'Full', n, n, ai, lda, work, n )
603 CALL slacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
604 abnorm =
slange(
'Fro', n, 2*n, work, n, work )
605
606
607
608 result( 1 ) = zero
609 CALL sget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
610 $ beta, work, result( 1 ) )
611 IF( result( 2 ).GT.thresh ) THEN
612 WRITE( nout, fmt = 9986 )'Left', 'SGGEVX', result( 2 ), n,
613 $ nptknt
614 END IF
615
616 result( 2 ) = zero
617 CALL sget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
618 $ beta, work, result( 2 ) )
619 IF( result( 3 ).GT.thresh ) THEN
620 WRITE( nout, fmt = 9986 )'Right', 'SGGEVX', result( 3 ), n,
621 $ nptknt
622 END IF
623
624
625
626 result( 3 ) = zero
627 DO 120 i = 1, n
628 IF( s( i ).EQ.zero ) THEN
629 IF( stru( i ).GT.abnorm*ulp )
630 $ result( 3 ) = ulpinv
631 ELSE IF( stru( i ).EQ.zero ) THEN
632 IF( s( i ).GT.abnorm*ulp )
633 $ result( 3 ) = ulpinv
634 ELSE
635 work( i ) = max( abs( stru( i ) / s( i ) ),
636 $ abs( s( i ) / stru( i ) ) )
637 result( 3 ) = max( result( 3 ), work( i ) )
638 END IF
639 120 CONTINUE
640
641
642
643 result( 4 ) = zero
644 IF( dif( 1 ).EQ.zero ) THEN
645 IF( diftru( 1 ).GT.abnorm*ulp )
646 $ result( 4 ) = ulpinv
647 ELSE IF( diftru( 1 ).EQ.zero ) THEN
648 IF( dif( 1 ).GT.abnorm*ulp )
649 $ result( 4 ) = ulpinv
650 ELSE IF( dif( 5 ).EQ.zero ) THEN
651 IF( diftru( 5 ).GT.abnorm*ulp )
652 $ result( 4 ) = ulpinv
653 ELSE IF( diftru( 5 ).EQ.zero ) THEN
654 IF( dif( 5 ).GT.abnorm*ulp )
655 $ result( 4 ) = ulpinv
656 ELSE
657 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
658 $ abs( dif( 1 ) / diftru( 1 ) ) )
659 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
660 $ abs( dif( 5 ) / diftru( 5 ) ) )
661 result( 4 ) = max( ratio1, ratio2 )
662 END IF
663
664 ntestt = ntestt + 4
665
666
667
668 DO 130 j = 1, 4
669 IF( result( j ).GE.thrsh2 ) THEN
670
671
672
673
674 IF( nerrs.EQ.0 ) THEN
675 WRITE( nout, fmt = 9997 )'SXV'
676
677
678
679
680
681 WRITE( nout, fmt = 9996 )
682
683
684
685 WRITE( nout, fmt = 9992 )'''', 'transpose', ''''
686
687 END IF
688 nerrs = nerrs + 1
689 IF( result( j ).LT.10000.0 ) THEN
690 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
691 ELSE
692 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
693 END IF
694 END IF
695 130 CONTINUE
696
697 140 CONTINUE
698
699 GO TO 90
700 150 CONTINUE
701
702
703
704 CALL alasvm(
'SXV', nout, nerrs, ntestt, 0 )
705
706 work( 1 ) = maxwrk
707
708 RETURN
709
710 9999 FORMAT( ' SDRGVX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
711 $ i6, ', JTYPE=', i6, ')' )
712
713 9998 FORMAT( ' SDRGVX: ', a, ' Eigenvectors from ', a, ' incorrectly ',
714 $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
715 $ 'N=', i6, ', JTYPE=', i6, ', IWA=', i5, ', IWB=', i5,
716 $ ', IWX=', i5, ', IWY=', i5 )
717
718 9997 FORMAT( / 1x, a3, ' -- Real Expert Eigenvalue/vector',
719 $ ' problem driver' )
720
721 9996 FORMAT( ' Input Example' )
722
723 9995 FORMAT( ' Matrix types: ', / )
724
725 9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
726 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
727 $ / ' YH and X are left and right eigenvectors. ', / )
728
729 9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
730 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
731 $ / ' YH and X are left and right eigenvectors. ', / )
732
733 9992 FORMAT( / ' Tests performed: ', / 4x,
734 $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
735 $ ' r is a right eigenvector and ', a, ' means ', a, '.',
736 $ / ' 1 = max | ( b A - a B )', a, ' l | / const.',
737 $ / ' 2 = max | ( b A - a B ) r | / const.',
738 $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
739 $ ' over all eigenvalues', /
740 $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
741 $ ' over the 1st and 5th eigenvectors', / )
742
743 9991 FORMAT( ' Type=', i2, ',', ' IWA=', i2, ', IWB=', i2, ', IWX=',
744 $ i2, ', IWY=', i2, ', result ', i2, ' is', 0p, f8.2 )
745 9990 FORMAT( ' Type=', i2, ',', ' IWA=', i2, ', IWB=', i2, ', IWX=',
746 $ i2, ', IWY=', i2, ', result ', i2, ' is', 1p, e10.3 )
747 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
748 $ ' result ', i2, ' is', 0p, f8.2 )
749 9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
750 $ ' result ', i2, ' is', 1p, e10.3 )
751 9987 FORMAT( ' SDRGVX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
752 $ i6, ', Input example #', i2, ')' )
753
754 9986 FORMAT( ' SDRGVX: ', a, ' Eigenvectors from ', a, ' incorrectly ',
755 $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
756 $ 'N=', i6, ', Input Example #', i2, ')' )
757
758
759
760
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine xerbla(srname, info)
subroutine sggevx(balanc, jobvl, jobvr, sense, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)
SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sget52(left, n, a, lda, b, ldb, e, lde, alphar, alphai, beta, work, result)
SGET52
subroutine slatm6(type, n, a, lda, b, x, ldx, y, ldy, alpha, beta, wx, wy, s, dif)
SLATM6