c	file inv10.f
c	Solves the inverse resistor problem.
c	Takes lambda, produces gamma.

	implicit double precision (a-h,o-z)
	dimension glam(80,80),v(20,20),gh(20,21),gv(21,20)
	dimension alph(20),c(20,20),work(20),iwork(20)
	dimension xi(20)
	dimension tarray(2)

	lda=20
	read(3,99)n
99	format(i2)
	n2=2*n
	n3=3*n
	n4=4*n
	no2=n/2
	nm12=(n-1)/2
	np12=(n+1)/2
	do 5 j=1,n4
	  do 6 i=1,n4
	    read(3,*)glam(i,j)
6	  continue
5	continue
	tm=dtime(tarray)
	print *, 'read time  =  ',tm
	do 7 j=1,n
	  do 8 i=1,n+1
	    gv(i,j)=1
8	  continue
7	continue
	do 71 j=1,n+1
	  do 81 i=1,n
	    gh(i,j)=1
81	  continue
71	continue
	do 10 j=1,n
	  do 11 i=1,n
	    c(i,j)=glam(n4-n+i,n+j)
	    write(8,*)c(i,j)
11	  continue
10	continue
	do 12 i=1,n
	 alph(i)=-glam(n4-n+i,n-n+1)
12	continue
	itask=1
	ii=n
	call dgefs(c,lda,ii,alph,itask,ind,work,iwork)
	write(9,*)ind
	do 15 i=1,n
	  write(9,*)alph(i)
15	continue
	write(9,*)
	do 20 i=1,n
	  a=glam(n-n+i,n-n+1)
	  do 19 j=1,n
	    a=a+glam(n-n+i,n+j)*alph(j)
19	  continue
	  write(10,*)a
	  xi(i)=a
20	continue
	do 50 i=1,n
	  v(i,i)=0
50	continue
	v(1,2)=-gv(1,1)/gh(1,2)
	write(12,*)v(1,2)
	write(12,*)
c	v(1,3)=v(1,2)+v(1,2)*(gv(1,2)+gh(1,2)+gv(2,2))/gh(1,3)
c	v(2,3)=-gv(2,2)*v(1,2)/gh(2,3)
c	write(12,*)
	do 102 j=2,np12-1
	  v(1,j+1)=v(1,j)+(gv(1,j)*v(1,j)+gv(2,j)*(v(1,j)-v(2,j))+
     1	  gh(1,j)*(v(1,j)-v(1,j-1)))/gh(1,j+1)
	  write(12,*)v(1,j+1)
	  do 101 i=2,j-1
	    v(i,j+1)=v(i,j)+(gv(i,j)*(v(i,j)-v(i-1,j))+
     1	    gv(i+1,j)*(v(i,j)-v(i+1,j))+
     2	    gh(i,j)*(v(i,j)-v(i,j-1)))/gh(i,j+1)
	    write(12,*)v(i,j+1)
101	  continue
	  v(j,j+1)=-v(j-1,j)*gv(j,j)/gh(j,j+1)
	  write(12,*)v(j,j+1)
	  write(12,*)
102	continue
	do 104 j=np12,n-2
	  v(1,j+1)=v(1,j)+(gv(1,j)*v(1,j)+gv(2,j)*(v(1,j)-v(2,j))+
     1	  gh(1,j)*(v(1,j)-v(1,j-1)))/gh(1,j+1)
	  write(12,*)v(1,j+1)
	  do 103 i=2,n-j-1
	    v(i,j+1)=v(i,j)+(gv(i,j)*(v(i,j)-v(i-1,j))+
     1	    gv(i+1,j)*(v(i,j)-v(i+1,j))+
     2	    gh(i,j)*(v(i,j)-v(i,j-1)))/gh(i,j+1)
	    write(12,*)v(i,j+1)
103	  continue
	  write(12,*)
104	continue
	v(1,n)=-xi(n)/gv(1,n)
	write(12,*)v(1,n)
	v(2,n-1)=v(1,n-1)+(gh(1,n-1)*(v(1,n-1)-v(1,n-2))+
     1	gh(1,n)*(v(1,n-1)-v(1,n))+
     2	gv(1,n-1)*v(1,n-1))/gv(2,n-1)
	write(12,*)v(2,n-1)
	do 110 i=3,no2
	  v(i,n+1-i)=v(i-1,n+1-i)+(gh(i-1,n+1-i)*(v(i-1,n+1-i)-
     1	  v(i-1,n-i))+gh(i-1,n+2-i)*(v(i-1,n+1-i)-v(i-1,n+2-i))+
     2	  gv(i-1,n+1-i)*(v(i-1,n+1-i)-v(i-2,n+1-i)))/gv(i,n+1-i)
	write(12,*)v(i,n+1-i)
110	continue
	tm=dtime(tarray)
	print *, 'compute and write time  =  ',tm
	stop
	end
