c	file dton2.f
c	computes the dirichlet to neumann map for
c	an n by n resistor network
c	Rounds the answer to two digits.

	implicit double precision (a-h,o-y)
	dimension b(4,20),gv(21,20),gh(20,21)
	dimension xb(400),abd(21,400),a(400,400)
	dimension glam(80,80),tarray(2)

	lda=21
	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,*)(gv(i,j),j=1,n)
10	continue
	do 11 i=1,n
	  read(1,*)(gh(i,j),j=1,np1)
11	continue
	do 16 i=1,nsq
	  j=(i-1)/n +1
	  k=mod(i-1,n)+1
	  a(i,i)=gh(j,k)+gh(j,k+1)+gv(j,k)+gv(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
	  a(i,i+1)=-gh(j,k+1)
	  a(i+1,i)=a(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
	  a(i,i+n)=-gv(j+1,k)
	  a(i+n,i)=a(i,i+n)
18	continue
	m=n
	do 20 j=1,nsq
	  i1=max(1,j-m)
	  do 19 i=i1,j
	    k=i-j+m+1
	    abd(k,j)=a(i,j)
19	  continue
20	continue
	do 102 i=1,np1
	  write(2,*)(gv(i,j),j=1,n)
102	continue
	write(2,*)
	do 110 i=1,n
	  write(2,*)(gh(i,j),j=1,np1)
110	continue
	write(2,*)
	do 104 i=1,nsq
	  write(2,*)(a(i,j),j=1,nsq)
104	continue
	write(2,*)
	do 105 i=1,np1
	  write(2,*)(abd(i,j),j=1,nsq)
105	continue
	tm=dtime(tarray)
	print *,'read, write, and, set up time = ',tm
	write(2,*)
	call dpbfa(abd,lda,nsq,m,info)
	write(2,*)'info=',info
	write(2,*)
	do 106 i=1,np1
	  write(2,*)(abd(i,j),j=1,nsq)
106	continue

	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
	  xb(i)=0
41	continue
	xb(1)=gv(1,1)*b(1,1)+gh(1,1)*b(4,n)
	xb(n)=gv(1,n)*b(1,n)+gh(1,np1)*b(2,1)
	xb(nsq-n+1)=gh(n,1)*b(4,1)+gv(np1,1)*b(3,n)
	xb(nsq)=gv(np1,n)*b(3,1)+gh(n,np1)*b(2,n)
	do 13 i=2,nm1
	  xb(i)=gv(1,i)*b(1,i)
13	continue
	do 14 i=1,n-2
	  xb(i*n+1)=gh(i+1,1)*b(4,n-i)
	  xb(i*n+n)=gh(i+1,np1)*b(2,i+1)
14	continue
	nn=nsq-n
	do 15 i=2,nm1
	  xb(nsq-n+i)=gv(np1,i)*b(3,np1-i)
15	continue

	write(2,*)
	call dpbsl(abd,lda,nsq,m,xb)

	do 201 k=1,n
	  glam(kk,k)=(b(1,k)-xb(k))*gv(1,k)
201	continue
	do 202 k=1,n
	  glam(kk,n+k)=(b(2,k)-xb(n*k))*gh(k,np1)
202	continue
	do 203 k=1,n
	  glam(kk,n2+k)=(b(3,k)-xb(nsq-k+1))*gv(np1,np1-k)
203	continue
	do 204 k=1,n
	  glam(kk,n3+k)=(b(4,k)-xb(n*(n-k)+1))*gh(np1-k,1)
204	continue
	do 220 k=1,n4
	  write(3,999)glam(kk,k)
220	continue
200	continue
999	format(e10.3)
	tm=dtime(tarray)
	print *,'compute and write time = ',tm
	stop
	end
