subroutine bndchl( a,lda ,n,m ,info )
      real a(lda,n)
c driver for high-speed banded-cholesky decomposition.
c algorith is the same as dongarras matrix times vector method for dense
c matrices -only this is applied to banded systems.
c
      real rowi(1000)
c
c  the first m rows we use ordinary banded decomp
c  (we could decompose this into 2 matrix times vector operations
c   where one matrix was rectangular and one was triangular)
c  -but i didn't bother to.
c
         do 50 i=1,m
c
c        --- get rowi ---
c
         j1=max0(m+2-i,1)
            do 10 j=j1,m
10          rowi(j)=-a(j,i)
c
         nrows=min0(m+1,n+1-i)
            do 30 jj=j1,m
cdir$ ivdep
            do 30 ii=1,min0(jj,nrows)
30          a(m+2-ii,i+ii-1)=a(m+2-ii,i+ii-1)
     *      +rowi(jj)*a(jj+1-ii,i-1+ii)
c
      info=i
      if(a(m+1,i).le.0. ) go to 60
         a(m+1,i)=1./sqrt(a(m+1,i))
            do 40 ii=1,nrows-1
40          a(m+1-ii,i+ii)=a(m+1-ii,i+ii)*a(m+1,i)
c
50       continue
c
c     central rows of the matrix we have only triangular matrix times
c     times vector operations
c
      j1=1
      nrows=     m+1
         do 150 i=m+1,n-m
c
c        --- get rowi ---
c
            do 110 j=j1,m
110         rowi(j)=-a(j,i)
c
c  do 130 accomplishes the same result as call mxmup
c           do 130 jj=j1,m
c           do 130 ii=1,nrows
c30         if ( jj+1-ii .gt. 0 )
c    *      a(m+2-ii,i+ii-1)=a(m+2-ii,i+ii-1)
c    *      +rowi(jj)*a(jj+1-ii,i-1+ii)
         call mxmup( a(1,i),lda-1 ,rowi,1 ,a(m+1,i),lda-1 ,m )
c
      info=i
      if(a(m+1,i).le.0. ) go to 60
         a(m+1,i)=1./sqrt(a(m+1,i))
            do 140 ii=1,nrows-1
140         a(m+1-ii,i+ii)=a(m+1-ii,i+ii)*a(m+1,i)
c
150      continue
c
c  again do 250 could be accomplished by matrix times vector operations
c  where we have both rectangular and triangular matrices
c
         do 250 i=n+1-m,n
c
c        --- get rowi ---
c
         j1=max0(m+2-i,1)
            do 210 j=j1,m
210         rowi(j)=-a(j,i)
c
         nrows=min0(m+1,n+1-i)
            do 230 jj=j1,m
cdir$ ivdep
            do 230 ii=1,min0(jj,nrows)   ] nrows
230      ]  if ( jj+1-ii .gt. 0 )
     *      a(m+2-ii,i+ii-1)=a(m+2-ii,i+ii-1)
     *      +rowi(jj)*a(jj+1-ii,i-1+ii)
c
      info=i
      if(a(m+1,i).le.0. ) go to 60
         a(m+1,i)=1./sqrt(a(m+1,i))
            do 240 ii=1,nrows-1
240         a(m+1-ii,i+ii)=a(m+1-ii,i+ii)*a(m+1,i)
c
250      continue
      info=0
60    return
      end
c
c   --- a fortran version (no attempt at efficiency) of the kernal
      subroutine mxmupf( a,lda ,r,ldr ,y,ldy ,m )
      real a(lda,m),r(ldr,m),y(ldy,m)
c
         do 10 i=1,m
         do 10 j=1,m
10       if ( j.ge.i ) y(1,i)=y(1,i)+ r(1,j)*a(j,i)
c
      return
      end
c
      subroutine bndchls( a,lda ,n,m ,x,b )
      real a(lda,n) ,x(n),b(n)
c
c   solve ax=b : where a is cholesky decomposition (output from bndchl)
c   valid for the case where x and b occupy the same memory storage
c   as well as the case where x ,and b as disjoint (in memory)
c
      if ( loc(x).ne.loc(b) ) then
         do 1 j=1,n
1        x(j)=b(j)
      endif
c
         do 20 k=1,n
         x(k)=x(k)*a(m+1,k)
cdir$ ivdep
            do 10 i=k+1,min0(k+m,n)
10          x(i)=x(i)-x(k)*a(m+k+1-i,i)
20       continue
c
         do 40 k=n,1,-1
         x(k)=x(k)*a(m+1,k)
cdir$ ivdep
            do 30 i=1,min0(k-1,m)
30          x(k-i)=x(k-i)-x(k)*a(m+1-i,k)
40       continue
c
      return
      end

.