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