subroutine arnoldi(a, ia, ja, lda, na, b, n, work, iwork)
      parameter(m = 200, nn = 3000)
      implicit double precision (a-h,o-z)
      double precision  a(lda), b(n), work(1), r(nn), h(m+1,m+1),
     c                  v(nn,m+1), av(nn), vy(nn), vv(m), w(nn)
      integer ja(lda), ia(lda), iwork(1)
c
      maxit = 10
      crit  = 0.000000001d0
      do 10 i=1,n
        r(i)  = b(i)
        b(i)  = 0.d0
 10   continue
      do 20 i=1,maxit
        do 30 j=1,m+1
          do 30 k=1,m+1
            h(j,k)=0.d0
 30     continue
        sum1 = 0.d0
        do 40 j=1,n
          sum1 = sum1+r(j)*r(j)
 40     continue
        sum1 = dsqrt(sum1)
        do 50 j=1,n
          v(j,1)=r(j)/sum1
 50     continue
        do 60 j=1,m
          do 70 k=1,n
            av(k) = 0.d0
 70       continue 
          do 80 k=1,na
            av(ia(k)) = av(ia(k))+ a(k)*v(ja(k),j)
 80       continue
          do 90 k=1,n
            w(k)  =0.d0
 90       continue
          do 100 k=1,j
            sum2 = 0.d0
            do 110 l=1,n
              sum2 = sum2+av(l)*v(l,k)
 110        continue
            h(k,j) = sum2
            do 120 l=1,n
              w(l) = w(l)+h(k,j)*v(l,k)
 120        continue
 100      continue
          do 130 k=1,n
            w(k) = av(k)-w(k)
 130      continue
          sum2 = 0.d0
          do 140 k=1,n
            sum2 = sum2+w(k)*w(k)
 140      continue
          h(j+1,j) = dsqrt(sum2)
          do 150 k=1,n
            v(k,j+1) = w(k)/h(j+1,j)
 150      continue
 60     continue
        do 160 k=1,m
          vv(k) = 0.d0
 160    continue
        vv(1) = sum1
        itask = 1
        call dgefs (h,m+1,m,vv,itask,ind1,work,iwork,rcond)
        do 180 k=1,n
          vy(k) = 0.d0
          do 180 l=1,m
            vy(k) = vy(k)+v(k,l)*vv(l)
 180    continue
        do 190 k=1,n
          b(k)=b(k)+vy(k)
 190    continue
        do 200 k=1,na
          r(ia(k)) = r(ia(k)) - a(k)*vy(ja(k))
 200    continue
        sum2 = 0.d0
        do 210 k=1,n 
          sum2 = sum2 + r(k)*r(k)
 210    continue
        write(*,*) sum2
        if (dsqrt(sum2)/float(n).lt.crit) go to 999
 20   continue
 999  continue
      return
      end