home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / reduc2.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  117 lines
  1.       subroutine reduc2(nm,n,a,b,dl,ierr)
  2. c
  3.       integer i,j,k,n,i1,j1,nm,nn,ierr
  4.       double precision a(nm,n),b(nm,n),dl(n)
  5.       double precision x,y
  6. c
  7. c     this subroutine is a translation of the algol procedure reduc2,
  8. c     num. math. 11, 99-110(1968) by martin and wilkinson.
  9. c     handbook for auto. comp., vol.ii-linear algebra, 303-314(1971).
  10. c
  11. c     this subroutine reduces the generalized symmetric eigenproblems
  12. c     abx=(lambda)x or bay=(lambda)y, where b is positive definite,
  13. c     to the standard symmetric eigenproblem using the cholesky
  14. c     factorization of b.
  15. c
  16. c     on input
  17. c
  18. c        nm must be set to the row dimension of two-dimensional
  19. c          array parameters as declared in the calling program
  20. c          dimension statement.
  21. c
  22. c        n is the order of the matrices a and b.  if the cholesky
  23. c          factor l of b is already available, n should be prefixed
  24. c          with a minus sign.
  25. c
  26. c        a and b contain the real symmetric input matrices.  only the
  27. c          full upper triangles of the matrices need be supplied.  if
  28. c          n is negative, the strict lower triangle of b contains,
  29. c          instead, the strict lower triangle of its cholesky factor l.
  30. c
  31. c        dl contains, if n is negative, the diagonal elements of l.
  32. c
  33. c     on output
  34. c
  35. c        a contains in its full lower triangle the full lower triangle
  36. c          of the symmetric matrix derived from the reduction to the
  37. c          standard form.  the strict upper triangle of a is unaltered.
  38. c
  39. c        b contains in its strict lower triangle the strict lower
  40. c          triangle of its cholesky factor l.  the full upper
  41. c          triangle of b is unaltered.
  42. c
  43. c        dl contains the diagonal elements of l.
  44. c
  45. c        ierr is set to
  46. c          zero       for normal return,
  47. c          7*n+1      if b is not positive definite.
  48. c
  49. c     questions and comments should be directed to burton s. garbow,
  50. c     mathematics and computer science div, argonne national laboratory
  51. c
  52. c     this version dated august 1983.
  53. c
  54. c     ------------------------------------------------------------------
  55. c
  56.       ierr = 0
  57.       nn = iabs(n)
  58.       if (n .lt. 0) go to 100
  59. c     .......... form l in the arrays b and dl ..........
  60.       do 80 i = 1, n
  61.          i1 = i - 1
  62. c
  63.          do 80 j = i, n
  64.             x = b(i,j)
  65.             if (i .eq. 1) go to 40
  66. c
  67.             do 20 k = 1, i1
  68.    20       x = x - b(i,k) * b(j,k)
  69. c
  70.    40       if (j .ne. i) go to 60
  71.             if (x .le. 0.0d0) go to 1000
  72.             y = dsqrt(x)
  73.             dl(i) = y
  74.             go to 80
  75.    60       b(j,i) = x / y
  76.    80 continue
  77. c     .......... form the lower triangle of a*l
  78. c                in the lower triangle of the array a ..........
  79.   100 do 200 i = 1, nn
  80.          i1 = i + 1
  81. c
  82.          do 200 j = 1, i
  83.             x = a(j,i) * dl(j)
  84.             if (j .eq. i) go to 140
  85.             j1 = j + 1
  86. c
  87.             do 120 k = j1, i
  88.   120       x = x + a(k,i) * b(k,j)
  89. c
  90.   140       if (i .eq. nn) go to 180
  91. c
  92.             do 160 k = i1, nn
  93.   160       x = x + a(i,k) * b(k,j)
  94. c
  95.   180       a(i,j) = x
  96.   200 continue
  97. c     .......... pre-multiply by transpose(l) and overwrite ..........
  98.       do 300 i = 1, nn
  99.          i1 = i + 1
  100.          y = dl(i)
  101. c
  102.          do 300 j = 1, i
  103.             x = y * a(i,j)
  104.             if (i .eq. nn) go to 280
  105. c
  106.             do 260 k = i1, nn
  107.   260       x = x + a(k,j) * b(k,i)
  108. c
  109.   280       a(i,j) = x
  110.   300 continue
  111. c
  112.       go to 1001
  113. c     .......... set error -- b is not positive definite ..........
  114.  1000 ierr = 7 * n + 1
  115.  1001 return
  116.       end
  117.