home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / cinvit.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  253 lines
  1.       subroutine cinvit(nm,n,ar,ai,wr,wi,select,mm,m,zr,zi,
  2.      x                  ierr,rm1,rm2,rv1,rv2)
  3. c
  4.       integer i,j,k,m,n,s,ii,mm,mp,nm,uk,ip1,its,km1,ierr
  5.       double precision ar(nm,n),ai(nm,n),wr(n),wi(n),zr(nm,mm),
  6.      x       zi(nm,mm),rm1(n,n),rm2(n,n),rv1(n),rv2(n)
  7.       double precision x,y,eps3,norm,normv,epslon,growto,ilambd,pythag,
  8.      x       rlambd,ukroot
  9.       logical select(n)
  10. c
  11. c     this subroutine is a translation of the algol procedure cx invit
  12. c     by peters and wilkinson.
  13. c     handbook for auto. comp. vol.ii-linear algebra, 418-439(1971).
  14. c
  15. c     this subroutine finds those eigenvectors of a complex upper
  16. c     hessenberg matrix corresponding to specified eigenvalues,
  17. c     using inverse iteration.
  18. c
  19. c     on input
  20. c
  21. c        nm must be set to the row dimension of two-dimensional
  22. c          array parameters as declared in the calling program
  23. c          dimension statement.
  24. c
  25. c        n is the order of the matrix.
  26. c
  27. c        ar and ai contain the real and imaginary parts,
  28. c          respectively, of the hessenberg matrix.
  29. c
  30. c        wr and wi contain the real and imaginary parts, respectively,
  31. c          of the eigenvalues of the matrix.  the eigenvalues must be
  32. c          stored in a manner identical to that of subroutine  comlr,
  33. c          which recognizes possible splitting of the matrix.
  34. c
  35. c        select specifies the eigenvectors to be found.  the
  36. c          eigenvector corresponding to the j-th eigenvalue is
  37. c          specified by setting select(j) to .true..
  38. c
  39. c        mm should be set to an upper bound for the number of
  40. c          eigenvectors to be found.
  41. c
  42. c     on output
  43. c
  44. c        ar, ai, wi, and select are unaltered.
  45. c
  46. c        wr may have been altered since close eigenvalues are perturbed
  47. c          slightly in searching for independent eigenvectors.
  48. c
  49. c        m is the number of eigenvectors actually found.
  50. c
  51. c        zr and zi contain the real and imaginary parts, respectively,
  52. c          of the eigenvectors.  the eigenvectors are normalized
  53. c          so that the component of largest magnitude is 1.
  54. c          any vector which fails the acceptance test is set to zero.
  55. c
  56. c        ierr is set to
  57. c          zero       for normal return,
  58. c          -(2*n+1)   if more than mm eigenvectors have been specified,
  59. c          -k         if the iteration corresponding to the k-th
  60. c                     value fails,
  61. c          -(n+k)     if both error situations occur.
  62. c
  63. c        rm1, rm2, rv1, and rv2 are temporary storage arrays.
  64. c
  65. c     the algol procedure guessvec appears in cinvit in line.
  66. c
  67. c     calls cdiv for complex division.
  68. c     calls pythag for  dsqrt(a*a + b*b) .
  69. c
  70. c     questions and comments should be directed to burton s. garbow,
  71. c     mathematics and computer science div, argonne national laboratory
  72. c
  73. c     this version dated august 1983.
  74. c
  75. c     ------------------------------------------------------------------
  76. c
  77.       ierr = 0
  78.       uk = 0
  79.       s = 1
  80. c
  81.       do 980 k = 1, n
  82.          if (.not. select(k)) go to 980
  83.          if (s .gt. mm) go to 1000
  84.          if (uk .ge. k) go to 200
  85. c     .......... check for possible splitting ..........
  86.          do 120 uk = k, n
  87.             if (uk .eq. n) go to 140
  88.             if (ar(uk+1,uk) .eq. 0.0d0 .and. ai(uk+1,uk) .eq. 0.0d0)
  89.      x         go to 140
  90.   120    continue
  91. c     .......... compute infinity norm of leading uk by uk
  92. c                (hessenberg) matrix ..........
  93.   140    norm = 0.0d0
  94.          mp = 1
  95. c
  96.          do 180 i = 1, uk
  97.             x = 0.0d0
  98. c
  99.             do 160 j = mp, uk
  100.   160       x = x + pythag(ar(i,j),ai(i,j))
  101. c
  102.             if (x .gt. norm) norm = x
  103.             mp = i
  104.   180    continue
  105. c     .......... eps3 replaces zero pivot in decomposition
  106. c                and close roots are modified by eps3 ..........
  107.          if (norm .eq. 0.0d0) norm = 1.0d0
  108.          eps3 = epslon(norm)
  109. c     .......... growto is the criterion for growth ..........
  110.          ukroot = uk
  111.          ukroot = dsqrt(ukroot)
  112.          growto = 0.1d0 / ukroot
  113.   200    rlambd = wr(k)
  114.          ilambd = wi(k)
  115.          if (k .eq. 1) go to 280
  116.          km1 = k - 1
  117.          go to 240
  118. c     .......... perturb eigenvalue if it is close
  119. c                to any previous eigenvalue ..........
  120.   220    rlambd = rlambd + eps3
  121. c     .......... for i=k-1 step -1 until 1 do -- ..........
  122.   240    do 260 ii = 1, km1
  123.             i = k - ii
  124.             if (select(i) .and. dabs(wr(i)-rlambd) .lt. eps3 .and.
  125.      x         dabs(wi(i)-ilambd) .lt. eps3) go to 220
  126.   260    continue
  127. c
  128.          wr(k) = rlambd
  129. c     .......... form upper hessenberg (ar,ai)-(rlambd,ilambd)*i
  130. c                and initial complex vector ..........
  131.   280    mp = 1
  132. c
  133.          do 320 i = 1, uk
  134. c
  135.             do 300 j = mp, uk
  136.                rm1(i,j) = ar(i,j)
  137.                rm2(i,j) = ai(i,j)
  138.   300       continue
  139. c
  140.             rm1(i,i) = rm1(i,i) - rlambd
  141.             rm2(i,i) = rm2(i,i) - ilambd
  142.             mp = i
  143.             rv1(i) = eps3
  144.   320    continue
  145. c     .......... triangular decomposition with interchanges,
  146. c                replacing zero pivots by eps3 ..........
  147.          if (uk .eq. 1) go to 420
  148. c
  149.          do 400 i = 2, uk
  150.             mp = i - 1
  151.             if (pythag(rm1(i,mp),rm2(i,mp)) .le.
  152.      x          pythag(rm1(mp,mp),rm2(mp,mp))) go to 360
  153. c
  154.             do 340 j = mp, uk
  155.                y = rm1(i,j)
  156.                rm1(i,j) = rm1(mp,j)
  157.                rm1(mp,j) = y
  158.                y = rm2(i,j)
  159.                rm2(i,j) = rm2(mp,j)
  160.                rm2(mp,j) = y
  161.   340       continue
  162. c
  163.   360       if (rm1(mp,mp) .eq. 0.0d0 .and. rm2(mp,mp) .eq. 0.0d0)
  164.      x         rm1(mp,mp) = eps3
  165.             call cdiv(rm1(i,mp),rm2(i,mp),rm1(mp,mp),rm2(mp,mp),x,y)
  166.             if (x .eq. 0.0d0 .and. y .eq. 0.0d0) go to 400
  167. c
  168.             do 380 j = i, uk
  169.                rm1(i,j) = rm1(i,j) - x * rm1(mp,j) + y * rm2(mp,j)
  170.                rm2(i,j) = rm2(i,j) - x * rm2(mp,j) - y * rm1(mp,j)
  171.   380       continue
  172. c
  173.   400    continue
  174. c
  175.   420    if (rm1(uk,uk) .eq. 0.0d0 .and. rm2(uk,uk) .eq. 0.0d0)
  176.      x      rm1(uk,uk) = eps3
  177.          its = 0
  178. c     .......... back substitution
  179. c                for i=uk step -1 until 1 do -- ..........
  180.   660    do 720 ii = 1, uk
  181.             i = uk + 1 - ii
  182.             x = rv1(i)
  183.             y = 0.0d0
  184.             if (i .eq. uk) go to 700
  185.             ip1 = i + 1
  186. c
  187.             do 680 j = ip1, uk
  188.                x = x - rm1(i,j) * rv1(j) + rm2(i,j) * rv2(j)
  189.                y = y - rm1(i,j) * rv2(j) - rm2(i,j) * rv1(j)
  190.   680       continue
  191. c
  192.   700       call cdiv(x,y,rm1(i,i),rm2(i,i),rv1(i),rv2(i))
  193.   720    continue
  194. c     .......... acceptance test for eigenvector
  195. c                and normalization ..........
  196.          its = its + 1
  197.          norm = 0.0d0
  198.          normv = 0.0d0
  199. c
  200.          do 780 i = 1, uk
  201.             x = pythag(rv1(i),rv2(i))
  202.             if (normv .ge. x) go to 760
  203.             normv = x
  204.             j = i
  205.   760       norm = norm + x
  206.   780    continue
  207. c
  208.          if (norm .lt. growto) go to 840
  209. c     .......... accept vector ..........
  210.          x = rv1(j)
  211.          y = rv2(j)
  212. c
  213.          do 820 i = 1, uk
  214.             call cdiv(rv1(i),rv2(i),x,y,zr(i,s),zi(i,s))
  215.   820    continue
  216. c
  217.          if (uk .eq. n) go to 940
  218.          j = uk + 1
  219.          go to 900
  220. c     .......... in-line procedure for choosing
  221. c                a new starting vector ..........
  222.   840    if (its .ge. uk) go to 880
  223.          x = ukroot
  224.          y = eps3 / (x + 1.0d0)
  225.          rv1(1) = eps3
  226. c
  227.          do 860 i = 2, uk
  228.   860    rv1(i) = y
  229. c
  230.          j = uk - its + 1
  231.          rv1(j) = rv1(j) - eps3 * x
  232.          go to 660
  233. c     .......... set error -- unaccepted eigenvector ..........
  234.   880    j = 1
  235.          ierr = -k
  236. c     .......... set remaining vector components to zero ..........
  237.   900    do 920 i = j, n
  238.             zr(i,s) = 0.0d0
  239.             zi(i,s) = 0.0d0
  240.   920    continue
  241. c
  242.   940    s = s + 1
  243.   980 continue
  244. c
  245.       go to 1001
  246. c     .......... set error -- underestimate of eigenvector
  247. c                space required ..........
  248.  1000 if (ierr .ne. 0) ierr = ierr - n
  249.       if (ierr .eq. 0) ierr = -(2 * n + 1)
  250.  1001 m = s - 1
  251.       return
  252.       end
  253.