219
220
221
222
223
224
225 CHARACTER HOWMNY, SIDE
226 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
227
228
229 LOGICAL SELECT( * )
230 DOUBLE PRECISION RWORK( * )
231 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
232 $ VR( LDVR, * ), WORK( * )
233
234
235
236
237
238
239 DOUBLE PRECISION ZERO, ONE
240 parameter( zero = 0.0d+0, one = 1.0d+0 )
241 COMPLEX*16 CZERO, CONE
242 parameter( czero = ( 0.0d+0, 0.0d+0 ),
243 $ cone = ( 1.0d+0, 0.0d+0 ) )
244
245
246 LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
247 $ LSA, LSB
248 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
249 $ J, JE, JR
250 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
251 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
252 $ SCALE, SMALL, TEMP, ULP, XMAX
253 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
254
255
256 LOGICAL LSAME
257 DOUBLE PRECISION DLAMCH
258 COMPLEX*16 ZLADIV
260
261
263
264
265 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
266
267
268 DOUBLE PRECISION ABS1
269
270
271 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
272
273
274
275
276
277 IF(
lsame( howmny,
'A' ) )
THEN
278 ihwmny = 1
279 ilall = .true.
280 ilback = .false.
281 ELSE IF(
lsame( howmny,
'S' ) )
THEN
282 ihwmny = 2
283 ilall = .false.
284 ilback = .false.
285 ELSE IF(
lsame( howmny,
'B' ) )
THEN
286 ihwmny = 3
287 ilall = .true.
288 ilback = .true.
289 ELSE
290 ihwmny = -1
291 END IF
292
293 IF(
lsame( side,
'R' ) )
THEN
294 iside = 1
295 compl = .false.
296 compr = .true.
297 ELSE IF(
lsame( side,
'L' ) )
THEN
298 iside = 2
299 compl = .true.
300 compr = .false.
301 ELSE IF(
lsame( side,
'B' ) )
THEN
302 iside = 3
303 compl = .true.
304 compr = .true.
305 ELSE
306 iside = -1
307 END IF
308
309 info = 0
310 IF( iside.LT.0 ) THEN
311 info = -1
312 ELSE IF( ihwmny.LT.0 ) THEN
313 info = -2
314 ELSE IF( n.LT.0 ) THEN
315 info = -4
316 ELSE IF( lds.LT.max( 1, n ) ) THEN
317 info = -6
318 ELSE IF( ldp.LT.max( 1, n ) ) THEN
319 info = -8
320 END IF
321 IF( info.NE.0 ) THEN
322 CALL xerbla(
'ZTGEVC', -info )
323 RETURN
324 END IF
325
326
327
328 IF( .NOT.ilall ) THEN
329 im = 0
330 DO 10 j = 1, n
331 IF( SELECT( j ) )
332 $ im = im + 1
333 10 CONTINUE
334 ELSE
335 im = n
336 END IF
337
338
339
340 ilbbad = .false.
341 DO 20 j = 1, n
342 IF( dimag( p( j, j ) ).NE.zero )
343 $ ilbbad = .true.
344 20 CONTINUE
345
346 IF( ilbbad ) THEN
347 info = -7
348 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
349 info = -10
350 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
351 info = -12
352 ELSE IF( mm.LT.im ) THEN
353 info = -13
354 END IF
355 IF( info.NE.0 ) THEN
356 CALL xerbla(
'ZTGEVC', -info )
357 RETURN
358 END IF
359
360
361
362 m = im
363 IF( n.EQ.0 )
364 $ RETURN
365
366
367
368 safmin =
dlamch(
'Safe minimum' )
369 big = one / safmin
370 CALL dlabad( safmin, big )
372 small = safmin*n / ulp
373 big = one / small
374 bignum = one / ( safmin*n )
375
376
377
378
379
380 anorm = abs1( s( 1, 1 ) )
381 bnorm = abs1( p( 1, 1 ) )
382 rwork( 1 ) = zero
383 rwork( n+1 ) = zero
384 DO 40 j = 2, n
385 rwork( j ) = zero
386 rwork( n+j ) = zero
387 DO 30 i = 1, j - 1
388 rwork( j ) = rwork( j ) + abs1( s( i, j ) )
389 rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
390 30 CONTINUE
391 anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
392 bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
393 40 CONTINUE
394
395 ascale = one / max( anorm, safmin )
396 bscale = one / max( bnorm, safmin )
397
398
399
400 IF( compl ) THEN
401 ieig = 0
402
403
404
405 DO 140 je = 1, n
406 IF( ilall ) THEN
407 ilcomp = .true.
408 ELSE
409 ilcomp = SELECT( je )
410 END IF
411 IF( ilcomp ) THEN
412 ieig = ieig + 1
413
414 IF( abs1( s( je, je ) ).LE.safmin .AND.
415 $ abs( dble( p( je, je ) ) ).LE.safmin ) THEN
416
417
418
419 DO 50 jr = 1, n
420 vl( jr, ieig ) = czero
421 50 CONTINUE
422 vl( ieig, ieig ) = cone
423 GO TO 140
424 END IF
425
426
427
428
429
430
431 temp = one / max( abs1( s( je, je ) )*ascale,
432 $ abs( dble( p( je, je ) ) )*bscale, safmin )
433 salpha = ( temp*s( je, je ) )*ascale
434 sbeta = ( temp*dble( p( je, je ) ) )*bscale
435 acoeff = sbeta*ascale
436 bcoeff = salpha*bscale
437
438
439
440 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
441 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
442 $ small
443
444 scale = one
445 IF( lsa )
446 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
447 IF( lsb )
448 $ scale = max( scale, ( small / abs1( salpha ) )*
449 $ min( bnorm, big ) )
450 IF( lsa .OR. lsb ) THEN
451 scale = min( scale, one /
452 $ ( safmin*max( one, abs( acoeff ),
453 $ abs1( bcoeff ) ) ) )
454 IF( lsa ) THEN
455 acoeff = ascale*( scale*sbeta )
456 ELSE
457 acoeff = scale*acoeff
458 END IF
459 IF( lsb ) THEN
460 bcoeff = bscale*( scale*salpha )
461 ELSE
462 bcoeff = scale*bcoeff
463 END IF
464 END IF
465
466 acoefa = abs( acoeff )
467 bcoefa = abs1( bcoeff )
468 xmax = one
469 DO 60 jr = 1, n
470 work( jr ) = czero
471 60 CONTINUE
472 work( je ) = cone
473 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
474
475
476
477
478
479
480
481 DO 100 j = je + 1, n
482
483
484
485
486
487
488
489 temp = one / xmax
490 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
491 $ temp ) THEN
492 DO 70 jr = je, j - 1
493 work( jr ) = temp*work( jr )
494 70 CONTINUE
495 xmax = one
496 END IF
497 suma = czero
498 sumb = czero
499
500 DO 80 jr = je, j - 1
501 suma = suma + dconjg( s( jr, j ) )*work( jr )
502 sumb = sumb + dconjg( p( jr, j ) )*work( jr )
503 80 CONTINUE
504 sum = acoeff*suma - dconjg( bcoeff )*sumb
505
506
507
508
509
510 d = dconjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
511 IF( abs1( d ).LE.dmin )
512 $ d = dcmplx( dmin )
513
514 IF( abs1( d ).LT.one ) THEN
515 IF( abs1( sum ).GE.bignum*abs1( d ) ) THEN
516 temp = one / abs1( sum )
517 DO 90 jr = je, j - 1
518 work( jr ) = temp*work( jr )
519 90 CONTINUE
520 xmax = temp*xmax
521 sum = temp*sum
522 END IF
523 END IF
524 work( j ) =
zladiv( -sum, d )
525 xmax = max( xmax, abs1( work( j ) ) )
526 100 CONTINUE
527
528
529
530 IF( ilback ) THEN
531 CALL zgemv(
'N', n, n+1-je, cone, vl( 1, je ), ldvl,
532 $ work( je ), 1, czero, work( n+1 ), 1 )
533 isrc = 2
534 ibeg = 1
535 ELSE
536 isrc = 1
537 ibeg = je
538 END IF
539
540
541
542 xmax = zero
543 DO 110 jr = ibeg, n
544 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
545 110 CONTINUE
546
547 IF( xmax.GT.safmin ) THEN
548 temp = one / xmax
549 DO 120 jr = ibeg, n
550 vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
551 120 CONTINUE
552 ELSE
553 ibeg = n + 1
554 END IF
555
556 DO 130 jr = 1, ibeg - 1
557 vl( jr, ieig ) = czero
558 130 CONTINUE
559
560 END IF
561 140 CONTINUE
562 END IF
563
564
565
566 IF( compr ) THEN
567 ieig = im + 1
568
569
570
571 DO 250 je = n, 1, -1
572 IF( ilall ) THEN
573 ilcomp = .true.
574 ELSE
575 ilcomp = SELECT( je )
576 END IF
577 IF( ilcomp ) THEN
578 ieig = ieig - 1
579
580 IF( abs1( s( je, je ) ).LE.safmin .AND.
581 $ abs( dble( p( je, je ) ) ).LE.safmin ) THEN
582
583
584
585 DO 150 jr = 1, n
586 vr( jr, ieig ) = czero
587 150 CONTINUE
588 vr( ieig, ieig ) = cone
589 GO TO 250
590 END IF
591
592
593
594
595
596
597 temp = one / max( abs1( s( je, je ) )*ascale,
598 $ abs( dble( p( je, je ) ) )*bscale, safmin )
599 salpha = ( temp*s( je, je ) )*ascale
600 sbeta = ( temp*dble( p( je, je ) ) )*bscale
601 acoeff = sbeta*ascale
602 bcoeff = salpha*bscale
603
604
605
606 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
607 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
608 $ small
609
610 scale = one
611 IF( lsa )
612 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
613 IF( lsb )
614 $ scale = max( scale, ( small / abs1( salpha ) )*
615 $ min( bnorm, big ) )
616 IF( lsa .OR. lsb ) THEN
617 scale = min( scale, one /
618 $ ( safmin*max( one, abs( acoeff ),
619 $ abs1( bcoeff ) ) ) )
620 IF( lsa ) THEN
621 acoeff = ascale*( scale*sbeta )
622 ELSE
623 acoeff = scale*acoeff
624 END IF
625 IF( lsb ) THEN
626 bcoeff = bscale*( scale*salpha )
627 ELSE
628 bcoeff = scale*bcoeff
629 END IF
630 END IF
631
632 acoefa = abs( acoeff )
633 bcoefa = abs1( bcoeff )
634 xmax = one
635 DO 160 jr = 1, n
636 work( jr ) = czero
637 160 CONTINUE
638 work( je ) = cone
639 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
640
641
642
643
644
645
646 DO 170 jr = 1, je - 1
647 work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
648 170 CONTINUE
649 work( je ) = cone
650
651 DO 210 j = je - 1, 1, -1
652
653
654
655
656 d = acoeff*s( j, j ) - bcoeff*p( j, j )
657 IF( abs1( d ).LE.dmin )
658 $ d = dcmplx( dmin )
659
660 IF( abs1( d ).LT.one ) THEN
661 IF( abs1( work( j ) ).GE.bignum*abs1( d ) ) THEN
662 temp = one / abs1( work( j ) )
663 DO 180 jr = 1, je
664 work( jr ) = temp*work( jr )
665 180 CONTINUE
666 END IF
667 END IF
668
669 work( j ) =
zladiv( -work( j ), d )
670
671 IF( j.GT.1 ) THEN
672
673
674
675 IF( abs1( work( j ) ).GT.one ) THEN
676 temp = one / abs1( work( j ) )
677 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
678 $ bignum*temp ) THEN
679 DO 190 jr = 1, je
680 work( jr ) = temp*work( jr )
681 190 CONTINUE
682 END IF
683 END IF
684
685 ca = acoeff*work( j )
686 cb = bcoeff*work( j )
687 DO 200 jr = 1, j - 1
688 work( jr ) = work( jr ) + ca*s( jr, j ) -
689 $ cb*p( jr, j )
690 200 CONTINUE
691 END IF
692 210 CONTINUE
693
694
695
696 IF( ilback ) THEN
697 CALL zgemv(
'N', n, je, cone, vr, ldvr, work, 1,
698 $ czero, work( n+1 ), 1 )
699 isrc = 2
700 iend = n
701 ELSE
702 isrc = 1
703 iend = je
704 END IF
705
706
707
708 xmax = zero
709 DO 220 jr = 1, iend
710 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
711 220 CONTINUE
712
713 IF( xmax.GT.safmin ) THEN
714 temp = one / xmax
715 DO 230 jr = 1, iend
716 vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
717 230 CONTINUE
718 ELSE
719 iend = 0
720 END IF
721
722 DO 240 jr = iend + 1, n
723 vr( jr, ieig ) = czero
724 240 CONTINUE
725
726 END IF
727 250 CONTINUE
728 END IF
729
730 RETURN
731
732
733
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
complex *16 function zladiv(X, Y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.