home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / ex / rsgabr.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  126 lines
  1.       subroutine rsgabr(nm,n,m,a,b,w,z,norm,resdul,r)
  2. c
  3.       double precision norm(m),w(m),a(nm,n),z(nm,m),norma,tnorm,s,sum, 
  4.      x       sumz, suma, resdul, b(nm,n), normb, sumb, r(n)
  5. c
  6. c     this subroutine forms the 1-norm of the residual matrix
  7. c     a*b*z-z*diag(w)  where  a  is a real symmetric matrix,
  8. c     b  is a real symmetric matrix ,  w  is a vector which
  9. c     contains  m  eigenvalues of the eigenproblem, a*b*z-z*diag(w),
  10. c     and  z  is an array which contains the  m  corresponding
  11. c     eigenvectors of the eigenproblem. all norms appearing in
  12. c     the comments below are 1-norms.
  13. c
  14. c     this subroutine is catalogued as eispdrv4(rsgabr).
  15. c
  16. c     input.
  17. c
  18. c        nm is the row dimension of two-dimensional array parameters
  19. c           as declared in the calling program dimension statement;
  20. c
  21. c        m is the number of eigenvectors for which residuals are
  22. c           desired;
  23. c
  24. c        n is the order of the matrix  a;
  25. c
  26. c        a(nm,n) is a real symmetric matrix.  only the full upper
  27. c           triangle need be supplied;
  28. c
  29. c        b(nm,n) is a real symmetric matrix.  only the full upper
  30. c           triangle need be supplied;
  31. c
  32. c        z(nm,m) is a real matrix whose first  m  columns contain the
  33. c          approximate eigenvectors of the eigenproblem;
  34. c
  35. c        w(m) is a vector whose first  m   components contain the
  36. c           approximate eigenvalues of the eigenproblem.  w(i) is
  37. c           associated with the i-th  column of  z.
  38. c
  39. c
  40. c     output.
  41. c
  42. c        z(nm,m) is an array which contains  m  normalized
  43. c           approximate eigenvectors of the eigenproblem. the
  44. c           eigenvectors are normalized using the 1-norm in such a way
  45. c           that the first element whose magnitude is larger than the
  46. c           norm of the eigenvector divided by  n  is positive;
  47. c
  48. c        a(nm,n) is altered by making the lower triangle of  a
  49. c           symmetric with the upper triangle of  a;
  50. c
  51. c        b(nm,n) is altered by making the lower triangle of  b
  52. c           symmetric with the upper triangle of  b;
  53. c
  54. c        norm(m) is an array such that for each  k,
  55. c           norm(k) = !!a*b*z(k)-z(k)*w(k)!!/(!!a!!*!!b!!*!!z(k)!!)
  56. c           where  z(k)  is the k-th eigenvector;
  57. c
  58. c        resdul is the real number
  59. c           !!a*b*z-z*diag(w)!!/(!!a!!*!!b!!*!!z!!);
  60. c
  61. c        r(n) is a temporary storage array used to store the product
  62. c           b*z.
  63. c
  64. c     ----------------------------------------------------------------
  65. c
  66.       resdul = 0.0d0
  67.       if( m .eq. 0 ) return
  68.       norma = 0.0d0
  69.       normb = 0.0d0
  70. c
  71.       do 40 i=1,n
  72.          sumb = 0.0d0
  73.          suma = 0.0d0
  74.          if(i .eq. 1) go to 20
  75. c
  76.          do 10 l=2,i
  77.            a(i,l-1) =a(l-1,i)
  78.            b(i,l-1) =b(l-1,i)
  79.            sumb =sumb + dabs(b(l-1,i))
  80.    10      suma =suma + dabs(a(l-1,i))
  81. c
  82.    20    do 30 l=i,n
  83.            suma =suma + dabs(a(i,l))
  84.    30      sumb =sumb + dabs(b(i,l))
  85. c
  86.          norma = dmax1(norma,suma)
  87.    40    normb = dmax1(normb,sumb)
  88. c
  89.       norma = norma*normb
  90.       if(norma .eq. 0.0d0) norma = 1.0d0
  91. c
  92.       do 100 i=1,m
  93.          s = 0.0d0
  94.          sumz = 0.0d0
  95.          do 55 l = 1,n
  96.            sum = 0.0d0
  97.            do 50 k = 1,n
  98.    50        sum = sum + b(l,k)*z(k,i)
  99.            sumz = sumz + dabs(z(l,i))
  100.    55      r(l) = sum
  101. c
  102.          do 65 l = 1,n
  103.            sum = - w(i)*z(l,i)
  104.            do 60 k = 1,n
  105.    60        sum = sum + r(k)*a(l,k)
  106.    65      s = s + dabs(sum)
  107.          norm(i) = sumz
  108.          if(sumz .eq. 0.0d0) go to 100
  109. c        ..........this loop will never be completed since there
  110. c                  will always exist an element in the vector z(i)
  111. c                  larger than !!z(i)!!/n..........
  112.          do 70 l=1,n
  113.             if(dabs(z(l,i)) .ge. norm(i)/dfloat(n)) go to 80
  114.    70       continue
  115. c
  116.    80    tnorm = dsign(norm(i),z(l,i))
  117. c
  118.          do 90 l=1,n
  119.    90       z(l,i) = z(l,i)/tnorm
  120. c
  121.          norm(i) = s/(norm(i)*norma)
  122.   100    resdul = dmax1(norm(i), resdul)
  123. c
  124.       return
  125.       end
  126.