home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / reduc.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  116 lines
  1.       subroutine reduc(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 reduc1,
  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 eigenproblem
  12. c     ax=(lambda)bx, where b is positive definite, to the standard
  13. c     symmetric eigenproblem using the cholesky factorization of b.
  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 matrices a and b.  if the cholesky
  22. c          factor l of b is already available, n should be prefixed
  23. c          with a minus sign.
  24. c
  25. c        a and b contain the real symmetric input matrices.  only the
  26. c          full upper triangles of the matrices need be supplied.  if
  27. c          n is negative, the strict lower triangle of b contains,
  28. c          instead, the strict lower triangle of its cholesky factor l.
  29. c
  30. c        dl contains, if n is negative, the diagonal elements of l.
  31. c
  32. c     on output
  33. c
  34. c        a contains in its full lower triangle the full lower triangle
  35. c          of the symmetric matrix derived from the reduction to the
  36. c          standard form.  the strict upper triangle of a is unaltered.
  37. c
  38. c        b contains in its strict lower triangle the strict lower
  39. c          triangle of its cholesky factor l.  the full upper
  40. c          triangle of b is unaltered.
  41. c
  42. c        dl contains the diagonal elements of l.
  43. c
  44. c        ierr is set to
  45. c          zero       for normal return,
  46. c          7*n+1      if b is not positive definite.
  47. c
  48. c     questions and comments should be directed to burton s. garbow,
  49. c     mathematics and computer science div, argonne national laboratory
  50. c
  51. c     this version dated august 1983.
  52. c
  53. c     ------------------------------------------------------------------
  54. c
  55.       ierr = 0
  56.       nn = iabs(n)
  57.       if (n .lt. 0) go to 100
  58. c     .......... form l in the arrays b and dl ..........
  59.       do 80 i = 1, n
  60.          i1 = i - 1
  61. c
  62.          do 80 j = i, n
  63.             x = b(i,j)
  64.             if (i .eq. 1) go to 40
  65. c
  66.             do 20 k = 1, i1
  67.    20       x = x - b(i,k) * b(j,k)
  68. c
  69.    40       if (j .ne. i) go to 60
  70.             if (x .le. 0.0d0) go to 1000
  71.             y = dsqrt(x)
  72.             dl(i) = y
  73.             go to 80
  74.    60       b(j,i) = x / y
  75.    80 continue
  76. c     .......... form the transpose of the upper triangle of inv(l)*a
  77. c                in the lower triangle of the array a ..........
  78.   100 do 200 i = 1, nn
  79.          i1 = i - 1
  80.          y = dl(i)
  81. c
  82.          do 200 j = i, nn
  83.             x = a(i,j)
  84.             if (i .eq. 1) go to 180
  85. c
  86.             do 160 k = 1, i1
  87.   160       x = x - b(i,k) * a(j,k)
  88. c
  89.   180       a(j,i) = x / y
  90.   200 continue
  91. c     .......... pre-multiply by inv(l) and overwrite ..........
  92.       do 300 j = 1, nn
  93.          j1 = j - 1
  94. c
  95.          do 300 i = j, nn
  96.             x = a(i,j)
  97.             if (i .eq. j) go to 240
  98.             i1 = i - 1
  99. c
  100.             do 220 k = j, i1
  101.   220       x = x - a(k,j) * b(i,k)
  102. c
  103.   240       if (j .eq. 1) go to 280
  104. c
  105.             do 260 k = 1, j1
  106.   260       x = x - a(j,k) * b(i,k)
  107. c
  108.   280       a(i,j) = x / dl(i)
  109.   300 continue
  110. c
  111.       go to 1001
  112. c     .......... set error -- b is not positive definite ..........
  113.  1000 ierr = 7 * n + 1
  114.  1001 return
  115.       end
  116.