1      subroutine split(a, v, n, l, e1, e2, na, nv)
2c
3c!purpose
4c
5c     given the upper hessenberg matrix a with a 2x2 block
6c     starting at a(l,l), split determines if the
7c     corresponding eigenvalues are real or complex, if they
8c     are real, a rotation is determined that reduces the
9c     block to upper triangular form with the eigenvalue
10c     of largest absolute value appearing first.  the
11c     rotation is accumulated in v. the eigenvalues (real
12c     or complex) are returned in e1 and e2.
13c!calling sequence
14c
15c     subroutine split(a, v, n, l, e1, e2, na, nv)
16c
17c     double precision a,v,e1,e2
18c     integer n,l,na,nv
19c     dimension a(na,n),v(nv,n)
20c
21c     starred parameters are  altered by the subroutine
22c
23c        *a        the upper hessenberg matrix whose 2x2
24c                  block is to be dsplit.
25c        *v        the array in which the dsplitting trans-
26c                  formation is to be accumulated.
27c         n        the order of the matrix a.
28c         l        the position of the 2x2 block.
29c        *e1       on return if the eigenvalues are complex
30c        *e2       e1 contains their common real part and
31c                  e2 contains the positive imaginary part.
32c                  if the eigenvalues are real. e1 contains
33c                  the one largest in absolute value and f2
34c                  contains the other one.
35c        na        the first dimension of the array a.
36c        nv        the first dimension of the array v.
37c!auxiliary routines
38c     abs sqrt (fortran)
39c!
40c originator
41c
42      double precision a,v,e1,e2
43      integer n,l,na,nv
44      dimension a(na,n),v(nv,n)
45c     internal variables
46c
47c     internal variables
48      double precision p,q,r,t,u,w,x,y,z,zero,two
49      integer i,j,l1
50      data zero, two /0.0d+0,2.0d+0/
51      l1 = l + 1
52c
53      x = a(l1,l1)
54      y = a(l,l)
55      w = a(l,l1)*a(l1,l)
56      p = (y-x)/two
57      q = p**2 + w
58      if (q.ge.zero) go to 10
59c
60c       complex eigenvalue.
61c
62      e1 = p + x
63      e2 = sqrt(-q)
64      return
65   10 continue
66c
67c     two real eigenvalues. set up transformation.
68c
69      z = sqrt(q)
70      if (p.lt.zero) go to 20
71      z = p + z
72      go to 30
73   20 continue
74      z = p - z
75   30 continue
76      if (z.eq.zero) go to 40
77      r = -w/z
78      go to 50
79   40 continue
80      r = zero
81   50 continue
82      if (abs(x+z).ge.abs(x+r)) z = r
83      y = y - x - z
84      x = -z
85      t = a(l,l1)
86      u = a(l1,l)
87      if (abs(y)+abs(u).le.abs(t)+abs(x)) go to 60
88      q = u
89      p = y
90      go to 70
91   60 continue
92      q = x
93      p = t
94   70 continue
95      r = sqrt(p**2+q**2)
96      if (r.gt.zero) go to 80
97      e1 = a(l,l)
98      e2 = a(l1,l1)
99      a(l1,l) = zero
100      return
101   80 continue
102      p = p/r
103      q = q/r
104c
105c     premultiply.
106c
107      do 90 j=l,n
108         z = a(l,j)
109         a(l,j) = p*z + q*a(l1,j)
110         a(l1,j) = p*a(l1,j) - q*z
111   90 continue
112c
113c     postmultiply.
114c
115      do 100 i=1,l1
116         z = a(i,l)
117         a(i,l) = p*z + q*a(i,l1)
118         a(i,l1) = p*a(i,l1) - q*z
119  100 continue
120c
121c     accumulate the transformation in v.
122c
123      do 110 i=1,n
124         z = v(i,l)
125         v(i,l) = p*z + q*v(i,l1)
126         v(i,l1) = p*v(i,l1) - q*z
127  110 continue
128      a(l1,l) = zero
129      e1 = a(l,l)
130      e2 = a(l1,l1)
131      return
132      end
133
134