home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tests / rblres.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  156 lines
  1.       subroutine rblres(nm,n,mbw,ip,a,b,x,norm,e21)
  2. c
  3.       double precision a(nm,mbw),b(nm,ip),x(nm,ip),norm(n)
  4.       double precision norma,normx,normr,sum,suma,e21
  5. c
  6. c     this subroutine forms the 1-norm of the residual matrix
  7. c     a*x-b  where  a  is a real band matrix, which may be symmetric,
  8. c     x  is a real  n  by  ip  matrix, and   b is a real
  9. c     n  by  ip  matrix.
  10. c
  11. c     this subroutine is catalogued as eispdrv4(rblres).
  12. c
  13. c     input.
  14. c        nm is the row dimension of two-dimensional array parameters
  15. c           as declared in the calling program dimension statement;
  16. c
  17. c        n is the row dimension of the matrix  a;
  18. c
  19. c        mbw is the number of columns of the array a used to store the
  20. c           band matrix.  if the matrix is symmetric, mbw is its (half)
  21. c           band width, denoted mb and defined as the number of adjacent
  22. c           diagonals, including the principal diagonal, required to
  23. c           specify the non-zero portion of the lower triangle of the
  24. c           matrix.  if the matrix is non-symmetric it must have the sam
  25. c           number of adjacent diagonals above the main diagonal as
  26. c           below.  in this case  mbw = 2*mb-1;
  27. c
  28. c        ip  is the column dimension of  b  and  x;
  29. c
  30. c        a(n,mbw) is an array which contains in its columns the
  31. c           subdiagonals and diagonal of the symmetric band
  32. c           matrix. if  a  is non-symmetric it also contains
  33. c           the super-diagonals;
  34. c
  35. c        b  is a real full  n  by  ip  matrix;
  36. c
  37. c        x  is a real full  n  by  ip  matrix;
  38. c
  39. c        e21 tells if the matrix is symmetric.  if e21 equals one the
  40. c           matrix is symmetric and if e21 equals minus one the
  41. c           matrix is non-symmetric.
  42. c
  43. c     output.
  44. c
  45. c        norm(n)  is an array such that for each k,
  46. c               norm(k) = !!a*x(k)-b(k)!!/(!!a!!*!!x(k)!!) .
  47. c
  48. c     ------------------------------------------------------------------
  49. c
  50.       if( e21 .eq. 1.0d0 ) go to 75
  51.       mbwt = mbw
  52.       mb = mbw/2 + 1
  53.       norma = 0.0d0
  54. c
  55.       do 35 i = 1,n
  56.          l = max0(mb - i + 1,1)
  57.          m = min0(mb - i + n,mbwt)
  58.          suma = 0.0d0
  59. c
  60.          do 34 j = l,m
  61.             suma = suma + dabs(a(i,j))
  62.    34    continue
  63. c
  64.          norma = dmax1(norma,suma)
  65.    35 continue
  66. c
  67.       do 70 ipp = 1,ip
  68.          normr = 0.0d0
  69.          normx = 0.0d0
  70.          do 50 i = 1,n
  71.             l = max0(mb - i + 1,1)
  72.             m = min0(mb - i + n,mbwt)
  73.             k = max0(-mb + 1 + i,1)
  74.             sum = -b(i,ipp)
  75.             suma = 0.0d0
  76. c
  77.             do 40 j = l,m
  78.                sum = sum + a(i,j)*x(k,ipp)
  79.                k = k + 1
  80.                suma = suma + dabs(a(i,j))
  81.    40       continue
  82. c
  83.             normr = normr + dabs(sum)
  84.             normx = normx + dabs(x(i,ipp))
  85.    50    continue
  86. c
  87.          if( normx .eq. 0.0d0 ) normx = 1.0d0
  88.          if( norma .eq. 0.0d0 ) norma = 1.0d0
  89.          norm(ipp) = normr/(norma*normx)
  90.    70 continue
  91. c
  92.       return
  93. c
  94.    75 mb = mbw
  95.       norma = 0.0d0
  96.       mb1 = mb - 1
  97. c
  98.       do 110 i = 1,n
  99.          suma = 0.0d0
  100.          if( i .eq. 1 ) go to 90
  101.          lstart = max0(1,mb + 1 - i)
  102. c
  103.          do 80 l = lstart,mb1
  104.             suma = suma + dabs(a(i,l))
  105.    80    continue
  106. c
  107.    90    lstop = min0(mb,n + 1 - i)
  108.          j = i
  109. c
  110.          do 100 l = 1,lstop
  111.             l1 = mb + 1 - l
  112.             suma = suma + dabs(a(j,l1))
  113.             j = j + 1
  114.   100    continue
  115. c
  116.          norma = dmax1(norma,suma)
  117.   110 continue
  118. c
  119.       do 170 i = 1,ip
  120.          normr = 0.0d0
  121.          normx = 0.0d0
  122. c
  123.          do 150 l = 1,n
  124.             sum = -b(l,i)
  125.             j = max0(0,l - mb)
  126.             if( l .eq. 1 ) go to 130
  127.             kstart = max0(1,mb + 1 - l)
  128. c
  129.             do 120 k = kstart,mb1
  130.                j = j + 1
  131.                sum = sum + a(l,k)*x(j,i)
  132.   120       continue
  133. c
  134.   130       kstop = min0(mb,n + 1 - l)
  135. c
  136.             do 140 k = 1,kstop
  137.                j = j + 1
  138.                k1 = mb + 1 - k
  139.                sum = sum + a(j,k1)*x(j,i)
  140.   140       continue
  141. c
  142.             normr = normr + dabs(sum)
  143.   150    continue
  144. c
  145.          do 160 k = 1,n
  146.             normx = normx + dabs(x(k,i))
  147.   160    continue
  148. c
  149.          if( normx .eq. 0.0d0 ) normx = 1.0d0
  150.          if( norma .eq. 0.0d0 ) norma = 1.0d0
  151.          norm(i) = normr/(norma*normx)
  152.   170 continue
  153.       return
  154. c
  155.       end
  156.