home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / bandv.f < prev    next >
Text File  |  1996-09-28  |  11KB  |  322 lines
  1.       subroutine bandv(nm,n,mbw,a,e21,m,w,z,ierr,nv,rv,rv6)
  2. c
  3.       integer i,j,k,m,n,r,ii,ij,jj,kj,mb,m1,nm,nv,ij1,its,kj1,mbw,m21,
  4.      x        ierr,maxj,maxk,group
  5.       double precision a(nm,mbw),w(m),z(nm,m),rv(nv),rv6(n)
  6.       double precision u,v,uk,xu,x0,x1,e21,eps2,eps3,eps4,norm,order,
  7.      x       epslon,pythag
  8. c
  9. c     this subroutine finds those eigenvectors of a real symmetric
  10. c     band matrix corresponding to specified eigenvalues, using inverse
  11. c     iteration.  the subroutine may also be used to solve systems
  12. c     of linear equations with a symmetric or non-symmetric band
  13. c     coefficient matrix.
  14. c
  15. c     on input
  16. c
  17. c        nm must be set to the row dimension of two-dimensional
  18. c          array parameters as declared in the calling program
  19. c          dimension statement.
  20. c
  21. c        n is the order of the matrix.
  22. c
  23. c        mbw is the number of columns of the array a used to store the
  24. c          band matrix.  if the matrix is symmetric, mbw is its (half)
  25. c          band width, denoted mb and defined as the number of adjacent
  26. c          diagonals, including the principal diagonal, required to
  27. c          specify the non-zero portion of the lower triangle of the
  28. c          matrix.  if the subroutine is being used to solve systems
  29. c          of linear equations and the coefficient matrix is not
  30. c          symmetric, it must however have the same number of adjacent
  31. c          diagonals above the main diagonal as below, and in this
  32. c          case, mbw=2*mb-1.
  33. c
  34. c        a contains the lower triangle of the symmetric band input
  35. c          matrix stored as an n by mb array.  its lowest subdiagonal
  36. c          is stored in the last n+1-mb positions of the first column,
  37. c          its next subdiagonal in the last n+2-mb positions of the
  38. c          second column, further subdiagonals similarly, and finally
  39. c          its principal diagonal in the n positions of column mb.
  40. c          if the subroutine is being used to solve systems of linear
  41. c          equations and the coefficient matrix is not symmetric, a is
  42. c          n by 2*mb-1 instead with lower triangle as above and with
  43. c          its first superdiagonal stored in the first n-1 positions of
  44. c          column mb+1, its second superdiagonal in the first n-2
  45. c          positions of column mb+2, further superdiagonals similarly,
  46. c          and finally its highest superdiagonal in the first n+1-mb
  47. c          positions of the last column.
  48. c          contents of storages not part of the matrix are arbitrary.
  49. c
  50. c        e21 specifies the ordering of the eigenvalues and contains
  51. c            0.0d0 if the eigenvalues are in ascending order, or
  52. c            2.0d0 if the eigenvalues are in descending order.
  53. c          if the subroutine is being used to solve systems of linear
  54. c          equations, e21 should be set to 1.0d0 if the coefficient
  55. c          matrix is symmetric and to -1.0d0 if not.
  56. c
  57. c        m is the number of specified eigenvalues or the number of
  58. c          systems of linear equations.
  59. c
  60. c        w contains the m eigenvalues in ascending or descending order.
  61. c          if the subroutine is being used to solve systems of linear
  62. c          equations (a-w(r)*i)*x(r)=b(r), where i is the identity
  63. c          matrix, w(r) should be set accordingly, for r=1,2,...,m.
  64. c
  65. c        z contains the constant matrix columns (b(r),r=1,2,...,m), if
  66. c          the subroutine is used to solve systems of linear equations.
  67. c
  68. c        nv must be set to the dimension of the array parameter rv
  69. c          as declared in the calling program dimension statement.
  70. c
  71. c     on output
  72. c
  73. c        a and w are unaltered.
  74. c
  75. c        z contains the associated set of orthogonal eigenvectors.
  76. c          any vector which fails to converge is set to zero.  if the
  77. c          subroutine is used to solve systems of linear equations,
  78. c          z contains the solution matrix columns (x(r),r=1,2,...,m).
  79. c
  80. c        ierr is set to
  81. c          zero       for normal return,
  82. c          -r         if the eigenvector corresponding to the r-th
  83. c                     eigenvalue fails to converge, or if the r-th
  84. c                     system of linear equations is nearly singular.
  85. c
  86. c        rv and rv6 are temporary storage arrays.  note that rv is
  87. c          of dimension at least n*(2*mb-1).  if the subroutine
  88. c          is being used to solve systems of linear equations, the
  89. c          determinant (up to sign) of a-w(m)*i is available, upon
  90. c          return, as the product of the first n elements of rv.
  91. c
  92. c     calls pythag for  dsqrt(a*a + b*b) .
  93. c
  94. c     questions and comments should be directed to burton s. garbow,
  95. c     mathematics and computer science div, argonne national laboratory
  96. c
  97. c     this version dated august 1983.
  98. c
  99. c     ------------------------------------------------------------------
  100. c
  101.       ierr = 0
  102.       if (m .eq. 0) go to 1001
  103.       mb = mbw
  104.       if (e21 .lt. 0.0d0) mb = (mbw + 1) / 2
  105.       m1 = mb - 1
  106.       m21 = m1 + mb
  107.       order = 1.0d0 - dabs(e21)
  108. c     .......... find vectors by inverse iteration ..........
  109.       do 920 r = 1, m
  110.          its = 1
  111.          x1 = w(r)
  112.          if (r .ne. 1) go to 100
  113. c     .......... compute norm of matrix ..........
  114.          norm = 0.0d0
  115. c
  116.          do 60 j = 1, mb
  117.             jj = mb + 1 - j
  118.             kj = jj + m1
  119.             ij = 1
  120.             v = 0.0d0
  121. c
  122.             do 40 i = jj, n
  123.                v = v + dabs(a(i,j))
  124.                if (e21 .ge. 0.0d0) go to 40
  125.                v = v + dabs(a(ij,kj))
  126.                ij = ij + 1
  127.    40       continue
  128. c
  129.             norm = dmax1(norm,v)
  130.    60    continue
  131. c
  132.          if (e21 .lt. 0.0d0) norm = 0.5d0 * norm
  133. c     .......... eps2 is the criterion for grouping,
  134. c                eps3 replaces zero pivots and equal
  135. c                roots are modified by eps3,
  136. c                eps4 is taken very small to avoid overflow ..........
  137.          if (norm .eq. 0.0d0) norm = 1.0d0
  138.          eps2 = 1.0d-3 * norm * dabs(order)
  139.          eps3 = epslon(norm)
  140.          uk = n
  141.          uk = dsqrt(uk)
  142.          eps4 = uk * eps3
  143.    80    group = 0
  144.          go to 120
  145. c     .......... look for close or coincident roots ..........
  146.   100    if (dabs(x1-x0) .ge. eps2) go to 80
  147.          group = group + 1
  148.          if (order * (x1 - x0) .le. 0.0d0) x1 = x0 + order * eps3
  149. c     .......... expand matrix, subtract eigenvalue,
  150. c                and initialize vector ..........
  151.   120    do 200 i = 1, n
  152.             ij = i + min0(0,i-m1) * n
  153.             kj = ij + mb * n
  154.             ij1 = kj + m1 * n
  155.             if (m1 .eq. 0) go to 180
  156. c
  157.             do 150 j = 1, m1
  158.                if (ij .gt. m1) go to 125
  159.                if (ij .gt. 0) go to 130
  160.                rv(ij1) = 0.0d0
  161.                ij1 = ij1 + n
  162.                go to 130
  163.   125          rv(ij) = a(i,j)
  164.   130          ij = ij + n
  165.                ii = i + j
  166.                if (ii .gt. n) go to 150
  167.                jj = mb - j
  168.                if (e21 .ge. 0.0d0) go to 140
  169.                ii = i
  170.                jj = mb + j
  171.   140          rv(kj) = a(ii,jj)
  172.                kj = kj + n
  173.   150       continue
  174. c
  175.   180       rv(ij) = a(i,mb) - x1
  176.             rv6(i) = eps4
  177.             if (order .eq. 0.0d0) rv6(i) = z(i,r)
  178.   200    continue
  179. c
  180.          if (m1 .eq. 0) go to 600
  181. c     .......... elimination with interchanges ..........
  182.          do 580 i = 1, n
  183.             ii = i + 1
  184.             maxk = min0(i+m1-1,n)
  185.             maxj = min0(n-i,m21-2) * n
  186. c
  187.             do 360 k = i, maxk
  188.                kj1 = k
  189.                j = kj1 + n
  190.                jj = j + maxj
  191. c
  192.                do 340 kj = j, jj, n
  193.                   rv(kj1) = rv(kj)
  194.                   kj1 = kj
  195.   340          continue
  196. c
  197.                rv(kj1) = 0.0d0
  198.   360       continue
  199. c
  200.             if (i .eq. n) go to 580
  201.             u = 0.0d0
  202.             maxk = min0(i+m1,n)
  203.             maxj = min0(n-ii,m21-2) * n
  204. c
  205.             do 450 j = i, maxk
  206.                if (dabs(rv(j)) .lt. dabs(u)) go to 450
  207.                u = rv(j)
  208.                k = j
  209.   450       continue
  210. c
  211.             j = i + n
  212.             jj = j + maxj
  213.             if (k .eq. i) go to 520
  214.             kj = k
  215. c
  216.             do 500 ij = i, jj, n
  217.                v = rv(ij)
  218.                rv(ij) = rv(kj)
  219.                rv(kj) = v
  220.                kj = kj + n
  221.   500       continue
  222. c
  223.             if (order .ne. 0.0d0) go to 520
  224.             v = rv6(i)
  225.             rv6(i) = rv6(k)
  226.             rv6(k) = v
  227.   520       if (u .eq. 0.0d0) go to 580
  228. c
  229.             do 560 k = ii, maxk
  230.                v = rv(k) / u
  231.                kj = k
  232. c
  233.                do 540 ij = j, jj, n
  234.                   kj = kj + n
  235.                   rv(kj) = rv(kj) - v * rv(ij)
  236.   540          continue
  237. c
  238.                if (order .eq. 0.0d0) rv6(k) = rv6(k) - v * rv6(i)
  239.   560       continue
  240. c
  241.   580    continue
  242. c     .......... back substitution
  243. c                for i=n step -1 until 1 do -- ..........
  244.   600    do 630 ii = 1, n
  245.             i = n + 1 - ii
  246.             maxj = min0(ii,m21)
  247.             if (maxj .eq. 1) go to 620
  248.             ij1 = i
  249.             j = ij1 + n
  250.             jj = j + (maxj - 2) * n
  251. c
  252.             do 610 ij = j, jj, n
  253.                ij1 = ij1 + 1
  254.                rv6(i) = rv6(i) - rv(ij) * rv6(ij1)
  255.   610       continue
  256. c
  257.   620       v = rv(i)
  258.             if (dabs(v) .ge. eps3) go to 625
  259. c     .......... set error -- nearly singular linear system ..........
  260.             if (order .eq. 0.0d0) ierr = -r
  261.             v = dsign(eps3,v)
  262.   625       rv6(i) = rv6(i) / v
  263.   630    continue
  264. c
  265.          xu = 1.0d0
  266.          if (order .eq. 0.0d0) go to 870
  267. c     .......... orthogonalize with respect to previous
  268. c                members of group ..........
  269.          if (group .eq. 0) go to 700
  270. c
  271.          do 680 jj = 1, group
  272.             j = r - group - 1 + jj
  273.             xu = 0.0d0
  274. c
  275.             do 640 i = 1, n
  276.   640       xu = xu + rv6(i) * z(i,j)
  277. c
  278.             do 660 i = 1, n
  279.   660       rv6(i) = rv6(i) - xu * z(i,j)
  280. c
  281.   680    continue
  282. c
  283.   700    norm = 0.0d0
  284. c
  285.          do 720 i = 1, n
  286.   720    norm = norm + dabs(rv6(i))
  287. c
  288.          if (norm .ge. 0.1d0) go to 840
  289. c     .......... in-line procedure for choosing
  290. c                a new starting vector ..........
  291.          if (its .ge. n) go to 830
  292.          its = its + 1
  293.          xu = eps4 / (uk + 1.0d0)
  294.          rv6(1) = eps4
  295. c
  296.          do 760 i = 2, n
  297.   760    rv6(i) = xu
  298. c
  299.          rv6(its) = rv6(its) - eps4 * uk
  300.          go to 600
  301. c     .......... set error -- non-converged eigenvector ..........
  302.   830    ierr = -r
  303.          xu = 0.0d0
  304.          go to 870
  305. c     .......... normalize so that sum of squares is
  306. c                1 and expand to full order ..........
  307.   840    u = 0.0d0
  308. c
  309.          do 860 i = 1, n
  310.   860    u = pythag(u,rv6(i))
  311. c
  312.          xu = 1.0d0 / u
  313. c
  314.   870    do 900 i = 1, n
  315.   900    z(i,r) = rv6(i) * xu
  316. c
  317.          x0 = x1
  318.   920 continue
  319. c
  320.  1001 return
  321.       end
  322.