subroutine getrow( a,lda ,n,m ,i ,row ,l1,l2)
      real a(lda,n) ,row(n)

c     --- gets the ith row of a symmetric banded matrix ---

c      if ( numarg().eq.8 ) then
         k1=l1
         k2=l2
c      else
c         k1=1
c         k2=n
c      endif

         do 7 j=k1,k2
7        row(j)=0.

         do 2 j=0,min0(m,i-1)  
2        row(i-j)=a(m+1-j,i)


         do 3 j=1,min0(m,n-i)  
3        row(i+j)=a(m+1-j,i+j)

      return
      end

      subroutine putval( a,lda ,n,m ,i,j ,val )
      real a(lda,n)

c     --- stores an element into a symmetric band store matrix ---

**?   if ( iabs(i-j).gt.m ) return
**?   if ( i.lt.1 .or. j.lt.1 .or. i.gt.n .or. j.gt.n ) return
      ii=max0(i,j)
c     jj=iabs(i-j)
      a(m+1-iabs(i-j),ii)=val
      return
      end

      subroutine symmsq( a,lda,n,m ,as,ldas,ms )
      real a(lda,n) ,as(ldas,n)
     *  ,r(10000),c(10000)
      common /qqq/ aa(200,200),bb(200,200)

c   --- squares the symmetric matrix a ---

      ms=2*m
         write(6,*)'debug: symmsq: do 1'
         do 1 j=1,ms+1
         do 1 i=1,n
1        as(j,i)=0.

         write(6,*)'debug: symmsq: do 10'
         do 10 i=1,n
         call getrow( a,lda,n,m ,i,r )
      k1=max(1,i-m-1)
      k2=min(n,i+m+1)

         do 10 j=max0(1,i-ms),i
         call getrow( a,lda,n,m ,j,c ,k1,k2)

         s=0.
         write(6,*)'debug: symmsq: do 5'
            do 5 k=k1,min(k2,j+m+1)
cccc        do 5 k=k1,k2
5           s=s+r(k)*c(k)

         call putval( as,ldas,n,ms ,i,j ,s )
10       continue

      return
      end
c   file bandc forfor

      subroutine bandc ( a,lda ,n,m ,info )
      real a(lda,n)
*dir$ vfunction sxmups ]fortran only version.
      real rowi(1000)

c     --- matrix decomposition for banded positive definite matrices
c     --- ( this routine, and it's partner bndchls can be used to
c     --- accomplish the same thing as linpack's spbfa,spbsl -note
c     --- however that bandc&bndchls store d**(-1/2) on the diagnol
c     --- whereas linpack stores d**(1/2)
c     --- note: data is most frequently accessed stride (lda-1)
c     --- avoid lda = l*2**q+1 : q >2

         do 50 i=1,m

c        --- get rowi ---

         j1=max0(m+2-i,1)
            do 10 j=j1,m
10          rowi(j)=-a(j,i)

         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)

      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)

50       continue

      j1=1
      nrows=     m+1

      if ( m.le.64 ) then
         do 100 i=m+1,n-m
        a(m+1,i)=sxmupf( a(1,i),a(1,i),a(m+1,i),lda-1,m )
c         a(m+1,i)=sxmups( loc(a(1,i)),loc(a(1,i)),loc(a(m+1,i))
c     *                    ,lda-1,m )
         info=i
         if ( a(m+1,i).le.0. ) goto 60
         a(m+1,i)=1./a(m+1,i)
cdir$ ivdep
            do 100 ii=1,m
100         a(m+1-ii,i+ii)=a(m+1-ii,i+ii)*a(m+1,i)
      else
         do 150 i=m+1,n-m
c        --- get rowi ---

            do 110 j=j1,m
110         rowi(j)=-a(j,i)

c            do 130 jj=j1,m
c            do 130 ii=1,nrows
c130         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 )

      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)

150      continue
      endif

         do 250 i=n+1-m,n

c        --- get rowi ---

         j1=max0(m+2-i,1)
            do 210 j=j1,m
210         rowi(j)=-a(j,i)

         nrows=min0(m+1,n+1-i)
            do 230 jj=j1,m
cdir$ ivdep
            do 230 ii=1,min0(jj,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)

      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)

250      continue
      info=0
60    return
      end
      subroutine bndchls( a,lda ,n,m ,x,b )
      real a(lda,n) ,x(n),b(n)

      if ( loc(x).ne.loc(b) ) then
         do 1 j=1,n
1        x(j)=b(j)
      endif

         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

         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

      return
      end

c /// fortran versions of assembly language kernals
c      function sxmups ( pq,pr,py,ldy ,m )
c      integer pq,pr,py
c      real p(*),r(*),q(*)
c      pointer ( pq,q ),( pr,r ),( py,y )
c      sxmups = sxmupf ( q,r,y,ldy ,m )
c      return
c      end
      function sxmupf( q,r,y,ldy ,m )
      real q(ldy,m),r(m),y(ldy,m)
c /// specialized m in range (1,64) only (incorrect results otherwise
         do 10 j=m,1,-1
         do 10 i=1,j
10       y(1,i)=y(1,i) -r(j)*q(j,i)
      sxmupf=y(1,1)
      if ( sxmupf.gt.0 ) sxmupf=sqrt(sxmupf)
      return
      end
      subroutine mxmup ( a,lda ,r,ldr ,y,ldy ,m )
      real a(lda,m),r(ldr,m),y(ldy,m)

         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)

      return
      end
      function ranf()
      real ranf
c
      data init /1325/
            init = mod(3125*init,65536)
            ranf = (init - 32768.0)/16384.0
      return
      end