@process directive('" (') subroutine cinvit(nm,n,ar,ai,wr,wi,select,mm,m,zr,zi, x ierr,rm1,rm2,rv1,rv2) c integer i,j,k,m,n,s,mm,mp,nm,uk,its,km1,ierr real ar(nm,n),ai(nm,n),wr(n),wi(n),zr(nm,mm), x zi(nm,mm),rm1(n,n),rm2(n,n),rv1(n),rv2(n) real x,y,eps3,norm,normv,epslon,growto,ilambd,pythag, x rlambd,ukroot,tmp,foo logical select(n) c c this subroutine is a translation of the algol procedure cx invit c by peters and wilkinson. c handbook for auto. comp. vol.ii-linear algebra, 418-439(1971). c c this subroutine finds those eigenvectors of a complex upper c hessenberg matrix corresponding to specified eigenvalues, c using inverse iteration. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain the real and imaginary parts, c respectively, of the hessenberg matrix. c c wr and wi contain the real and imaginary parts, respectively, c of the eigenvalues of the matrix. the eigenvalues must be c stored in a manner identical to that of subroutine comlr, c which recognizes possible splitting of the matrix. c c select specifies the eigenvectors to be found. the c eigenvector corresponding to the j-th eigenvalue is c specified by setting select(j) to .true.. c c mm should be set to an upper bound for the number of c eigenvectors to be found. c c on output c c ar, ai, wi, and select are unaltered. c c wr may have been altered since close eigenvalues are perturbed c slightly in searching for independent eigenvectors. c c m is the number of eigenvectors actually found. c c zr and zi contain the real and imaginary parts, respectively, c of the eigenvectors. the eigenvectors are normalized c so that the component of largest magnitude is 1. c any vector which fails the acceptance test is set to zero. c c ierr is set to c zero for normal return, c -(2*n+1) if more than mm eigenvectors have been specified, c -k if the iteration corresponding to the k-th c value fails, c -(n+k) if both error situations occur. c c rm1, rm2, rv1, and rv2 are temporary storage arrays. c c the algol procedure guessvec appears in cinvit in line. c c calls cdiv for complex division. c calls pythag for dsqrt(a*a + b*b) . c c Questions and comments should be directed to Alan K. Cline, c Pleasant Valley Software, 8603 Altus Cove, Austin, TX 78759. c Electronic mail to cline@cs.utexas.edu. c c this version dated january 1989. (for the IBM 3090vf) c c ------------------------------------------------------------------ c call xuflow(0) ierr = 0 uk = 0 s = 1 c do 980 k = 1, n if (.not. select(k)) go to 980 if (s .gt. mm) go to 1000 if (uk .ge. k) go to 200 c .......... check for possible splitting .......... do 120 uk = k, n if (uk .eq. n) go to 140 if (ar(uk+1,uk) .eq. 0.0e0 .and. ai(uk+1,uk) .eq. 0.0e0) x go to 140 120 continue c .......... compute infinity norm of leading uk by uk c (hessenberg) matrix .......... 140 norm = 0.0e0 do 180 i = 1, uk x = 0.0e0 c do 160 j = max0(1,i-1), uk c-- tmp = pythag(ar(i,j),ai(i,j)) if (abs(ar(i,j)).gt.abs(ai(i,j))) then tmp = abs(ar(i,j))*sqrt(1e0+(ai(i,j)/ar(i,j))**2) else if (ai(i,j).ne.0e0) then tmp = abs(ai(i,j))*sqrt((ar(i,j)/ai(i,j))**2+1e0) else tmp = abs(ar(i,j)) endif x = x + tmp 160 continue c norm = amax1(norm,x) 180 continue c .......... eps3 replaces zero pivot in decomposition c and close roots are modified by eps3 .......... if (norm .eq. 0.0e0) norm = 1.0e0 eps3 = epslon(norm) c .......... growto is the criterion for growth .......... ukroot = uk ukroot = sqrt(ukroot) growto = 0.1e0 / ukroot 200 rlambd = wr(k) ilambd = wi(k) if (k .eq. 1) go to 280 km1 = k - 1 go to 240 c .......... perturb eigenvalue if it is close c to any previous eigenvalue .......... 220 rlambd = rlambd + eps3 c .......... for i=k-1 step -1 until 1 do -- .......... 240 do 260 i = k-1,1,-1 if (select(i) .and. abs(wr(i)-rlambd) .lt. eps3 .and. x abs(wi(i)-ilambd) .lt. eps3) go to 220 260 continue c wr(k) = rlambd c .......... form upper hessenberg (ar,ai)-(rlambd,ilambd)*i c and initial complex vector .......... 280 do 320 i = 1, uk do 300 j = max0(1,i-1), uk rm1(i,j) = ar(i,j) rm2(i,j) = ai(i,j) 300 continue c rm1(i,i) = rm1(i,i) - rlambd rm2(i,i) = rm2(i,i) - ilambd rv1(i) = eps3 320 continue c .......... triangular decomposition with interchanges, c replacing zero pivots by eps3 .......... if (uk .eq. 1) go to 420 c do 400 i = 2, uk mp = i - 1 if (pythag(rm1(i,mp),rm2(i,mp)) .le. x pythag(rm1(mp,mp),rm2(mp,mp))) go to 360 c c" ( prefer vector do 340 j = mp, uk foo = rm1(i,j) rm1(i,j) = rm1(mp,j) rm1(mp,j) = foo foo = rm2(i,j) rm2(i,j) = rm2(mp,j) rm2(mp,j) = foo 340 continue c 360 if (rm1(mp,mp) .eq. 0.0e0 .and. rm2(mp,mp) .eq. 0.0e0) x rm1(mp,mp) = eps3 call cdiv(rm1(i,mp),rm2(i,mp),rm1(mp,mp),rm2(mp,mp),x,y) if (x .eq. 0.0e0 .and. y .eq. 0.0e0) go to 400 c x = -x c" ( prefer vector do 380 j = i, uk rm1(i,j) = rm1(i,j) + x * rm1(mp,j) + y * rm2(mp,j) rm2(i,j) = rm2(i,j) + x * rm2(mp,j) - y * rm1(mp,j) 380 continue c 400 continue c 420 if (rm1(uk,uk) .eq. 0.0e0 .and. rm2(uk,uk) .eq. 0.0e0) x rm1(uk,uk) = eps3 its = 0 c .......... back substitution c for i=uk step -1 until 1 do -- .......... 660 do 720 i = uk, 1, -1 x = rv1(i) y = 0.0e0 c do 680 j = i+1, uk x = x - rm1(i,j) * rv1(j) + rm2(i,j) * rv2(j) y = y - rm1(i,j) * rv2(j) - rm2(i,j) * rv1(j) 680 continue c call cdiv(x,y,rm1(i,i),rm2(i,i),rv1(i),rv2(i)) 720 continue c .......... acceptance test for eigenvector c and normalization .......... its = its + 1 norm = 0.0e0 normv = 0.0e0 c do 780 i = 1, uk c-- foo = pythag(rv1(i),rv2(i)) if (abs(rv1(i)).gt.abs(rv2(i))) then foo = abs(rv1(i))*sqrt(1e0+(rv2(i)/rv1(i))**2) else if (rv2(i).ne.0e0) then foo = abs(rv2(i))*sqrt((rv1(i)/rv2(i))**2+1e0) else foo = abs(rv1(i)) endif if (normv .ge. foo) go to 760 normv = foo j = i 760 norm = norm + foo 780 continue c if (norm .lt. growto) go to 840 c .......... accept vector .......... x = rv1(j) y = rv2(j) c do 820 i = 1, uk call cdiv(rv1(i),rv2(i),x,y,zr(i,s),zi(i,s)) 820 continue c if (uk .eq. n) go to 940 j = uk + 1 go to 900 c .......... in-line procedure for choosing c a new starting vector .......... 840 if (its .ge. uk) go to 880 x = ukroot y = eps3 / (x + 1.0e0) rv1(1) = eps3 c do 860 i = 2, uk 860 rv1(i) = y c j = uk - its + 1 rv1(j) = rv1(j) - eps3 * x go to 660 c .......... set error -- unaccepted eigenvector .......... 880 j = 1 ierr = -k c .......... set remaining vector components to zero .......... 900 do 920 i = j, n zr(i,s) = 0.0e0 zi(i,s) = 0.0e0 920 continue c 940 s = s + 1 980 continue c go to 1001 c .......... set error -- underestimate of eigenvector c space required .......... 1000 if (ierr .ne. 0) ierr = ierr - n if (ierr .eq. 0) ierr = -(2 * n + 1) 1001 m = s - 1 return end