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