program testudu
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      parameter(neq=99,ib=2)
      dimension a(neq,ib), b(neq)

      do j=1,ib
         do i=1,neq
            a(i,j)=0.0
	end do
      end do	

c ... setup matrix
      do i=1,neq
        a(i,1)= 2.0
        a(i,2)=-1.0
        b(i)=0.0
      end do
      a(neq,2)=0.0
      dx=1.0/float(neq+1)
      b(1)=1.0
      b(neq)=0.0
      call udu(a,b,1,neq,ib,neq,ib,6)
c     call udu(a,b,0,neq,ib,neq,ib,6)
      x=0.0
      do i=1,neq
        x=x+dx
        write(*,1000) i,x,b(i)
      end do
      stop
 1000 format('i,x,b(i)=',i4,1x,e10.4,1x,e10.4)
      end