dpbfa.f

   1       subroutine dpbfa(abd,lda,n,m,info)
   2
   3       integer lda,n,m,info
   4       double precision abd(lda,n)
   5 c
   6 c     dpbfa factors a double precision symmetric positive definite
   7 c     matrix stored in band form.
   8 c
   9 c     dpbfa is usually called by dpbco, but it can be called
  10 c     directly with a saving in time if  rcond  is not needed.
  11 c
  12 c     on entry
  13 c
  14 c        abd     double precision(lda, n)
  15 c                the matrix to be factored.  the columns of the upper
  16 c                triangle are stored in the columns of abd and the
  17 c                diagonals of the upper triangle are stored in the
  18 c                rows of abd .  see the comments below for details.
  19 c
  20 c        lda     integer
  21 c                the leading dimension of the array  abd .
  22 c                lda must be .ge. m + 1 .
  23 c
  24 c        n       integer
  25 c                the order of the matrix  a .
  26 c
  27 c        m       integer
  28 c                the number of diagonals above the main diagonal.
  29 c                0 .le. m .lt. n .
  30 c
  31 c     on return
  32 c
  33 c        abd     an upper triangular matrix  r , stored in band
  34 c                form, so that  a = trans(r)*r .
  35 c
  36 c        info    integer
  37 c                = 0  for normal return.
  38 c                = k  if the leading minor of order  k  is not
  39 c                     positive definite.
  40 c
  41 c     band storage
  42 c
  43 c           if  a  is a symmetric positive definite band matrix,
  44 c           the following program segment will set up the input.
  45 c
  46 c                   m = (band width above diagonal)
  47 c                   do 20 j = 1, n
  48 c                      i1 = max0(1, j-m)
  49 c                      do 10 i = i1, j
  50 c                         k = i-j+m+1
  51 c                         abd(k,j) = a(i,j)
  52 c                10    continue
  53 c                20 continue
  54 c
  55 c     linpack.  this version dated 08/14/78 .
  56 c     cleve moler, university of new mexico, argonne national lab.
  57 c
  58 c     subroutines and functions
  59 c
  60 c     blas ddot
  61 c     fortran max0,sqrt
  62 c
  63 c     internal variables
  64 c
  65       double precision ddot,t, a(10), b(10), temp
  66       double precision s
  67       integer ik,j,jk,k,mu
  68 c     begin block with ...exits to 40
  69 c
  70 c
  71            do i = 1, 10
  72          a(i) = i
  73          b(i) = 2*i
  74 c         PRINT *, 'a = ', a(i), 'i = ', i
  75 c         PRINT *, 'b = ', b(i)
  76       end do
  77
  78           temp = ddot(10,a,1,b,1)
  79
  80          do 30 j = 1, n
  81             info = j
  82             s = 0.0d0
  83             ik = m + 1
  84             jk = max0(j-m,1)
  85             mu = max0(m+2-j,1)
  86             if (m .lt. mu) go to 20
  87             do 10 k = mu, m
  88
  89                t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1)
  90                t = t/abd(m+1,jk)
  91                abd(k,j) = t
  92                s = s + t*t
  93                ik = ik - 1
  94                jk = jk + 1
  95    10       continue
  96    20       continue
  97
  98             s = abd(m+1,j) - s
  99 c     ......exit
 100             if (s .le. 0.0d0) go to 40
 101
 102             abd(m+1,j) = sqrt(s)
 103
 104    30    continue
 105          info = 0
 106    40 continue
 107       return
 108       end
 109