home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / invit.f < prev    next >
Text File  |  1996-09-28  |  13KB  |  386 lines
  1.       subroutine invit(nm,n,a,wr,wi,select,mm,m,z,ierr,rm1,rv1,rv2)
  2. c
  3.       integer i,j,k,l,m,n,s,ii,ip,mm,mp,nm,ns,n1,uk,ip1,its,km1,ierr
  4.       double precision a(nm,n),wr(n),wi(n),z(nm,mm),rm1(n,n),
  5.      x       rv1(n),rv2(n)
  6.       double precision t,w,x,y,eps3,norm,normv,epslon,growto,ilambd,
  7.      x       pythag,rlambd,ukroot
  8.       logical select(n)
  9. c
  10. c     this subroutine is a translation of the algol procedure invit
  11. c     by peters and wilkinson.
  12. c     handbook for auto. comp., vol.ii-linear algebra, 418-439(1971).
  13. c
  14. c     this subroutine finds those eigenvectors of a real upper
  15. c     hessenberg matrix corresponding to specified eigenvalues,
  16. c     using inverse iteration.
  17. c
  18. c     on input
  19. c
  20. c        nm must be set to the row dimension of two-dimensional
  21. c          array parameters as declared in the calling program
  22. c          dimension statement.
  23. c
  24. c        n is the order of the matrix.
  25. c
  26. c        a contains the hessenberg matrix.
  27. c
  28. c        wr and wi contain the real and imaginary parts, respectively,
  29. c          of the eigenvalues of the matrix.  the eigenvalues must be
  30. c          stored in a manner identical to that of subroutine  hqr,
  31. c          which recognizes possible splitting of the matrix.
  32. c
  33. c        select specifies the eigenvectors to be found. the
  34. c          eigenvector corresponding to the j-th eigenvalue is
  35. c          specified by setting select(j) to .true..
  36. c
  37. c        mm should be set to an upper bound for the number of
  38. c          columns required to store the eigenvectors to be found.
  39. c          note that two columns are required to store the
  40. c          eigenvector corresponding to a complex eigenvalue.
  41. c
  42. c     on output
  43. c
  44. c        a and wi 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        select may have been altered.  if the elements corresponding
  50. c          to a pair of conjugate complex eigenvalues were each
  51. c          initially set to .true., the program resets the second of
  52. c          the two elements to .false..
  53. c
  54. c        m is the number of columns actually used to store
  55. c          the eigenvectors.
  56. c
  57. c        z contains the real and imaginary parts of the eigenvectors.
  58. c          if the next selected eigenvalue is real, the next column
  59. c          of z contains its eigenvector.  if the eigenvalue is
  60. c          complex, the next two columns of z contain the real and
  61. c          imaginary parts of its eigenvector.  the eigenvectors are
  62. c          normalized so that the component of largest magnitude is 1.
  63. c          any vector which fails the acceptance test is set to zero.
  64. c
  65. c        ierr is set to
  66. c          zero       for normal return,
  67. c          -(2*n+1)   if more than mm columns of z are necessary
  68. c                     to store the eigenvectors corresponding to
  69. c                     the specified eigenvalues.
  70. c          -k         if the iteration corresponding to the k-th
  71. c                     value fails,
  72. c          -(n+k)     if both error situations occur.
  73. c
  74. c        rm1, rv1, and rv2 are temporary storage arrays.  note that rm1
  75. c          is square of dimension n by n and, augmented by two columns
  76. c          of z, is the transpose of the corresponding algol b array.
  77. c
  78. c     the algol procedure guessvec appears in invit in line.
  79. c
  80. c     calls cdiv for complex division.
  81. c     calls pythag for  dsqrt(a*a + b*b) .
  82. c
  83. c     questions and comments should be directed to burton s. garbow,
  84. c     mathematics and computer science div, argonne national laboratory
  85. c
  86. c     this version dated august 1983.
  87. c
  88. c     ------------------------------------------------------------------
  89. c
  90.       ierr = 0
  91.       uk = 0
  92.       s = 1
  93. c     .......... ip = 0, real eigenvalue
  94. c                     1, first of conjugate complex pair
  95. c                    -1, second of conjugate complex pair ..........
  96.       ip = 0
  97.       n1 = n - 1
  98. c
  99.       do 980 k = 1, n
  100.          if (wi(k) .eq. 0.0d0 .or. ip .lt. 0) go to 100
  101.          ip = 1
  102.          if (select(k) .and. select(k+1)) select(k+1) = .false.
  103.   100    if (.not. select(k)) go to 960
  104.          if (wi(k) .ne. 0.0d0) s = s + 1
  105.          if (s .gt. mm) go to 1000
  106.          if (uk .ge. k) go to 200
  107. c     .......... check for possible splitting ..........
  108.          do 120 uk = k, n
  109.             if (uk .eq. n) go to 140
  110.             if (a(uk+1,uk) .eq. 0.0d0) go to 140
  111.   120    continue
  112. c     .......... compute infinity norm of leading uk by uk
  113. c                (hessenberg) matrix ..........
  114.   140    norm = 0.0d0
  115.          mp = 1
  116. c
  117.          do 180 i = 1, uk
  118.             x = 0.0d0
  119. c
  120.             do 160 j = mp, uk
  121.   160       x = x + dabs(a(i,j))
  122. c
  123.             if (x .gt. norm) norm = x
  124.             mp = i
  125.   180    continue
  126. c     .......... eps3 replaces zero pivot in decomposition
  127. c                and close roots are modified by eps3 ..........
  128.          if (norm .eq. 0.0d0) norm = 1.0d0
  129.          eps3 = epslon(norm)
  130. c     .......... growto is the criterion for the growth ..........
  131.          ukroot = uk
  132.          ukroot = dsqrt(ukroot)
  133.          growto = 0.1d0 / ukroot
  134.   200    rlambd = wr(k)
  135.          ilambd = wi(k)
  136.          if (k .eq. 1) go to 280
  137.          km1 = k - 1
  138.          go to 240
  139. c     .......... perturb eigenvalue if it is close
  140. c                to any previous eigenvalue ..........
  141.   220    rlambd = rlambd + eps3
  142. c     .......... for i=k-1 step -1 until 1 do -- ..........
  143.   240    do 260 ii = 1, km1
  144.             i = k - ii
  145.             if (select(i) .and. dabs(wr(i)-rlambd) .lt. eps3 .and.
  146.      x         dabs(wi(i)-ilambd) .lt. eps3) go to 220
  147.   260    continue
  148. c
  149.          wr(k) = rlambd
  150. c     .......... perturb conjugate eigenvalue to match ..........
  151.          ip1 = k + ip
  152.          wr(ip1) = rlambd
  153. c     .......... form upper hessenberg a-rlambd*i (transposed)
  154. c                and initial real vector ..........
  155.   280    mp = 1
  156. c
  157.          do 320 i = 1, uk
  158. c
  159.             do 300 j = mp, uk
  160.   300       rm1(j,i) = a(i,j)
  161. c
  162.             rm1(i,i) = rm1(i,i) - rlambd
  163.             mp = i
  164.             rv1(i) = eps3
  165.   320    continue
  166. c
  167.          its = 0
  168.          if (ilambd .ne. 0.0d0) go to 520
  169. c     .......... real eigenvalue.
  170. c                triangular decomposition with interchanges,
  171. c                replacing zero pivots by eps3 ..........
  172.          if (uk .eq. 1) go to 420
  173. c
  174.          do 400 i = 2, uk
  175.             mp = i - 1
  176.             if (dabs(rm1(mp,i)) .le. dabs(rm1(mp,mp))) go to 360
  177. c
  178.             do 340 j = mp, uk
  179.                y = rm1(j,i)
  180.                rm1(j,i) = rm1(j,mp)
  181.                rm1(j,mp) = y
  182.   340       continue
  183. c
  184.   360       if (rm1(mp,mp) .eq. 0.0d0) rm1(mp,mp) = eps3
  185.             x = rm1(mp,i) / rm1(mp,mp)
  186.             if (x .eq. 0.0d0) go to 400
  187. c
  188.             do 380 j = i, uk
  189.   380       rm1(j,i) = rm1(j,i) - x * rm1(j,mp)
  190. c
  191.   400    continue
  192. c
  193.   420    if (rm1(uk,uk) .eq. 0.0d0) rm1(uk,uk) = eps3
  194. c     .......... back substitution for real vector
  195. c                for i=uk step -1 until 1 do -- ..........
  196.   440    do 500 ii = 1, uk
  197.             i = uk + 1 - ii
  198.             y = rv1(i)
  199.             if (i .eq. uk) go to 480
  200.             ip1 = i + 1
  201. c
  202.             do 460 j = ip1, uk
  203.   460       y = y - rm1(j,i) * rv1(j)
  204. c
  205.   480       rv1(i) = y / rm1(i,i)
  206.   500    continue
  207. c
  208.          go to 740
  209. c     .......... complex eigenvalue.
  210. c                triangular decomposition with interchanges,
  211. c                replacing zero pivots by eps3.  store imaginary
  212. c                parts in upper triangle starting at (1,3) ..........
  213.   520    ns = n - s
  214.          z(1,s-1) = -ilambd
  215.          z(1,s) = 0.0d0
  216.          if (n .eq. 2) go to 550
  217.          rm1(1,3) = -ilambd
  218.          z(1,s-1) = 0.0d0
  219.          if (n .eq. 3) go to 550
  220. c
  221.          do 540 i = 4, n
  222.   540    rm1(1,i) = 0.0d0
  223. c
  224.   550    do 640 i = 2, uk
  225.             mp = i - 1
  226.             w = rm1(mp,i)
  227.             if (i .lt. n) t = rm1(mp,i+1)
  228.             if (i .eq. n) t = z(mp,s-1)
  229.             x = rm1(mp,mp) * rm1(mp,mp) + t * t
  230.             if (w * w .le. x) go to 580
  231.             x = rm1(mp,mp) / w
  232.             y = t / w
  233.             rm1(mp,mp) = w
  234.             if (i .lt. n) rm1(mp,i+1) = 0.0d0
  235.             if (i .eq. n) z(mp,s-1) = 0.0d0
  236. c
  237.             do 560 j = i, uk
  238.                w = rm1(j,i)
  239.                rm1(j,i) = rm1(j,mp) - x * w
  240.                rm1(j,mp) = w
  241.                if (j .lt. n1) go to 555
  242.                l = j - ns
  243.                z(i,l) = z(mp,l) - y * w
  244.                z(mp,l) = 0.0d0
  245.                go to 560
  246.   555          rm1(i,j+2) = rm1(mp,j+2) - y * w
  247.                rm1(mp,j+2) = 0.0d0
  248.   560       continue
  249. c
  250.             rm1(i,i) = rm1(i,i) - y * ilambd
  251.             if (i .lt. n1) go to 570
  252.             l = i - ns
  253.             z(mp,l) = -ilambd
  254.             z(i,l) = z(i,l) + x * ilambd
  255.             go to 640
  256.   570       rm1(mp,i+2) = -ilambd
  257.             rm1(i,i+2) = rm1(i,i+2) + x * ilambd
  258.             go to 640
  259.   580       if (x .ne. 0.0d0) go to 600
  260.             rm1(mp,mp) = eps3
  261.             if (i .lt. n) rm1(mp,i+1) = 0.0d0
  262.             if (i .eq. n) z(mp,s-1) = 0.0d0
  263.             t = 0.0d0
  264.             x = eps3 * eps3
  265.   600       w = w / x
  266.             x = rm1(mp,mp) * w
  267.             y = -t * w
  268. c
  269.             do 620 j = i, uk
  270.                if (j .lt. n1) go to 610
  271.                l = j - ns
  272.                t = z(mp,l)
  273.                z(i,l) = -x * t - y * rm1(j,mp)
  274.                go to 615
  275.   610          t = rm1(mp,j+2)
  276.                rm1(i,j+2) = -x * t - y * rm1(j,mp)
  277.   615          rm1(j,i) = rm1(j,i) - x * rm1(j,mp) + y * t
  278.   620       continue
  279. c
  280.             if (i .lt. n1) go to 630
  281.             l = i - ns
  282.             z(i,l) = z(i,l) - ilambd
  283.             go to 640
  284.   630       rm1(i,i+2) = rm1(i,i+2) - ilambd
  285.   640    continue
  286. c
  287.          if (uk .lt. n1) go to 650
  288.          l = uk - ns
  289.          t = z(uk,l)
  290.          go to 655
  291.   650    t = rm1(uk,uk+2)
  292.   655    if (rm1(uk,uk) .eq. 0.0d0 .and. t .eq. 0.0d0) rm1(uk,uk) = eps3
  293. c     .......... back substitution for complex vector
  294. c                for i=uk step -1 until 1 do -- ..........
  295.   660    do 720 ii = 1, uk
  296.             i = uk + 1 - ii
  297.             x = rv1(i)
  298.             y = 0.0d0
  299.             if (i .eq. uk) go to 700
  300.             ip1 = i + 1
  301. c
  302.             do 680 j = ip1, uk
  303.                if (j .lt. n1) go to 670
  304.                l = j - ns
  305.                t = z(i,l)
  306.                go to 675
  307.   670          t = rm1(i,j+2)
  308.   675          x = x - rm1(j,i) * rv1(j) + t * rv2(j)
  309.                y = y - rm1(j,i) * rv2(j) - t * rv1(j)
  310.   680       continue
  311. c
  312.   700       if (i .lt. n1) go to 710
  313.             l = i - ns
  314.             t = z(i,l)
  315.             go to 715
  316.   710       t = rm1(i,i+2)
  317.   715       call cdiv(x,y,rm1(i,i),t,rv1(i),rv2(i))
  318.   720    continue
  319. c     .......... acceptance test for real or complex
  320. c                eigenvector and normalization ..........
  321.   740    its = its + 1
  322.          norm = 0.0d0
  323.          normv = 0.0d0
  324. c
  325.          do 780 i = 1, uk
  326.             if (ilambd .eq. 0.0d0) x = dabs(rv1(i))
  327.             if (ilambd .ne. 0.0d0) x = pythag(rv1(i),rv2(i))
  328.             if (normv .ge. x) go to 760
  329.             normv = x
  330.             j = i
  331.   760       norm = norm + x
  332.   780    continue
  333. c
  334.          if (norm .lt. growto) go to 840
  335. c     .......... accept vector ..........
  336.          x = rv1(j)
  337.          if (ilambd .eq. 0.0d0) x = 1.0d0 / x
  338.          if (ilambd .ne. 0.0d0) y = rv2(j)
  339. c
  340.          do 820 i = 1, uk
  341.             if (ilambd .ne. 0.0d0) go to 800
  342.             z(i,s) = rv1(i) * x
  343.             go to 820
  344.   800       call cdiv(rv1(i),rv2(i),x,y,z(i,s-1),z(i,s))
  345.   820    continue
  346. c
  347.          if (uk .eq. n) go to 940
  348.          j = uk + 1
  349.          go to 900
  350. c     .......... in-line procedure for choosing
  351. c                a new starting vector ..........
  352.   840    if (its .ge. uk) go to 880
  353.          x = ukroot
  354.          y = eps3 / (x + 1.0d0)
  355.          rv1(1) = eps3
  356. c
  357.          do 860 i = 2, uk
  358.   860    rv1(i) = y
  359. c
  360.          j = uk - its + 1
  361.          rv1(j) = rv1(j) - eps3 * x
  362.          if (ilambd .eq. 0.0d0) go to 440
  363.          go to 660
  364. c     .......... set error -- unaccepted eigenvector ..........
  365.   880    j = 1
  366.          ierr = -k
  367. c     .......... set remaining vector components to zero ..........
  368.   900    do 920 i = j, n
  369.             z(i,s) = 0.0d0
  370.             if (ilambd .ne. 0.0d0) z(i,s-1) = 0.0d0
  371.   920    continue
  372. c
  373.   940    s = s + 1
  374.   960    if (ip .eq. (-1)) ip = 0
  375.          if (ip .eq. 1) ip = -1
  376.   980 continue
  377. c
  378.       go to 1001
  379. c     .......... set error -- underestimate of eigenvector
  380. c                space required ..........
  381.  1000 if (ierr .ne. 0) ierr = ierr - n
  382.       if (ierr .eq. 0) ierr = -(2 * n + 1)
  383.  1001 m = s - 1 - iabs(ip)
  384.       return
  385.       end
  386.