c	file inv16.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),av(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
	print *, 'unit number for output?'
	read(5,99)iout
	do 300 jc=n,1,-1
	  nc=n-jc+1
	  kc=(n-jc)/2
	  do 10 j=1,kc
	   do 11 i=1,nc
	    c(i,j)=glam(n4+1-i,j+jc)
11	   continue
10	  continue
	  do 12 j=1,nc-kc
	   do 13 i=1,nc
	    c(i,kc+j)=glam(n4+1-i,n+kc+j)
13	   continue
12	  continue
	do 14 i=1,n
	 alph(i)=-glam(n4+1-i,jc)
14	continue
	itask=1
	call dgefs(c,lda,nc,alph,itask,ind,work,iwork)

	av(jc)=1
	write(11,*)av(jc)
	do 17 i=1,kc
	  av(jc+i)=alph(i)
	write(11,*)av(jc+i)
17	continue
	do 18 i=jc+kc+1,n
	  av(i)=0
	write(11,*)av(i)
18	continue
	write(11,*)
	do 21 i=jc,n
	  a=glam(i,jc)
	  do 19 j=1,kc
	    a=a+glam(i,j+jc)*alph(j)
19	  continue
	  do 20 j=1,nc-kc
	    a=a+glam(i,n+kc+j)*alph(kc+j)
20	  continue
	  write(10,*)a
	  xi(i)=a
21	continue
	  write(10,*)
	do 100 j=jc+1,n
	  v(1,j)=av(j)-xi(j)/gv(1,j)
	  write(13,*)v(1,j)
100	continue
	  write(13,*)
	do 102 j=jc+2,n-1
	  v(2,j)=v(1,j)+(gv(1,j)*(v(1,j)-av(j))+
     1	  gh(1,j)*(v(1,j)-v(1,j-1))+
     2	  gh(1,j+1)*(v(1,j)-v(1,j+1)))/gv(2,j)
	  write(13,*)v(2,j)
102	continue
	  write(13,*)
	do 103 i=3,nc/2
	  do 104 j=jc+i,n+1-i
	   v(i,j)=v(i-1,j)-(gv(i-1,j)*(v(i-2,j)-v(i-1,j))+
     1	   gh(i-1,j)*(v(i-1,j-1)-v(i-1,j))+ 
     2	   gh(i-1,j+1)*(v(i-1,j+1)-v(i-1,j)))/gv(i,j)
	  write(13,*)v(i,j)
104	  continue
	  write(13,*)
103	continue
	
	gv(1,jc)=xi(jc)
	write(iout,*)gv(1,jc)
	ai=xi(jc)
	if(kc .ge. 1) then
	  gh(1,jc+1)=-ai/v(1,jc+1)
	  write(iout,*)gh(1,jc+1)
	  ai=gv(1,jc+1)*(av(jc+1)-v(1,jc+1))
     1	  -gh(1,jc+1)*v(1,jc+1)+
     2	  gh(1,jc+2)*(v(1,jc+2)-v(1,jc+1))
	  gv(2,jc+1)=ai/v(1,jc+1)
	  write(iout,*)gv(2,jc+1)
	endif
	do 200 j=2,kc
	  gh(j,jc+j)=-ai/v(j,jc+j)
	  write(iout,*)gh(j,jc+j)
	  ai=gv(j,jc+j)*(v(j-1,jc+j)-v(j,jc+j))
     1	  -gh(j,jc+j)*v(j,jc+j)+
     2	  gh(j,jc+j+1)*(v(j,jc+j+1)-v(j,jc+j))
	  gv(j+1,jc+j)=ai/v(j,jc+j)
	  write(iout,*)gv(j+1,jc+j)
200	continue
	if(mod(n-jc,2) .eq. 1) then
	  gh(kc+1,jc+kc+1)=-ai/v(kc+1,jc+kc+1)
	  write(iout,*)gh(kc+1,jc+kc+1)
	endif
	write(iout,*)
300	continue
	tm=dtime(tarray)
	print *, 'compute and write time  =  ',tm
	stop
	end
