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 REAL RWORK( * )
231 COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
232 $ VR( LDVR, * ), WORK( * )
233
234
235
236
237
238
239 REAL ZERO, ONE
240 parameter( zero = 0.0e+0, one = 1.0e+0 )
241 COMPLEX CZERO, CONE
242 parameter( czero = ( 0.0e+0, 0.0e+0 ),
243 $ cone = ( 1.0e+0, 0.0e+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 REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
251 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
252 $ SCALE, SMALL, TEMP, ULP, XMAX
253 COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
254
255
256 LOGICAL LSAME
257 REAL SLAMCH
258 COMPLEX CLADIV
260
261
263
264
265 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
266
267
268 REAL ABS1
269
270
271 abs1( x ) = abs( real( x ) ) + abs( aimag( 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(
'CTGEVC', -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( aimag( 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(
'CTGEVC', -info )
357 RETURN
358 END IF
359
360
361
362 m = im
363 IF( n.EQ.0 )
364 $ RETURN
365
366
367
368 safmin =
slamch(
'Safe minimum' )
369 big = one / safmin
371 small = safmin*n / ulp
372 big = one / small
373 bignum = one / ( safmin*n )
374
375
376
377
378
379 anorm = abs1( s( 1, 1 ) )
380 bnorm = abs1( p( 1, 1 ) )
381 rwork( 1 ) = zero
382 rwork( n+1 ) = zero
383 DO 40 j = 2, n
384 rwork( j ) = zero
385 rwork( n+j ) = zero
386 DO 30 i = 1, j - 1
387 rwork( j ) = rwork( j ) + abs1( s( i, j ) )
388 rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
389 30 CONTINUE
390 anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
391 bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
392 40 CONTINUE
393
394 ascale = one / max( anorm, safmin )
395 bscale = one / max( bnorm, safmin )
396
397
398
399 IF( compl ) THEN
400 ieig = 0
401
402
403
404 DO 140 je = 1, n
405 IF( ilall ) THEN
406 ilcomp = .true.
407 ELSE
408 ilcomp = SELECT( je )
409 END IF
410 IF( ilcomp ) THEN
411 ieig = ieig + 1
412
413 IF( abs1( s( je, je ) ).LE.safmin .AND.
414 $ abs( real( p( je, je ) ) ).LE.safmin ) THEN
415
416
417
418 DO 50 jr = 1, n
419 vl( jr, ieig ) = czero
420 50 CONTINUE
421 vl( ieig, ieig ) = cone
422 GO TO 140
423 END IF
424
425
426
427
428
429
430 temp = one / max( abs1( s( je, je ) )*ascale,
431 $ abs( real( p( je, je ) ) )*bscale, safmin )
432 salpha = ( temp*s( je, je ) )*ascale
433 sbeta = ( temp*real( p( je, je ) ) )*bscale
434 acoeff = sbeta*ascale
435 bcoeff = salpha*bscale
436
437
438
439 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
440 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
441 $ small
442
443 scale = one
444 IF( lsa )
445 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
446 IF( lsb )
447 $ scale = max( scale, ( small / abs1( salpha ) )*
448 $ min( bnorm, big ) )
449 IF( lsa .OR. lsb ) THEN
450 scale = min( scale, one /
451 $ ( safmin*max( one, abs( acoeff ),
452 $ abs1( bcoeff ) ) ) )
453 IF( lsa ) THEN
454 acoeff = ascale*( scale*sbeta )
455 ELSE
456 acoeff = scale*acoeff
457 END IF
458 IF( lsb ) THEN
459 bcoeff = bscale*( scale*salpha )
460 ELSE
461 bcoeff = scale*bcoeff
462 END IF
463 END IF
464
465 acoefa = abs( acoeff )
466 bcoefa = abs1( bcoeff )
467 xmax = one
468 DO 60 jr = 1, n
469 work( jr ) = czero
470 60 CONTINUE
471 work( je ) = cone
472 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
473
474
475
476
477
478
479
480 DO 100 j = je + 1, n
481
482
483
484
485
486
487
488 temp = one / xmax
489 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
490 $ temp ) THEN
491 DO 70 jr = je, j - 1
492 work( jr ) = temp*work( jr )
493 70 CONTINUE
494 xmax = one
495 END IF
496 suma = czero
497 sumb = czero
498
499 DO 80 jr = je, j - 1
500 suma = suma + conjg( s( jr, j ) )*work( jr )
501 sumb = sumb + conjg( p( jr, j ) )*work( jr )
502 80 CONTINUE
503 sum = acoeff*suma - conjg( bcoeff )*sumb
504
505
506
507
508
509 d = conjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
510 IF( abs1( d ).LE.dmin )
511 $ d = cmplx( dmin )
512
513 IF( abs1( d ).LT.one ) THEN
514 IF( abs1( sum ).GE.bignum*abs1( d ) ) THEN
515 temp = one / abs1( sum )
516 DO 90 jr = je, j - 1
517 work( jr ) = temp*work( jr )
518 90 CONTINUE
519 xmax = temp*xmax
520 sum = temp*sum
521 END IF
522 END IF
523 work( j ) =
cladiv( -sum, d )
524 xmax = max( xmax, abs1( work( j ) ) )
525 100 CONTINUE
526
527
528
529 IF( ilback ) THEN
530 CALL cgemv(
'N', n, n+1-je, cone, vl( 1, je ), ldvl,
531 $ work( je ), 1, czero, work( n+1 ), 1 )
532 isrc = 2
533 ibeg = 1
534 ELSE
535 isrc = 1
536 ibeg = je
537 END IF
538
539
540
541 xmax = zero
542 DO 110 jr = ibeg, n
543 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
544 110 CONTINUE
545
546 IF( xmax.GT.safmin ) THEN
547 temp = one / xmax
548 DO 120 jr = ibeg, n
549 vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
550 120 CONTINUE
551 ELSE
552 ibeg = n + 1
553 END IF
554
555 DO 130 jr = 1, ibeg - 1
556 vl( jr, ieig ) = czero
557 130 CONTINUE
558
559 END IF
560 140 CONTINUE
561 END IF
562
563
564
565 IF( compr ) THEN
566 ieig = im + 1
567
568
569
570 DO 250 je = n, 1, -1
571 IF( ilall ) THEN
572 ilcomp = .true.
573 ELSE
574 ilcomp = SELECT( je )
575 END IF
576 IF( ilcomp ) THEN
577 ieig = ieig - 1
578
579 IF( abs1( s( je, je ) ).LE.safmin .AND.
580 $ abs( real( p( je, je ) ) ).LE.safmin ) THEN
581
582
583
584 DO 150 jr = 1, n
585 vr( jr, ieig ) = czero
586 150 CONTINUE
587 vr( ieig, ieig ) = cone
588 GO TO 250
589 END IF
590
591
592
593
594
595
596 temp = one / max( abs1( s( je, je ) )*ascale,
597 $ abs( real( p( je, je ) ) )*bscale, safmin )
598 salpha = ( temp*s( je, je ) )*ascale
599 sbeta = ( temp*real( p( je, je ) ) )*bscale
600 acoeff = sbeta*ascale
601 bcoeff = salpha*bscale
602
603
604
605 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
606 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
607 $ small
608
609 scale = one
610 IF( lsa )
611 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
612 IF( lsb )
613 $ scale = max( scale, ( small / abs1( salpha ) )*
614 $ min( bnorm, big ) )
615 IF( lsa .OR. lsb ) THEN
616 scale = min( scale, one /
617 $ ( safmin*max( one, abs( acoeff ),
618 $ abs1( bcoeff ) ) ) )
619 IF( lsa ) THEN
620 acoeff = ascale*( scale*sbeta )
621 ELSE
622 acoeff = scale*acoeff
623 END IF
624 IF( lsb ) THEN
625 bcoeff = bscale*( scale*salpha )
626 ELSE
627 bcoeff = scale*bcoeff
628 END IF
629 END IF
630
631 acoefa = abs( acoeff )
632 bcoefa = abs1( bcoeff )
633 xmax = one
634 DO 160 jr = 1, n
635 work( jr ) = czero
636 160 CONTINUE
637 work( je ) = cone
638 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
639
640
641
642
643
644
645 DO 170 jr = 1, je - 1
646 work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
647 170 CONTINUE
648 work( je ) = cone
649
650 DO 210 j = je - 1, 1, -1
651
652
653
654
655 d = acoeff*s( j, j ) - bcoeff*p( j, j )
656 IF( abs1( d ).LE.dmin )
657 $ d = cmplx( dmin )
658
659 IF( abs1( d ).LT.one ) THEN
660 IF( abs1( work( j ) ).GE.bignum*abs1( d ) ) THEN
661 temp = one / abs1( work( j ) )
662 DO 180 jr = 1, je
663 work( jr ) = temp*work( jr )
664 180 CONTINUE
665 END IF
666 END IF
667
668 work( j ) =
cladiv( -work( j ), d )
669
670 IF( j.GT.1 ) THEN
671
672
673
674 IF( abs1( work( j ) ).GT.one ) THEN
675 temp = one / abs1( work( j ) )
676 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
677 $ bignum*temp ) THEN
678 DO 190 jr = 1, je
679 work( jr ) = temp*work( jr )
680 190 CONTINUE
681 END IF
682 END IF
683
684 ca = acoeff*work( j )
685 cb = bcoeff*work( j )
686 DO 200 jr = 1, j - 1
687 work( jr ) = work( jr ) + ca*s( jr, j ) -
688 $ cb*p( jr, j )
689 200 CONTINUE
690 END IF
691 210 CONTINUE
692
693
694
695 IF( ilback ) THEN
696 CALL cgemv(
'N', n, je, cone, vr, ldvr, work, 1,
697 $ czero, work( n+1 ), 1 )
698 isrc = 2
699 iend = n
700 ELSE
701 isrc = 1
702 iend = je
703 END IF
704
705
706
707 xmax = zero
708 DO 220 jr = 1, iend
709 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
710 220 CONTINUE
711
712 IF( xmax.GT.safmin ) THEN
713 temp = one / xmax
714 DO 230 jr = 1, iend
715 vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
716 230 CONTINUE
717 ELSE
718 iend = 0
719 END IF
720
721 DO 240 jr = iend + 1, n
722 vr( jr, ieig ) = czero
723 240 CONTINUE
724
725 END IF
726 250 CONTINUE
727 END IF
728
729 RETURN
730
731
732
subroutine xerbla(srname, info)
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
complex function cladiv(x, y)
CLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
real function slamch(cmach)
SLAMCH
logical function lsame(ca, cb)
LSAME