c	file dirp3.f
c	computes the dirichlet to neumann map for
c	an n by n resistor network

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

	read(1,*)n
	np1=n+1
	nm1=n-1
	nsq=n*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 12 i=1,4
	  read(1,*)(b(i,j),j=1,n)
12	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
	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

	ml=n
	mu=n
	m=2*n+1
	lda=3*n+1
	do 20 j=1,nsq
	  i1=max(1,j-mu)
	  i2=min(nsq,j+ml)
	  do 19 i=i1,i2
	    k=i-j+m
	    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 103 i=1,4
	  write(2,*)(b(i,j),j=1,n)
103	continue
	write(2,*)
	do 104 i=1,nsq
	  write(2,*)(a(i,j),j=1,nsq)
104	continue
	write(2,*)
	do 105 i=1,lda
	  write(2,*)(abd(i,j),j=1,nsq)
105	continue
	write(2,*)
	write(2,*)(xb(i),i=1,nsq)
	write(2,*)
	call dgbfa(abd,lda,nsq,ml,mu,ipvt,info)
	write(2,*)info
	write(2,*)
	do 106 i=1,lda
	  write(2,*)(abd(i,j),j=1,nsq)
106	continue
c	call dpbsl(abd,lda,nsq,m,xb)
	write(2,*)(xb(i),i=1,nsq)

	stop
	end
