subroutine eiginv(a,valu,p,v,resid,n,eps,limit,delmax)
c
c     numerical solution of the inverse eigenvalue problem
c
      double precision a(20,20,23),valu(20),p(20),v(20),resid(20)
      double precision eps,delmax,temp
      integer n,limit,i,j,k,l,ierr,it,info,ipvt(20)
c
      delmax = eps
      it = 0
c
c     determine eigensystem of next approximation (with values ordered)
c
   50 do 100 i = 1, n
         do 90 j = 1, i
            a(i,j,n+2) = -a(i,j,n+1)
   90       continue
  100    continue
      do 200 i = 1, n
         do 190 j = 1, i
            do 180 k = 1, n
               a(i,j,n+2) = a(i,j,n+2) - p(k)*a(i,j,k)
  180          continue
  190       continue
  200    continue
c
      call rs(20,n,a(1,1,n+2),v,1,a(1,1,n+2),a(1,1,n+3),a(1,1,n+3),ierr)
      do 210 i = 1, n
         v(i) = -v(i)
  210    continue
      if (delmax .lt. eps .or. it .eq. limit) return
c
c     form and solve linear system derived in next newton iteration
c
      do 230 i = 1, n
         resid(i) = valu(i) - v(i)
  230    continue
c
      do 300 i = 1, n
         do 290 k = 1, n
            a(i,k,n+3) = 0.0
            do 250 j = 1, n
               v(j) = 0.0
               do 240 l = 1, n
                  v(j) = v(j) + a(j,l,k)*a(l,i,n+2)
  240             continue
  250          continue
            do 280 l = 1, n
               a(i,k,n+3) = a(i,k,n+3) + v(l)*a(l,i,n+2)
  280          continue
  290       continue
  300    continue
c
         call dgefa(a(1,1,n+3),20,n,ipvt,info)
         call dgesl(a(1,1,n+3),20,n,ipvt,resid,0)
c
      delmax = 0.0
      do 400 k = 1, n
         p(k) = p(k) + resid(k)
         temp = abs(resid(k))
         if (temp .gt. delmax) delmax = temp
  400    continue
c
      it = it + 1
      go to 50
c
      end