c	file cdton.f
c	computes the complex dirichlet to neumann map for
c	an n by n resistor network
c
c	variables: zgv=vertical conductivities
c		zgh=horizontal conductivities
c		za=matrix for solving the Dirichlet problem
c		zabe=banded version of za sent to cnbco
c		zxb=right hand side of equation for solving
c			the Dirchlet problem including the boundary
c			conditions
c		zglam=Dirichlet to Neumann map

	implicit complex (z)
	dimension b(4,20)
	dimension zgv(21,20),zgh(20,21)
	dimension zxb(400),zabe(400,71),za(400,400)
	dimension z(400)
	dimension zglam(80,80),tarray(2)
	dimension ipvt(400)

	lda=400
	read(1,*)n
	np1=n+1
	nm1=n-1
	nsq=n*n
	n2=2*n
	n3=3*n
	n4=4*n
	do 10 i=1,np1
	  read(1,*)(zgv(i,j),j=1,n)
10	continue
	do 11 i=1,n
	  read(1,*)(zgh(i,j),j=1,np1)
11	continue
	do 16 i=1,nsq
	  j=(i-1)/n +1
	  k=mod(i-1,n)+1
	  za(i,i)=zgh(j,k)+zgh(j,k+1)+zgv(j,k)+zgv(j+1,k)
16	continue
	jj=1
	do 30 ii=1,n
	 do 17 i=jj,jj+n-2
	  j=(i-1)/n +1
	  k=mod(i-1,n)+1
	  za(i,i+1)=-zgh(j,k+1)
	  za(i+1,i)=za(i,i+1)
17	 continue
	 jj=jj+n
30	continue
	do 18 i=1,nsq-n
	  j=(i-1)/n +1
	  k=mod(i-1,n)+1
	  za(i,i+n)=-zgv(j+1,k)
	  za(i+n,i)=za(i,i+n)
18	continue
	
	ml=n
	mu=n
	nc=2*ml+mu+1
	do 21 i=1,nsq
	  j1=max0(1,i-ml)
	  j2=min0(nsq,i+mu)
	  do 20 j=j1,j2
	    k=j-i+ml+1
	    zabe(i,k)=za(i,j)
20	  continue
21	continue
	do 102 i=1,np1
	  write(2,999)(real(zgv(i,j)),aimag(zgv(i,j)),j=1,n)
102	continue
	write(2,*)
	do 110 i=1,n
	  write(2,999)(real(zgh(i,j)),aimag(zgh(i,j)),j=1,np1)
110	continue
	write(2,*)'beginning of za matrix'
	do 104 i=1,nsq
	  write(2,999)(real(za(i,j)),aimag(za(i,j)),j=1,nsq)
	  write(2,*)
104	continue
	write(2,*)'beginning of zabe matrix'
	do 105 i=1,nsq
	  write(2,*)(real(zabe(i,j)),j=1,nc)
	  write(2,*)
105	continue
999	format(2f10.5)
	tm=dtime(tarray)
	print *,'read, write, and, set up time = ',tm
	write(2,*)
	call cnbco(zabe,lda,nsq,ml,mu,ipvt,rcond,z)
	write(2,*)'rcond=',rcond
	write(2,*)'new zabe'
	do 106 i=1,np1
	  write(2,999)(real(zabe(i,j)), aimag(zabe(i,j)),j=1,nc)
106	continue

c	beginning of main loop to compute zglam
	do 200 kk=1,4*n
	do 12 i=1,4
	do 40 j=1,n
	  b(i,j)=0
	  ii=(kk-1)/n+1
	  jj=mod(kk-1,n)+1
	  if(ii .eq. i .and. jj .eq. j)b(i,j)=1
40	continue
12	continue
	do 41 i=1,nsq
	  zxb(i)=0
41	continue
	zxb(1)=zgv(1,1)*b(1,1)+zgh(1,1)*b(4,n)
	zxb(n)=zgv(1,n)*b(1,n)+zgh(1,np1)*b(2,1)
	zxb(nsq-n+1)=zgh(n,1)*b(4,1)+zgv(np1,1)*b(3,n)
	zxb(nsq)=zgv(np1,n)*b(3,1)+zgh(n,np1)*b(2,n)
	do 13 i=2,nm1
	  zxb(i)=zgv(1,i)*b(1,i)
13	continue
	do 14 i=1,n-2
	  zxb(i*n+1)=zgh(i+1,1)*b(4,n-i)
	  zxb(i*n+n)=zgh(i+1,np1)*b(2,i+1)
14	continue
	nn=nsq-n
	do 15 i=2,nm1
	  zxb(nsq-n+i)=zgv(np1,i)*b(3,np1-i)
15	continue

	write(2,*)
	call cnbsl(zabe,lda,nsq,ml,mu,ipvt,zxb,job)

	do 201 k=1,n
	  zglam(kk,k)=(b(1,k)-zxb(k))*zgv(1,k)
201	continue
	do 202 k=1,n
	  zglam(kk,n+k)=(b(2,k)-zxb(n*k))*zgh(k,np1)
202	continue
	do 203 k=1,n
	  zglam(kk,n2+k)=(b(3,k)-zxb(nsq-k+1))*zgv(np1,np1-k)
203	continue
	do 204 k=1,n
	  zglam(kk,n3+k)=(b(4,k)-zxb(n*(n-k)+1))*zgh(np1-k,1)
204	continue
c	write(2,*)(zglam(kk,k),k=1,n4)
	do 220 k=1,n4
	  write(3,*)zglam(kk,k)
220	continue
200	continue
	tm=dtime(tarray)
	print *,'compute and write time = ',tm
	stop
	end
