dpbsl.f

   1       subroutine dpbsl(abd,lda,n,m,b)
   2
   3       integer lda,n,m
   4       double precision abd(lda,n),b(n)
   5 c
   6 c     dpbsl solves the double precision symmetric positive definite
   7 c     band system  a*x = b
   8 c     using the factors computed by dpbco or dpbfa.
   9 c
  10 c     on entry
  11 c
  12 c        abd     double precision(lda, n)
  13 c                the output from dpbco or dpbfa.
  14 c
  15 c        lda     integer
  16 c                the leading dimension of the array  abd .
  17 c
  18 c        n       integer
  19 c                the order of the matrix  a .
  20 c
  21 c        m       integer
  22 c                the number of diagonals above the main diagonal.
  23 c
  24 c        b       double precision(n)
  25 c                the right hand side vector.
  26 c
  27 c     on return
  28 c
  29 c        b       the solution vector  x .
  30 c
  31 c     error condition
  32 c
  33 c        a division by zero will occur if the input factor contains
  34 c        a zero on the diagonal.  technically this indicates
  35 c        singularity but it is usually caused by improper subroutine
  36 c        arguments.  it will not occur if the subroutines are called
  37 c        correctly and  info .eq. 0 .
  38 c
  39 c     to compute  inverse(a) * c  where  c  is a matrix
  40 c     with  p  columns
  41 c           call dpbco(abd,lda,n,rcond,z,info)
  42 c           if (rcond is too small .or. info .ne. 0) go to ...
  43 c           do 10 j = 1, p
  44 c              call dpbsl(abd,lda,n,c(1,j))
  45 c        10 continue
  46 c
  47 c     linpack.  this version dated 08/14/78 .
  48 c     cleve moler, university of new mexico, argonne national lab.
  49 c
  50 c     subroutines and functions
  51 c
  52 c     blas daxpy,ddot
  53 c     fortran min0
  54 c
  55 c     internal variables
  56 c
  57       double precision ddot,t
  58       integer k,kb,la,lb,lm
  59 c
  60 c     solve trans(r)*y = b
  61 c
  62       do 10 k = 1, n
  63          lm = min0(k-1,m)
  64          la = m + 1 - lm
  65          lb = k - lm
  66          t = ddot(lm,abd(la,k),1,b(lb),1)
  67          b(k) = (b(k) - t)/abd(m+1,k)
  68    10 continue
  69 c
  70 c     solve r*x = y
  71 c
  72       do 20 kb = 1, n
  73          k = n + 1 - kb
  74          lm = min0(k-1,m)
  75          la = m + 1 - lm
  76          lb = k - lm
  77          b(k) = b(k)/abd(m+1,k)
  78          t = -b(k)
  79          call daxpy(lm,t,abd(la,k),1,b(lb),1)
  80    20 continue
  81       return
  82       end
  83