Revert to previous commit

[mothur.git] / dpbsl.f

diff --git a/dpbsl.f b/dpbsl.f
new file mode 100644 (file)
index 0000000..d910dee
--- /dev/null
+++ b/dpbsl.f
@@ -0,0 +1,83 @@
+      subroutine dpbsl(abd,lda,n,m,b)
+
+      integer lda,n,m
+      double precision abd(lda,n),b(n)
+c
+c     dpbsl solves the double precision symmetric positive definite
+c     band system  a*x = b
+c     using the factors computed by dpbco or dpbfa.
+c
+c     on entry
+c
+c        abd     double precision(lda, n)
+c                the output from dpbco or dpbfa.
+c
+c        lda     integer
+c                the leading dimension of the array  abd .
+c
+c        n       integer
+c                the order of the matrix  a .
+c
+c        m       integer
+c                the number of diagonals above the main diagonal.
+c
+c        b       double precision(n)
+c                the right hand side vector.
+c
+c     on return
+c
+c        b       the solution vector  x .
+c
+c     error condition
+c
+c        a division by zero will occur if the input factor contains
+c        a zero on the diagonal.  technically this indicates
+c        singularity but it is usually caused by improper subroutine
+c        arguments.  it will not occur if the subroutines are called
+c        correctly and  info .eq. 0 .
+c
+c     to compute  inverse(a) * c  where  c  is a matrix
+c     with  p  columns
+c           call dpbco(abd,lda,n,rcond,z,info)
+c           if (rcond is too small .or. info .ne. 0) go to ...
+c           do 10 j = 1, p
+c              call dpbsl(abd,lda,n,c(1,j))
+c        10 continue
+c
+c     linpack.  this version dated 08/14/78 .
+c     cleve moler, university of new mexico, argonne national lab.
+c
+c     subroutines and functions
+c
+c     blas daxpy,ddot
+c     fortran min0
+c
+c     internal variables
+c
+      double precision ddot,t
+      integer k,kb,la,lb,lm
+c
+c     solve trans(r)*y = b
+c
+      do 10 k = 1, n
+         lm = min0(k-1,m)
+         la = m + 1 - lm
+         lb = k - lm
+         t = ddot(lm,abd(la,k),1,b(lb),1)
+         b(k) = (b(k) - t)/abd(m+1,k)
+   10 continue
+c
+c     solve r*x = y
+c
+      do 20 kb = 1, n
+         k = n + 1 - kb
+         lm = min0(k-1,m)
+         la = m + 1 - lm
+         lb = k - lm
+         b(k) = b(k)/abd(m+1,k)
+         t = -b(k)
+         call daxpy(lm,t,abd(la,k),1,b(lb),1)
+   20 continue
+      return
+      end
+