1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
     89
     90
     91
     92
     93
     94
     95
     96
     97
      SUBROUTINE CLAQR1( N, H, LDH, S1, S2, V )
*
*  -- LAPACK auxiliary routine (version 3.2) --
*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      COMPLEX            S1, S2
      INTEGER            LDH, N
*     ..
*     .. Array Arguments ..
      COMPLEX            H( LDH, * ), V( * )
*     ..
*
*       Given a 2-by-2 or 3-by-3 matrix H, CLAQR1 sets v to a
*       scalar multiple of the first column of the product
*
*       (*)  K = (H - s1*I)*(H - s2*I)
*
*       scaling to avoid overflows and most underflows.
*
*       This is useful for starting double implicit shift bulges
*       in the QR algorithm.
*
*
*       N      (input) integer
*              Order of the matrix H. N must be either 2 or 3.
*
*       H      (input) COMPLEX array of dimension (LDH,N)
*              The 2-by-2 or 3-by-3 matrix H in (*).
*
*       LDH    (input) integer
*              The leading dimension of H as declared in
*              the calling procedure.  LDH.GE.N
*
*       S1     (input) COMPLEX
*       S2     S1 and S2 are the shifts defining K in (*) above.
*
*       V      (output) COMPLEX array of dimension N
*              A scalar multiple of the first column of the
*              matrix K in (*).
*
*     ================================================================
*     Based on contributions by
*        Karen Braman and Ralph Byers, Department of Mathematics,
*        University of Kansas, USA
*
*     ================================================================
*
*     .. Parameters ..
      COMPLEX            ZERO
      PARAMETER          ( ZERO = ( 0.0e00.0e0 ) )
      REAL               RZERO
      PARAMETER          ( RZERO = 0.0e0 )
*     ..
*     .. Local Scalars ..
      COMPLEX            CDUM, H21S, H31S
      REAL               S
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABSAIMAG, REAL
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      CABS1( CDUM ) = ABSREAL( CDUM ) ) + ABSAIMAG( CDUM ) )
*     ..
*     .. Executable Statements ..
      IF( N.EQ.2 ) THEN
         S = CABS1( H( 11 )-S2 ) + CABS1( H( 21 ) )
         IF( S.EQ.RZERO ) THEN
            V( 1 ) = ZERO
            V( 2 ) = ZERO
         ELSE
            H21S = H( 21 ) / S
            V( 1 ) = H21S*H( 12 ) + ( H( 11 )-S1 )*
     $               ( ( H( 11 )-S2 ) / S )
            V( 2 ) = H21S*( H( 11 )+H( 22 )-S1-S2 )
         END IF
      ELSE
         S = CABS1( H( 11 )-S2 ) + CABS1( H( 21 ) ) +
     $       CABS1( H( 31 ) )
         IF( S.EQ.ZERO ) THEN
            V( 1 ) = ZERO
            V( 2 ) = ZERO
            V( 3 ) = ZERO
         ELSE
            H21S = H( 21 ) / S
            H31S = H( 31 ) / S
            V( 1 ) = ( H( 11 )-S1 )*( ( H( 11 )-S2 ) / S ) +
     $               H( 12 )*H21S + H( 13 )*H31S
            V( 2 ) = H21S*( H( 11 )+H( 22 )-S1-S2 ) + H( 23 )*H31S
            V( 3 ) = H31S*( H( 11 )+H( 33 )-S1-S2 ) + H21S*H( 32 )
         END IF
      END IF
      END