1 subroutine ortbak(nm,low,igh,a,ort,m,z) 2c 3 integer i,j,m,la,mm,mp,nm,igh,kp1,low,mp1 4 double precision a(nm,igh),ort(igh),z(nm,m) 5 double precision g 6c 7c this subroutine is a translation of the algol procedure ortbak, 8c num. math. 12, 349-368(1968) by martin and wilkinson. 9c handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). 10c 11c this subroutine forms the eigenvectors of a real general 12c matrix by back transforming those of the corresponding 13c upper hessenberg matrix determined by orthes. 14c 15c on input 16c 17c nm must be set to the row dimension of two-dimensional 18c array parameters as declared in the calling program 19c dimension statement. 20c 21c low and igh are integers determined by the balancing 22c subroutine balanc. if balanc has not been used, 23c set low=1 and igh equal to the order of the matrix. 24c 25c a contains information about the orthogonal trans- 26c formations used in the reduction by orthes 27c in its strict lower triangle. 28c 29c ort contains further information about the trans- 30c formations used in the reduction by orthes. 31c only elements low through igh are used. 32c 33c m is the number of columns of z to be back transformed. 34c 35c z contains the real and imaginary parts of the eigen- 36c vectors to be back transformed in its first m columns. 37c 38c on output 39c 40c z contains the real and imaginary parts of the 41c transformed eigenvectors in its first m columns. 42c 43c ort has been altered. 44c 45c note that ortbak preserves vector euclidean norms. 46c 47c questions and comments should be directed to burton s. garbow, 48c mathematics and computer science div, argonne national laboratory 49c 50c this version dated august 1983. 51c 52c ------------------------------------------------------------------ 53c 54 if (m .eq. 0) go to 200 55 la = igh - 1 56 kp1 = low + 1 57 if (la .lt. kp1) go to 200 58c .......... for mp=igh-1 step -1 until low+1 do -- .......... 59 do 140 mm = kp1, la 60 mp = low + igh - mm 61 if (a(mp,mp-1) .eq. 0.0d0) go to 140 62 mp1 = mp + 1 63c 64 do 100 i = mp1, igh 65 100 ort(i) = a(i,mp-1) 66c 67 do 130 j = 1, m 68 g = 0.0d0 69c 70 do 110 i = mp, igh 71 110 g = g + ort(i) * z(i,j) 72c .......... divisor below is negative of h formed in orthes. 73c double division avoids possible underflow .......... 74 g = (g / ort(mp)) / a(mp,mp-1) 75c 76 do 120 i = mp, igh 77 120 z(i,j) = z(i,j) + g * ort(i) 78c 79 130 continue 80c 81 140 continue 82c 83 200 return 84 end 85