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