C ALGORITHM 732, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 20, NO. 3, SEPTEMBER, 1994, PP. 247-261.
cat > Readme < Makefile < testcapc.out
testdcapc > testdcapc.out
testcapcn > testcapcn.out
testdcapcn > testdcapcn.out
testreccn > testreccn.out
testdreccn > testdreccn.out
# testcapc
OBJECTS1= testcapc.o solvers.o poiss.o linpack.o
testcapc: $(OBJECTS1)
$(FC) -o testcapc $(OBJECTS1)
# testdcapc
OBJECTS2= testdcapc.o solvers.o poiss.o linpack.o
testdcapc: $(OBJECTS2)
$(FC) -o testdcapc $(OBJECTS2)
# testcapcn
OBJECTS3= testcapcn.o solvers.o poiss.o linpack.o
testcapcn: $(OBJECTS3)
$(FC) -o testcapcn $(OBJECTS3)
# testdcapcn
OBJECTS4= testdcapcn.o solvers.o poiss.o linpack.o
testdcapcn: $(OBJECTS4)
$(FC) -o testdcapcn $(OBJECTS4)
# testreccn
OBJECTS5= testreccn.o solvers.o poiss.o linpack.o
testreccn: $(OBJECTS5)
$(FC) -o testreccn $(OBJECTS5)
# testdreccn
OBJECTS6= testdreccn.o solvers.o poiss.o linpack.o
testdreccn: $(OBJECTS6)
$(FC) -o testdreccn $(OBJECTS6)
clean :
rm -f core $(ALL) *.o
C-END-OF-FILE
cat > linpack.f < poiss.f < solvers.f < 0 must hold.
c
c DY real constant
c The discretization interval in the Y direction.
c The condition DY > 0 must hold.
c
c BA real array (NY)
c The boundary values on the regular boundary
c at (I=1, J=2,NY-1).
c
c BB real array (NY)
c The boundary values on the regular boundary
c at (I=NX, J=2,NY-1).
c
c BC real array (nx)
c The boundary values on the regular boundary
c at (I=1,NX, J=1).
c
c BD real array (NX)
c The boundary values on the regular boundary
c at (I=1,NX, J=NY).
c
c VP real array (NP)
c The boundary values on the irregular boundary,
c defined by the (IP,JP) coordinates.
c
c IP integer array (NP)
c The coordinates in the X direction of the
c irregular boundary points. The condition,
c 1 < IP(N) < NX, must hold for all N.
c
c JP integer array (NP)
c The coordinates in the Y direction of the
c irregular boundary points. The condition,
c 1 < JP(N) < NY, must hold for all N.
c
c Additional notes on defining the IP and JP arrays are
c given below.
c
c Note: The order in which the irregular boundary grid points, defined
c by the (IP(N),JP(N)) pairs, are specified is immaterial.
c However, these grid points must all be distinct; there must
c not be a repeated (IP(N),JP(N)) pair. Failure to ensure this
c will result in the error condition IERROR=4, defined below.
c
c ICFLAG integer constant
c A flag: ICFLAG=0 - do preprocessing only and compute
c the matrix C and indices IPS.
c
c ICFLAG=1 - do preprocessing to compute the
c matrix C and indices IPS, and solve
c the Helmholtz problem.
c
c ICFLAG=2 - solve the Helmholtz problem only; this
c assumes that preprocessing has been
c done and that C and IPS are defined.
c
c Any other value for ICFLAG gives IERROR=5 and no attempt
c is made to do preprocessing or to solve. On return,
c ICFLAG=2.
c
c
c ******************* ON RETURN ********************
c
c U real array (NX,NY)
c The solution in the irregular domain, with
c U(IP(N),JP(N)) = VP(N) for N=1,NP on the irregular
c portion of the boundary and also with
c U(1,J)=BA(J) for J=2,NY-1
c U(NX,J)=BB(J) for J=2,NY-1
c U(I,1)=BC(I) for I=1,NX
c U(I,NY)=BD(I) for I=1,NX.
c
c C real array (NP,NP)
c The inverse capacitance matrix. This array must not
c be altered on successive calls to the routine with
c different source functions when ICFLAG=2.
c
c IPS integer array (NP)
c An integer vector of pivot indices. This array must
c not be alter on successive calls with different
c source functions when ICFLAG=2.
c
c TRIGS real array (2*(NX-1))
c An array of coefficients required by the FFT
c routine employed by the fast solver. This array
c is initialized on the first call to CAPC, or
c whenever ICFLAG=0 or 1. The array must not be
c altered between successive calls to the routine
c with different source functions when ICFLAG=2.
c
c IERROR integer constant
c An error flag on input variables. On return IERROR
c has the following meaning:
c
c IERROR=0 - normal status, no error detected.
c =1 - value of NX-1 or NY-1 is less than 8.
c =2 - illegal value for NX.
c =3 - illegal value for IP or JP.
c =4 - an error occurred in factoring the inverse
c capacitance matrix, i.e. the matrix is singular.
c This condition will occur if the user has
c defined duplicate irregular boundary points.
c =5 - illegal value for ICFLAG.
c =6 - illegal value for DX or DY.
c
c Note: Error checking is normally only done on the first call
c to the routine. The user may force the routine to do error
c checking at any time by setting IERROR=-1 on entry.
c
c
c ********* NOTES ON CALLING ROUTINE CAPC ********
c
c In applying the capacitance matrix method to the solution
c of the Helmholtz equation, routine CAPC works in two
c stages. The first stage is a preprocessing one in which
c the arrays C and IPS are determined. The second stage is
c the actual solution of the problem for a given source
c function. The flag ICFLAG determines the operation of the
c routine:
c
c If the routine is called with ICFLAG=0, only
c preprocessing is done, the arrays C and IPC
c are determined, and CAPC immediately returns to
c the calling program without solving. On return
c ICFLAG=2
c
c If the routine is called with ICFLAG=1,
c preprocessing is done, the arrays C and IPC
c are determined, and the solution for the
c given source function is determined. On return
c ICFLAG=2.
c
c If the routine is called with ICFLAG=2,
c the solution for the given source function
c is determined. This assumes that preprocessing
c has been done and that the C and IPS arrays
c are defined. On return ICFLAG=2.
c
c If the user wishes to solve the Helmholtz equation in a
c given domain repeatedly with different source functions,
c then the routine could be called with ICFLAG=1 initially.
c Afterwards this flag is set to 2 and need not be altered.
c If the user has stored the arrays C and IPS from a previous
c run, then these may be supplied as input data and the
c routine may be called initially with ICFLAG=2 to obtain
c the solution directly.
c
c Preprocessing need only be done once for a given domain
c geometry and Helmholtz coefficient. The user-specified
c variables listed below must not change between successive
c calls to CAPC with ICFLAG=2. If the user changes any
c one of the following variables:
c
c NX,
c NY,
c DX,
c DY,
c IP,
c JP,
c NP,
c HELM,
c
c then the preprocessing must be redone (ICFLAG=0 or 1). The
c arrays TRIGS, C, and IPS should also remain unchanged on
c successive calls with ICFLAG=2.
c
c
c ********* NOTES ON DEFINING THE SOLUTION DOMAIN *********
c
c The routine solves an elliptic equation in an irregular,
c polygonal domain. The domain is polygonal because its boundary
c consists of straight line segments parallel or diagonal to
c the discrete mesh. This domain is embedded in a rectangle
c dimensioned NX by NY. The user is free to position the irregular
c domain as desired within the embedding rectangle and the sides of
c the rectangle may be part of the boundary of the solution domain.
c
c The boundary of the irregular domain is described by a set of
c grid points which are classified as 'regular' boundary points
c and 'irregular' boundary points. A regular boundary point is one
c which lies on one of the sides of the embedding rectangle. An
c irregular boundary point is one which is not on a side of the
c embedding rectangle. The routine requires that the user provide
c the coordinates of the irregular boundary points through the
c one-dimensional arrays IP and JP, both dimensioned NP. The user
c must ensure that the specified coordinates of the irregular boundary
c points do not lie on a side of the embedding rectangle, i.e. the
c conditions 1 < IP(N) < NX and 1 < JP(N) < NY must hold for all N.
c Failure to ensure this will result in an error condition. In
c addition the routine requires that there be at least one irregular
c grid point ( NP greater than or equal to 1 ).
c
c Boundary values at the irregular boundary points are supplied
c through the one-dimensional array VP, while boundary values on
c the edge of the embedding rectangle are supplied through the BA,
c BB, BC, and BD arrays. For grid points on the edge of the embedding
c rectangle which are not regular boundary points (i.e. not part of
c the boundary of the solution domain), the corresponding boundary values
c specified through the BA, BB, BC, and BD arrays are arbitrary. For
c definiteness, they should be set to zero.
c
c Three examples follow. In each case the grid points over
c the rectangle are labelled according to the following
c convention: Interior Points - labelled 1,
c Regular Boundary Points - labelled 2,
c Irregular Boundary Points - labelled 3,
c Exterior Points - labelled 0.
c
c **** Example 1 - rectangle with a slot removed ****
c
c NX=17, NY=9, NP=10
c IP=(8,8,8,8,9,10,11,11,11,11)
c JP=(2,3,4,5,5, 5, 5, 4, 3, 2)
c
c 22222222222222222
c 21111111111111112
c 21111111111111112
c 21111111111111112
c 21111113333111112
c 21111113003111112
c 21111113003111112
c 21111113003111112
c 22222222002222222
c
c Boundary values for the problem are specified at the irregular
c boundary points through the VP array, and at the regular boundary
c points through the BA, BB, BC and BD arrays. The
c BC(I), I=1,NX array specifies boundary values on the lower edge
c (J=1) on the embedding rectangle. These values are arbitrary for
c I=9,10 and should be set to zero.
c
c **** Example 2 - rectangle with a triangular hole ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 22222222222222222
c 21111111111111112
c 21133333333333112
c 21113000000031112
c 21111300000311112
c 21111130003111112
c 21111113031111112
c 21111111311111112
c 22222222222222222
c
c **** Example 3 - triangular region ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 00000000000000000
c 00000000000000000
c 00033333333333000
c 00003111111130000
c 00000311111300000
c 00000031113000000
c 00000003130000000
c 00000000300000000
c 00000000000000000
c
c There are no regular boundary points in this case, and so the
c values of the BA, BB, BC and BD arrays are arbitrary. They
c should be set to zero for definiteness.
c
c Note that, in example 3, repositioning the triangular region
c so that one of its sides coincides with one of the sides of the
c embedding rectangle would reduce the number of irregular grid
c points and reduce time for preprocessing and solving.
c
c *********************************************************************
c
REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
REAL C(NP,NP),CHARGE(NP),VP(NP),HELM,HELM2,RHO
REAL BA(NY),BB(NY),BC(NX),BD(NX),DX,DY
INTEGER IPS(NP),IP(NP),JP(NP)
INTEGER IFAX(10)
SAVE LFIRST,NX1,NY1,IFAX
DATA LFIRST/1/
c initialize arrays on first call or for preprocessing
IF ((LFIRST.EQ.1).OR.(ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
IERROR=0
CALL CHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
NX1=NX-1
NY1=NY-1
CALL FAX(IFAX,NX1,4)
CALL FFTRIG(TRIGS,NX1,4)
LFIRST=0
END IF
c If required, execute the error checking routine
IF(IERROR.EQ.-1) THEN
IERROR=0
CALL CHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
END IF
HELM2=HELM*(DY**2)
RHO=DY/DX
c preprocessing: generate the Green's function matrix
IF(ICFLAG.EQ.2) GOTO 5000
IF((ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
ICFG=ICFLAG
CALL CPGN(C,U,WORK,TRIGS,IFAX,IPS,HELM2,
# RHO,NX,NY,IP,JP,NP,IERROR)
ICFLAG=2
ELSE
IERROR=5
RETURN
END IF
IF(ICFG.EQ.0) RETURN
5000 CONTINUE
DO 15 J=1,NY
DO 10 I=1,NX
U(I,J)=F(I,J)*(DY*DY)
10 CONTINUE
15 CONTINUE
c modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the regular boundary
DO 20 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
20 CONTINUE
DO 30 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
30 CONTINUE
c get the solution in the regular domain
CALL POISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)
c obtain the correction to the rhs at the irregular boundary points
CVD$ NOSYNC
CDIR$ IVDEP
DO 40 N1=1,NP
CHARGE(N1)=VP(N1)-U(IP(N1),JP(N1))
40 CONTINUE
CALL SGESL(C,NP,NP,IPS,CHARGE,0)
c reconstruct the rhs
DO 55 J=1,NY
DO 50 I=1,NX
U(I,J)=F(I,J)*(DY*DY)
50 CONTINUE
55 CONTINUE
DO 60 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
60 CONTINUE
DO 70 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
70 CONTINUE
c modify the irregular boundary points
CVD$ NOSYNC
CDIR$ IVDEP
DO 80 N1=1,NP
U(IP(N1),JP(N1))=U(IP(N1),JP(N1))+CHARGE(N1)
80 CONTINUE
c get the solution in the irregular domain
CALL POISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)
c set the boundary conditions on the regular boundary
DO 110 I=1,NX
U(I,1)=BC(I)
U(I,NY)=BD(I)
110 CONTINUE
DO 120 J=2,NY1
U(1,J)=BA(J)
U(NX,J)=BB(J)
120 CONTINUE
CVD$ NOSYNC
CDIR$ IVDEP
DO 130 N1=1,NP
U(IP(N1),JP(N1))=VP(N1)
130 CONTINUE
RETURN
END
c Routine to generate and factor the inverse capacitance matrix
SUBROUTINE CPGN(C,R,WORK,TRIGS,IFAX,IPS,HELM,
# RHO,NX,NY,IP,JP,NP,IERROR)
REAL C(NP,NP),R(NX,NY),WORK(NX,NY),TRIGS(*)
REAL HELM,RHO
INTEGER IP(NP),JP(NP),IPS(NP),IFAX(10)
CVD$ NOSELECT
DO 10 N1=1,NP
DO 25 J=1,NY
DO 20 I=1,NX
R(I,J)=0.
20 CONTINUE
25 CONTINUE
I1=IP(N1)
J1=JP(N1)
R(I1,J1)=1.
CALL POISS(R,HELM,RHO,WORK,TRIGS,IFAX,NX,NY)
DO 30 N2=1,NP
C(N2,N1)=R(IP(N2),JP(N2))
30 CONTINUE
10 CONTINUE
c Linpack routine to factor the matrix
CALL SGEFA(C,NP,NP,IPS,IER)
IF(IER.NE.0) THEN
IERROR=4
END IF
RETURN
END
c
c Routine to check input variables
c
SUBROUTINE CHKV (NX,NY,IP,JP,NP,DX,DY,IERROR)
INTEGER IP(NP), JP(NP)
DATA MAXCK/40/
NX1=NX-1
NY1=NY-1
c Check if the values of NX and NY are too small
IF (NX1.LT.8.OR.NY1.LT.8) THEN
IERROR=1
RETURN
END IF
c Check if the value of NX is admissible
LR=0
LP=0
LQ=0
N1=NX-1
DO 100 IT=1,MAXCK
MX=MOD(N1,2)
IF(MX.NE.0) GOTO 110
LP=LP+1
N1=N1/2
IF(N1.EQ.1) GOTO 401
100 CONTINUE
110 CONTINUE
IF(LP.EQ.0) THEN
IERROR=2
RETURN
END IF
DO 200 IT=1,MAXCK
MX=MOD(N1,3)
IF(MX.NE.0) GOTO 220
LQ=LQ+1
N1=N1/3
IF(N1.EQ.1) GOTO 401
200 CONTINUE
220 CONTINUE
DO 300 IT=1,MAXCK
MX=MOD(N1,5)
IF(MX.NE.0) GOTO 401
LR=LR+1
N1=N1/5
IF(N1.EQ.1) GOTO 401
300 CONTINUE
401 CONTINUE
NXX=(2**LP)*(3**LQ)*(5**LR) + 1
IF(NXX.NE.NX) THEN
IERROR=2
RETURN
END IF
c Check values of the irregular boundary points
DO 500 IN=1,NP
IF(IP(IN).LT.2.OR.IP(IN).GT.NX-1) IERROR=3
IF(JP(IN).LT.2.OR.JP(IN).GT.NY-1) IERROR=3
500 CONTINUE
c Check values of DX and DY
IF((DX.LE.0.).OR.(DY.LE.0.)) IERROR=6
RETURN
END
c***DECK : dcapc.f
SUBROUTINE DCAPC(F,U,WORK,TRIGS,NX,NY,HELM,DX,DY,
# C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
# BA,BB,BC,BD,IERROR)
c
c Programmer: P. Cummins
c Version 1.0
c
c *********************************************************************
c
c DCAPC is a double precision routine to solve the second-order
c accurate finite difference approximation to the Helmholtz equation
c
c (d/dX)(dU/dX) + (d/dY)(dU/dY) - H*U = F(X,Y)
c
c with Dirichlet boundary conditions in irregular polygonal
c two-dimensional domains in Cartesian coordinates using the
c capacitance matrix method.
c
c ******************* ON ENTRY ********************
c
c NX integer constant
c The number of mesh points in the X direction
c Note that the value of NX is restricted such
c that (NX-1) is of the form (2**I * 3**J * 5**K),
c where I is an integer greater than or equal to
c one, and J and K are integers greater than or
c equal to zero. In addition, NX-1 must be greater
c than or equal to 8.
c
c NY integer constant
c The number of mesh points in the Y direction. NY-1
c must be greater than or equal to 8.
c
c NP integer constant
c The number of irregular boundary points. NP must be
c greater than or equal to 1.
c
c F double precision array (NX,NY)
c The source function on the right hand side. The value
c of F is arbitrary at points in the rectangle which are
c outside the irregular solution domain. This array is
c unchanged on exit.
c
c HELM double precision constant
c The coefficient H in the Helmholtz equation given above;
c unchanged on exit.
c
c WORK double precision array (NX,NY)
c A work space array
c
c CHARGE double precision array (NP)
c A work space array
c
c DX double precision constant
c The discretization interval in the X direction.
c The condition DX > 0 must hold.
c
c DY double precision constant
c The discretization interval in the Y direction.
c The condition DY > 0 must hold.
c
c BA double precision array (NY)
c The boundary values on the regular boundary
c at (I=1, J=2,NY-1).
c
c BB double precision array (NY)
c The boundary values on the regular boundary
c at (I=NX, J=2,NY-1).
c
c BC double precision array (nx)
c The boundary values on the regular boundary
c at (I=1,NX, J=1).
c
c BD double precision array (NX)
c The boundary values on the regular boundary
c at (I=1,NX, J=NY).
c
c VP double precision array (NP)
c The boundary values on the irregular boundary,
c defined by the (IP,JP) coordinates.
c
c IP integer array (NP)
c The coordinates in the X direction of the
c irregular boundary points. The condition,
c 1 < IP(N) < NX, must hold for all N.
c
c JP integer array (NP)
c The coordinates in the Y direction of the
c irregular boundary points. The condition,
c 1 < JP(N) < NY, must hold for all N.
c
c Additional notes on defining the IP and JP arrays are
c given below.
c
c Note: The order in which the irregular boundary grid points, defined
c by the (IP(N),JP(N)) pairs, are specified is immaterial.
c However, these grid points must all be distinct; there must
c not be a repeated (IP(N),JP(N)) pair. Failure to ensure this
c will result in the error condition IERROR=4, defined below.
c
c ICFLAG integer constant
c A flag: ICFLAG=0 - do preprocessing only and compute
c the matrix C and indices IPS.
c
c ICFLAG=1 - do preprocessing to compute the
c matrix C and indices IPS, and solve
c the Helmholtz problem.
c
c ICFLAG=2 - solve the Helmholtz problem only; this
c assumes that preprocessing has been
c done and that C and IPS are defined.
c
c Any other value for ICFLAG gives IERROR=5 and no attempt
c is made to do preprocessing or to solve. On return,
c ICFLAG=2.
c
c
c ******************* ON RETURN ********************
c
c U double precision array (NX,NY)
c The solution in the irregular domain, with
c U(IP(N),JP(N)) = VP(N) for N=1,NP on the irregular
c portion of the boundary and also with
c U(1,J)=BA(J) for J=2,NY-1
c U(NX,J)=BB(J) for J=2,NY-1
c U(I,1)=BC(I) for I=1,NX
c U(I,NY)=BD(I) for I=1,NX.
c
c C double precision array (NP,NP)
c The inverse capacitance matrix. This array must not
c be altered on successive calls to the routine with
c different source functions when ICFLAG=2.
c
c IPS integer array (NP)
c An integer vector of pivot indices. This array must
c not be alter on successive calls with different
c source functions when ICFLAG=2.
c
c TRIGS double precision array (2*(NX-1))
c An array of coefficients required by the FFT
c routine employed by the fast solver. This array
c is initialized on the first call to DCAPC, or
c whenever ICFLAG=0 or 1. The array must not be
c altered between successive calls to the routine
c with different source functions when ICFLAG=2.
c
c IERROR integer constant
c An error flag on input variables. On return IERROR
c has the following meaning:
c
c IERROR=0 - normal status, no error detected.
c =1 - value of NX-1 or NY-1 is less than 8.
c =2 - illegal value for NX.
c =3 - illegal value for IP or JP.
c =4 - an error occurred in factoring the inverse
c capacitance matrix, i.e. the matrix is singular.
c This condition will occur if the user has
c defined duplicate irregular boundary points.
c =5 - illegal value for ICFLAG.
c =6 - illegal value for DX or DY.
c
c Note: Error checking is normally only done on the first call
c to the routine. The user may force the routine to do error
c checking at any time by setting IERROR=-1 on entry.
c
c
c ********* NOTES ON CALLING ROUTINE DCAPC ********
c
c In applying the capacitance matrix method to the solution
c of the Helmholtz equation, routine DCAPC works in two
c stages. The first stage is a preprocessing one in which
c the arrays C and IPS are determined. The second stage is
c the actual solution of the problem for a given source
c function. The flag ICFLAG determines the operation of the
c routine:
c
c If the routine is called with ICFLAG=0, only
c preprocessing is done, the arrays C and IPC
c are determined, and DCAPC immediately returns to
c the calling program without solving. On return
c ICFLAG=2
c
c If the routine is called with ICFLAG=1,
c preprocessing is done, the arrays C and IPC
c are determined, and the solution for the
c given source function is determined. On return
c ICFLAG=2.
c
c If the routine is called with ICFLAG=2,
c the solution for the given source function
c is determined. This assumes that preprocessing
c has been done and that the C and IPS arrays
c are defined. On return ICFLAG=2.
c
c If the user wishes to solve the Helmholtz equation in a
c given domain repeatedly with different source functions,
c then the routine could be called with ICFLAG=1 initially.
c Afterwards this flag is set to 2 and need not be altered.
c If the user has stored the arrays C and IPS from a previous
c run, then these may be supplied as input data and the
c routine may be called initially with ICFLAG=2 to obtain
c the solution directly.
c
c Preprocessing need only be done once for a given domain
c geometry and Helmholtz coefficient. The user-specified
c variables listed below must not change between successive
c calls to DCAPC with ICFLAG=2. If the user changes any
c one of the following variables:
c
c NX,
c NY,
c DX,
c DY,
c IP,
c JP,
c NP,
c HELM,
c
c then the preprocessing must be redone (ICFLAG=0 or 1). The
c arrays TRIGS, C, and IPS should also remain unchanged on
c successive calls with ICFLAG=2.
c
c
c ********* NOTES ON DEFINING THE SOLUTION DOMAIN *********
c
c The routine solves an elliptic equation in an irregular,
c polygonal domain. The domain is polygonal because its boundary
c consists of straight line segments parallel or diagonal to
c the discrete mesh. This domain is embedded in a rectangle
c dimensioned NX by NY. The user is free to position the irregular
c domain as desired within the embedding rectangle and the sides of
c the rectangle may be part of the boundary of the solution domain.
c
c The boundary of the irregular domain is described by a set of
c grid points which are classified as 'regular' boundary points
c and 'irregular' boundary points. A regular boundary point is one
c which lies on one of the sides of the embedding rectangle. An
c irregular boundary point is one which is not on a side of the
c embedding rectangle. The routine requires that the user provide
c the coordinates of the irregular boundary points through the
c one-dimensional arrays IP and JP, both dimensioned NP. The user
c must ensure that the specified coordinates of the irregular boundary
c points do not lie on a side of the embedding rectangle, i.e. the
c conditions 1 < IP(N) < NX and 1 < JP(N) < NY must hold for all N.
c Failure to ensure this will result in an error condition. In
c addition the routine requires that there be at least one irregular
c grid point ( NP greater than or equal to 1 ).
c
c Boundary values at the irregular boundary points are supplied
c through the one-dimensional array VP, while boundary values on
c the edge of the embedding rectangle are supplied through the BA,
c BB, BC, and BD arrays. For grid points on the edge of the embedding
c rectangle which are not regular boundary points (i.e. not part of
c the boundary of the solution domain), the corresponding boundary values
c specified through the BA, BB, BC, and BD arrays are arbitrary. For
c definiteness, they should be set to zero.
c
c Three examples follow. In each case the grid points over
c the rectangle are labelled according to the following
c convention: Interior Points - labelled 1,
c Regular Boundary Points - labelled 2,
c Irregular Boundary Points - labelled 3,
c Exterior Points - labelled 0.
c
c **** Example 1 - rectangle with a slot removed ****
c
c NX=17, NY=9, NP=10
c IP=(8,8,8,8,9,10,11,11,11,11)
c JP=(2,3,4,5,5, 5, 5, 4, 3, 2)
c
c 22222222222222222
c 21111111111111112
c 21111111111111112
c 21111111111111112
c 21111113333111112
c 21111113003111112
c 21111113003111112
c 21111113003111112
c 22222222002222222
c
c Boundary values for the problem are specified at the irregular
c boundary points through the VP array, and at the regular boundary
c points through the BA, BB, BC and BD arrays. The
c BC(I), I=1,NX array specifies boundary values on the lower edge
c (J=1) on the embedding rectangle. These values are arbitrary for
c I=9,10 and should be set to zero.
c
c **** Example 2 - rectangle with a triangular hole ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 22222222222222222
c 21111111111111112
c 21133333333333112
c 21113000000031112
c 21111300000311112
c 21111130003111112
c 21111113031111112
c 21111111311111112
c 22222222222222222
c
c **** Example 3 - triangular region ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 00000000000000000
c 00000000000000000
c 00033333333333000
c 00003111111130000
c 00000311111300000
c 00000031113000000
c 00000003130000000
c 00000000300000000
c 00000000000000000
c
c There are no regular boundary points in this case, and so the
c values of the BA, BB, BC and BD arrays are arbitrary. They
c should be set to zero for definiteness.
c
c Note that, in example 3, repositioning the triangular region
c so that one of its sides coincides with one of the sides of the
c embedding rectangle would reduce the number of irregular grid
c points and reduce time for preprocessing and solving.
c
c *********************************************************************
c
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
DOUBLE PRECISION C(NP,NP),CHARGE(NP),VP(NP),HELM,HELM2,RHO
DOUBLE PRECISION BA(NY),BB(NY),BC(NX),BD(NX),DX,DY
INTEGER IPS(NP),IP(NP),JP(NP)
INTEGER IFAX(10)
SAVE LFIRST,NX1,NY1,IFAX
DATA LFIRST/1/
c initialize arrays on first call or for preprocessing
IF ((LFIRST.EQ.1).OR.(ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
IERROR=0
CALL DCHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
NX1=NX-1
NY1=NY-1
CALL DFAX(IFAX,NX1,4)
CALL DFFTRG(TRIGS,NX1,4)
LFIRST=0
END IF
c If required, execute the error checking routine
IF(IERROR.EQ.-1) THEN
IERROR=0
CALL DCHKV(NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
END IF
HELM2=HELM*(DY**2)
RHO=DY/DX
c preprocessing: generate the Green's function matrix
IF(ICFLAG.EQ.2) GOTO 5000
IF((ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
ICFG=ICFLAG
CALL DCPGN(C,U,WORK,TRIGS,IFAX,IPS,HELM2,
# RHO,NX,NY,IP,JP,NP,IERROR)
ICFLAG=2
ELSE
IERROR=5
RETURN
END IF
IF(ICFG.EQ.0) RETURN
5000 CONTINUE
DO 15 J=1,NY
DO 10 I=1,NX
U(I,J)=F(I,J)*(DY*DY)
10 CONTINUE
15 CONTINUE
c modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the regular boundary
DO 20 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
20 CONTINUE
DO 30 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
30 CONTINUE
c get the solution in the regular domain
CALL DPOISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)
c obtain the correction to the rhs at the irregular boundary points
CVD$ NOSYNC
CDIR$ IVDEP
DO 40 N1=1,NP
CHARGE(N1)=VP(N1)-U(IP(N1),JP(N1))
40 CONTINUE
CALL DGESL(C,NP,NP,IPS,CHARGE,0)
c reconstruct the rhs
DO 55 J=1,NY
DO 50 I=1,NX
U(I,J)=F(I,J)*(DY*DY)
50 CONTINUE
55 CONTINUE
DO 60 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
60 CONTINUE
DO 70 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
70 CONTINUE
c modify the irregular boundary points
CVD$ NOSYNC
CDIR$ IVDEP
DO 80 N1=1,NP
U(IP(N1),JP(N1))=U(IP(N1),JP(N1))+CHARGE(N1)
80 CONTINUE
c get the solution in the irregular domain
CALL DPOISS(U,HELM2,RHO,WORK,TRIGS,IFAX,NX,NY)
c set the boundary conditions on the regular boundary
DO 110 I=1,NX
U(I,1)=BC(I)
U(I,NY)=BD(I)
110 CONTINUE
DO 120 J=2,NY1
U(1,J)=BA(J)
U(NX,J)=BB(J)
120 CONTINUE
CVD$ NOSYNC
CDIR$ IVDEP
DO 130 N1=1,NP
U(IP(N1),JP(N1))=VP(N1)
130 CONTINUE
RETURN
END
c Routine to generate and factor the inverse capacitance matrix
SUBROUTINE DCPGN(C,R,WORK,TRIGS,IFAX,IPS,HELM,
# RHO,NX,NY,IP,JP,NP,IERROR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION C(NP,NP),R(NX,NY),WORK(NX,NY),TRIGS(*)
DOUBLE PRECISION HELM,RHO
INTEGER IP(NP),JP(NP),IPS(NP),IFAX(10)
CVD$ NOSELECT
DO 10 N1=1,NP
DO 25 J=1,NY
DO 20 I=1,NX
R(I,J)=0.D0
20 CONTINUE
25 CONTINUE
I1=IP(N1)
J1=JP(N1)
R(I1,J1)=1.D0
CALL DPOISS(R,HELM,RHO,WORK,TRIGS,IFAX,NX,NY)
DO 30 N2=1,NP
C(N2,N1)=R(IP(N2),JP(N2))
30 CONTINUE
10 CONTINUE
c Linpack routine to factor the matrix
CALL DGEFA(C,NP,NP,IPS,IER)
IF(IER.NE.0) THEN
IERROR=4
END IF
RETURN
END
c
c Routine to check input variables
c
SUBROUTINE DCHKV (NX,NY,IP,JP,NP,DX,DY,IERROR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
INTEGER IP(NP), JP(NP)
DATA MAXCK/40/
NX1=NX-1
NY1=NY-1
c Check if the values of NX and NY are too small
IF (NX1.LT.8.OR.NY1.LT.8) THEN
IERROR=1
RETURN
END IF
c Check if the value of NX is admissible
LR=0
LP=0
LQ=0
N1=NX-1
DO 100 IT=1,MAXCK
MX=MOD(N1,2)
IF(MX.NE.0) GOTO 110
LP=LP+1
N1=N1/2
IF(N1.EQ.1) GOTO 401
100 CONTINUE
110 CONTINUE
IF(LP.EQ.0) THEN
IERROR=2
RETURN
END IF
DO 200 IT=1,MAXCK
MX=MOD(N1,3)
IF(MX.NE.0) GOTO 220
LQ=LQ+1
N1=N1/3
IF(N1.EQ.1) GOTO 401
200 CONTINUE
220 CONTINUE
DO 300 IT=1,MAXCK
MX=MOD(N1,5)
IF(MX.NE.0) GOTO 401
LR=LR+1
N1=N1/5
IF(N1.EQ.1) GOTO 401
300 CONTINUE
401 CONTINUE
NXX=(2**LP)*(3**LQ)*(5**LR) + 1
IF(NXX.NE.NX) THEN
IERROR=2
RETURN
END IF
c Check values of the irregular boundary points
DO 500 IN=1,NP
IF(IP(IN).LT.2.OR.IP(IN).GT.NX-1) IERROR=3
IF(JP(IN).LT.2.OR.JP(IN).GT.NY-1) IERROR=3
500 CONTINUE
c Check values of DX and DY
IF((DX.LE.0.D0).OR.(DY.LE.0.D0)) IERROR=6
RETURN
END
c***DECK : reccn.f
SUBROUTINE RECCN(F,G,U,P,WORK,TRIGS,NX,NY,DX,DY,
# BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
c Programmer: P. Cummins
c Version 1.0
c
c *********************************************************************
c
c RECCN is a routine to solve a second-order accurate finite
c difference approximation to the self-adjoint elliptic equation
c
c (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = F(X,Y)
c
c with Dirichlet boundary conditions in rectangular
c two-dimensional domains in Cartesian coordinates.
c
c Note that the diffusion function, G, is staggered with respect
c to the solution and the source function on the grid.
c
c | | |
c --U(I-1,J+1)-------U(I,J+1)----U(I+1,J+1)--
c | | |
c | | |
c | G(I-1,J) | G(I,J) |
c | | |
c | | |
c --U(I-1,J)---------U(I,J)------U(I+1,J)--
c | | |
c | | |
c | G(I-1,J-1) | G(I,J-1) |
c | | |
c | | |
c --U(I-1,J-1)-------U(I,J-1)----U(I-1,J+1)--
c | | |
c
c
c ******************* ON ENTRY ********************
c
c NX integer constant
c The number of mesh points in the X direction
c Note that the value of NX is restricted such
c that (NX-1) is of the form (2**I * 3**J * 5**K),
c where I is an integer greater than or equal to
c one, and J and K are integers greater than or
c equal to zero. In addition, NX-1 must be greater
c than or equal to 8.
c
c NY integer constant
c The number of mesh points in the Y direction. NY-1
c must be greater than or equal to 8.
c
c F real array (NX,NY)
c The source function on the right hand side.
c This array is unchanged on exit.
c
c P real array (NX,NY)
c The initial estimate to the solution; if the user does
c not have an initial estimate, this array should be set
c to zero. This array is overwritten on exit.
c
c G real array (NX-1,NY-1)
c The array containing values of the diffusion function.
c Note that this array is staggered with respect to the
c mesh on which the source function and the solution are
c defined. To ensure a well posed problem, the condition
c G(I,J) > 0 must hold for all points within the rectangle.
c Failure to ensure this will cause condition IERROR=3
c (see below).
c
c WORK real array (NX,NY)
c A work space array.
c
c DX real constant
c The discretization interval in the X direction.
c The condition DX > 0 must hold.
c
c DY real constant
c The discretization interval in the Y direction.
c The condition DY > 0 must hold.
c
c BA real array (NY)
c The boundary values at (I=1, J=2,NY-1).
c
c BB real array (NY)
c The boundary values at (I=NX, J=2,NY-1).
c
c BC real array (nx)
c The boundary values at (I=1,NX, J=1).
c
c BD real array (NX)
c The boundary values at (I=1,NX, J=NY).
c
c CCON real constant
c The convergence criterion. The routine will proceed
c with the iteration to reach a solution until the two
c following normalized conditions are satisfied
c
c (1) || U(k)-U(k-1) || / ||U(k)|| <= CCON
c
c (2) || L(U(k))-F || / ||F|| <= CCON,
c
c or until the number of iterations has reached the maximum
c allowable, given by parameter ITMAX. Note that ||U(k)||
c refers to the L2 norm of the kth iterate of U and that L
c refers to the discrete elliptic operator. If ||F||=0, then
c convergence is determined from condition (1) alone. CCON
c must be greater than zero or the routine returns without
c attempting a solution and IERROR=6.
c
c The minimum possible value of CCON is limited by
c the accumulation of round-off error and depends
c mainly on the machine precision and the size of
c the grid. The following table of approximate minimum
c values of CCON was constructed from solution of the
c elliptic equation, in a unit square domain with the
c source and diffusion functions given by a uniformly
c distributed random number, on an Alliant FX/40. This table
c is intended as a guide to the user; the actual minimum
c attainable value of CCON depends on the particular problem
c and on the computer used.
c
c (NX,NY) Single Precision (32 bit) Double Precision (64 bit)
c
c (17,17) 4.E-6 2.E-14
c
c (33,33) 3.E-5 5.E-14
c
c (65,65) 2.E-4 4.E-13
c
c (129,129) 2.E-3 4.E-12
c
c (257,257) 2.E-2 4.E-11
c
c (513,513) ----- 2.E-10
c
c (1025,1025) ----- 2.E-09
c
c
c Under normal operation, the solver will iterate until
c the convergence conditions are satisfied and then
c return. The convergence of the routine is monitored
c at each iteration and if a lack of convergence is
c detected, then the routine will return with
c IERROR=5. Specifically the ratio D(k+1)/D(k), where
c D(k)=|| U(k)-U(k-1) ||, is monitored for cases when
c it is greater than one. Should this occur repeatedly
c (three times over six iterations), then lack of
c convergence is assumed and the routine returns to
c the calling program with IERROR=5. This feature is
c invoked automatically, but may be disabled by
c setting ITMAX (see below) to a negative value.
c
c One situation which may lead to a lack of convergence
c arises if CCON is set to a value which is too small
c for a given problem. In such cases, the routine
c initially converges to the solution for a number of
c iterations and then fails to converge further with
c additional iterations. This may be monitored in the
c output to unit IWRITE (see below). The remedy is to
c increase the value of CCON or to use the double precision
c routine DRECCN.
c
c IWRITE integer constant
c Logical unit number to which the routine will write out
c the iteration number, k, and the convergence conditions,
c
c ANORM = || U(k)-U(k-1)|| / ||U(k)||
c
c RESIDUAL = || L(U(k))-F || / ||F||.
c
c This is useful for examining the iteration by iteration
c convergence of the solver. IWRITE=0 suppresses this output.
c If ||F||=0, then RESIDUAL is set to zero and convergence
c depends on ANORM.
c
c ITMAX integer constant
c The absolute value, |ITMAX|, gives the maximum number of
c iteration permitted. If the number of iterations reaches
c this maximum without obtaining a solution to the required
c accuracy, then IERROR=4 on output. The routine accepts
c positive or negative values for ITMAX, with negative
c values suppressing the iterative check on convergence
c described above under parameter CCON.
c
c ******************* ON RETURN ********************
c
c U real array (NX,NY)
c The solution in the rectangular domain, with
c boundary conditions
c U(1,J)=BA(J) for J=2,NY-1
c U(NX,J)=BB(J) for J=2,NY-1
c U(I,1)=BC(I) for I=1,NX
c U(I,NY)=BD(I) for I=1,NX.
c
c NITT integer constant
c The number of iterations required to get the
c solution within the prescribed margin of error
c given by CCON.
c
c TRIGS real array (2*(NX-1))
c An array of coefficients required by the FFT
c routine employed by the fast solver. This array
c is initialized on the first call to RECCN and
c must not be altered between successive calls to
c the routine with different source or diffusion
c functions.
c
c IERROR integer constant
c An error flag on input variables. On return IERROR
c has the following meaning:
c
c IERROR=0 - normal status, no error detected.
c =1 - value of NX-1 or NY-1 is less than 8.
c =2 - illegal value for NX.
c =3 - an illegal value of G detected.
c =4 - maximum number of iterations reached.
c =5 - nonconvergence of the iteration detected.
c =6 - CCON is less than or equal to zero.
c =7 - illegal value for DX or DY.
c
c Note: Error checking is normally only done on the first call
c to the routine. The user may force the routine to do
c error checking on subsequent calls by setting
c IERROR=-1 on entry.
c
c *********************************************************************
c
PARAMETER(MMAX=6,MMAX1=MMAX-1)
REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
REAL P(NX,NY),G(NX-1,NY-1),RHO
REAL BA(NY),BB(NY),BC(NX),BD(NX),AN1(MMAX)
INTEGER IFAX(10)
SAVE LFIRST,NX1,NY1,NXSAVE,NYSAVE,IFAX
DATA LFIRST/1/
c initialize arrays
IF ((LFIRST.EQ.1).OR.(NXSAVE.NE.NX).OR.(NYSAVE.NE.NY)) THEN
IERROR=0
CALL CHKV3(G,NX,NY,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
NX1=NX-1
NY1=NY-1
CALL FAX(IFAX,NX1,4)
CALL FFTRIG(TRIGS,NX1,4)
NXSAVE=NX
NYSAVE=NY
LFIRST=0
END IF
c If required, execute the error checking routine
IF(IERROR.EQ.-1) THEN
IERROR=0
CALL CHKV3(G,NX,NY,DX,DY,IERROR)
IF (IERROR.NE.0) RETURN
END IF
IF(CCON.LE.0.) THEN
IERROR=6
RETURN
END IF
NITT=1
DO 5 M=1,MMAX
AN1(M)=0.
5 CONTINUE
RHO=DY/DX
DO 10 I=1,NX
P(I,1)=BC(I)
P(I,NY)=BD(I)
10 CONTINUE
DO 20 J=2,NY1
P(1,J)=BA(J)
P(NX,J)=BB(J)
20 CONTINUE
1000 CONTINUE
c form the right hand side
DO 35 J=2,NY1
DO 30 I=2,NX1
U(I,J)=(4.*DY*DY*F(I,J)-RHO*RHO*(P(I+1,J)-P(I-1,J))*
# (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
# - (P(I,J+1)-P(I,J-1))*
# (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
# (G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
30 CONTINUE
35 CONTINUE
c modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the boundary
DO 40 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
40 CONTINUE
DO 50 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
50 CONTINUE
c get the solution in the domain
CALL POISS(U,0.,RHO,WORK,TRIGS,IFAX,NX,NY)
c set the boundary conditions
DO 120 I=1,NX
U(I,1)=BC(I)
U(I,NY)=BD(I)
120 CONTINUE
DO 130 J=2,NY1
U(1,J)=BA(J)
U(NX,J)=BB(J)
130 CONTINUE
c check for convergence of the solution
ADEN=0.
ADIFF=0.
DL2N=0.
DENM=0.
CVD$ ASSOC (D)
DO 145 J=2,NY1
DO 140 I=2,NX1
ADEN=ADEN+U(I,J)*U(I,J)
ADIFF=ADIFF+(U(I,J)-P(I,J))*(U(I,J)-P(I,J))
FD=0.5*((G(I,J)+G(I,J-1))*(U(I+1,J)-U(I,J))
# - (G(I-1,J)+G(I-1,J-1))*(U(I,J)-U(I-1,J)))/(DX*DX)
# +0.5*((G(I,J)+G(I-1,J))*(U(I,J+1)-U(I,J))
# - (G(I,J-1)+G(I-1,J-1))*(U(I,J)-U(I,J-1)))/(DY*DY)
DL2N=DL2N+(FD-F(I,J))**2
DENM=DENM+F(I,J)**2
140 CONTINUE
145 CONTINUE
IF(ADEN.NE.0.) THEN
ANORM=SQRT(ADIFF/ADEN)
ELSE
ANORM=SQRT(ADIFF)
END IF
IF(DENM.NE.0.) THEN
RRESD=SQRT(DL2N/DENM)
ELSE IF(ADEN.NE.0.) THEN
RRESD=SQRT(DL2N/ADEN)
ELSE
RRESD=SQRT(DL2N)
END IF
c If required, print out the iteration by iteration convergence
c of the solver.
IF(IWRITE.NE.0) THEN
WRITE(IWRITE,101) NITT,ANORM,RRESD
END IF
101 FORMAT(' ','ITERATION ',I4,' ANORM= ',1PG11.5,' RESIDUAL= ',
# 1PG11.5)
IF((ANORM.LE.CCON).AND.(RRESD.LE.CCON)) RETURN
IF(ITMAX.GT.0) THEN
DO 160 M=2,MMAX
MM=MMAX-M+2
AN1(MM)=AN1(MM-1)
160 CONTINUE
AN1(1)=SQRT(ADIFF)
RTSUM1=0.
DO 170 M=1,MMAX1
RAT1=0.
IF(AN1(M+1).NE.0.) THEN
RAT1=AN1(M)/AN1(M+1)
IF(RAT1.GT.1.) RTSUM1=RTSUM1+1.
END IF
170 CONTINUE
IF((RTSUM1.GT.2.99)) THEN
IERROR=5
RETURN
END IF
END IF
IF(NITT.GE.ABS(ITMAX)) THEN
IERROR=4
RETURN
END IF
DO 600 J=1,NY
DO 500 I=1,NX
P(I,J)=U(I,J)
500 CONTINUE
600 CONTINUE
NITT=NITT+1
GO TO 1000
END
c
c Routine to check input variables and the G array
c
SUBROUTINE CHKV3 (G,NX,NY,DX,DY,IERROR)
REAL G(NX-1,NY-1),DX,DY
DATA MAXCK/40/
NX1=NX-1
NY1=NY-1
c Check if the values of NX and NY are too small
IF (NX1.LT.8.OR.NY1.LT.8) THEN
IERROR=1
RETURN
END IF
c Check if the value of NX is admissible
LR=0
LP=0
LQ=0
N1=NX-1
DO 100 IT=1,MAXCK
MX=MOD(N1,2)
IF(MX.NE.0) GOTO 110
LP=LP+1
N1=N1/2
IF(N1.EQ.1) GOTO 401
100 CONTINUE
110 CONTINUE
IF(LP.EQ.0) THEN
IERROR=2
RETURN
END IF
DO 200 IT=1,MAXCK
MX=MOD(N1,3)
IF(MX.NE.0) GOTO 220
LQ=LQ+1
N1=N1/3
IF(N1.EQ.1) GOTO 401
200 CONTINUE
220 CONTINUE
DO 300 IT=1,MAXCK
MX=MOD(N1,5)
IF(MX.NE.0) GOTO 401
LR=LR+1
N1=N1/5
IF(N1.EQ.1) GOTO 401
300 CONTINUE
401 CONTINUE
NXX=(2**LP)*(3**LQ)*(5**LR) + 1
IF(NXX.NE.NX) THEN
IERROR=2
RETURN
END IF
IF(IERROR.NE.0) RETURN
c Check if the diffusion function is everywhere positive over
c the interior of the rectangular domain
DO 750 J=1,NY1
DO 700 I=1,NX1
IF(G(I,J).LE.0.) IERROR=3
700 CONTINUE
750 CONTINUE
c Check values of DX and DY
IF((DX.LE.0.).OR.(DY.LE.0.)) IERROR=7
RETURN
END
c***DECK : dreccn.f
SUBROUTINE DRECCN(F,G,U,P,WORK,TRIGS,NX,NY,DX,DY,
# BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
c Programmer: P. Cummins
c Version 1.0
c
c *********************************************************************
c
c DRECCN is a double precision routine to solve a second-order
c accurate finite difference approximation to the self-adjoint
c elliptic equation
c
c (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = F(X,Y)
c
c with Dirichlet boundary conditions in rectangular
c two-dimensional domains in Cartesian coordinates.
c
c Note that the diffusion function, G, is staggered with respect
c to the solution and the source function on the grid.
c
c | | |
c --U(I-1,J+1)-------U(I,J+1)----U(I+1,J+1)--
c | | |
c | | |
c | G(I-1,J) | G(I,J) |
c | | |
c | | |
c --U(I-1,J)---------U(I,J)------U(I+1,J)--
c | | |
c | | |
c | G(I-1,J-1) | G(I,J-1) |
c | | |
c | | |
c --U(I-1,J-1)-------U(I,J-1)----U(I-1,J+1)--
c | | |
c
c
c ******************* ON ENTRY ********************
c
c NX integer constant
c The number of mesh points in the X direction
c Note that the value of NX is restricted such
c that (NX-1) is of the form (2**I * 3**J * 5**K),
c where I is an integer greater than or equal to
c one, and J and K are integers greater than or
c equal to zero. In addition, NX-1 must be greater
c than or equal to 8.
c
c NY integer constant
c The number of mesh points in the Y direction. NY-1
c must be greater than or equal to 8.
c
c F double precision array (NX,NY)
c The source function on the right hand side.
c This array is unchanged on exit.
c
c P double precision array (NX,NY)
c The initial estimate to the solution; if the user does
c not have an initial estimate, this array should be set
c to zero. This array is overwritten on exit.
c
c G double precision array (NX-1,NY-1)
c The array containing values of the diffusion function.
c Note that this array is staggered with respect to the
c mesh on which the source function and the solution are
c defined. To ensure a well posed problem, the condition
c G(I,J) > 0 must hold for all points within the rectangle.
c Failure to ensure this will cause condition IERROR=3
c (see below).
c
c WORK double precision array (NX,NY)
c A work space array.
c
c DX double precision constant
c The discretization interval in the X direction.
c The condition DX > 0 must hold.
c
c DY double precision constant
c The discretization interval in the Y direction.
c The condition DY > 0 must hold.
c
c BA double precision array (NY)
c The boundary values at (I=1, J=2,NY-1).
c
c BB double precision array (NY)
c The boundary values at (I=NX, J=2,NY-1).
c
c BC double precision array (nx)
c The boundary values at (I=1,NX, J=1).
c
c BD double precision array (NX)
c The boundary values at (I=1,NX, J=NY).
c
c CCON double precision constant
c The convergence criterion. The routine will proceed
c with the iteration to reach a solution until the two
c following normalized conditions are satisfied
c
c (1) || U(k)-U(k-1) || / ||U(k)|| <= CCON
c
c (2) || L(U(k))-F || / ||F|| <= CCON,
c
c or until the number of iterations has reached the maximum
c allowable, given by parameter ITMAX. Note that ||U(k)||
c refers to the L2 norm of the kth iterate of U and that L
c refers to the discrete elliptic operator. If ||F||=0, then
c convergence is determined from condition (1) alone. CCON
c must be greater than zero or the routine returns without
c attempting a solution and IERROR=6.
c
c The minimum possible value of CCON is limited by
c the accumulation of round-off error and depends
c mainly on the machine precision and the size of
c the grid. The following table of approximate minimum
c values of CCON was constructed from solution of the
c elliptic equation, in a unit square domain with the
c source and diffusion functions given by a uniformly
c distributed random number, on an Alliant FX/40. This table
c is intended as a guide to the user; the actual minimum
c attainable value of CCON depends on the particular problem
c and on the computer used.
c
c (NX,NY) Single Precision (32 bit) Double Precision (64 bit)
c
c (17,17) 4.E-6 2.E-14
c
c (33,33) 3.E-5 5.E-14
c
c (65,65) 2.E-4 4.E-13
c
c (129,129) 2.E-3 4.E-12
c
c (257,257) 2.E-2 4.E-11
c
c (513,513) ----- 2.E-10
c
c (1025,1025) ----- 2.E-09
c
c
c Under normal operation, the solver will iterate until
c the convergence conditions are satisfied and then
c return. The convergence of the routine is monitored
c at each iteration and if a lack of convergence is
c detected, then the routine will return with
c IERROR=5. Specifically the ratio D(k+1)/D(k), where
c D(k)=|| U(k)-U(k-1) ||, is monitored for cases when
c it is greater than one. Should this occur repeatedly
c (three times over six iterations), then lack of
c convergence is assumed and the routine returns to
c the calling program with IERROR=5. This feature is
c invoked automatically, but may be disabled by
c setting ITMAX (see below) to a negative value.
c
c One situation which may lead to a lack of convergence
c arises if CCON is set to a value which is too small
c for a given problem. In such cases, the routine
c initially converges to the solution for a number of
c iterations and then fails to converge further with
c additional iterations. This may be monitored in the
c output to unit IWRITE (see below). The remedy is to
c increase the value of CCON or to use the double precision
c routine DRECCN.
c
c IWRITE integer constant
c Logical unit number to which the routine will write out
c the iteration number, k, and the convergence conditions,
c
c ANORM = || U(k)-U(k-1)|| / ||U(k)||
c
c RESIDUAL = || L(U(k))-F || / ||F||.
c
c This is useful for examining the iteration by iteration
c convergence of the solver. IWRITE=0 suppresses this output.
c If ||F||=0, then RESIDUAL is set to zero and convergence
c depends on ANORM.
c
c ITMAX integer constant
c The absolute value, |ITMAX|, gives the maximum number of
c iteration permitted. If the number of iterations reaches
c this maximum without obtaining a solution to the required
c accuracy, then IERROR=4 on output. The routine accepts
c positive or negative values for ITMAX, with negative
c values suppressing the iterative check on convergence
c described above under parameter CCON.
c
c ******************* ON RETURN ********************
c
c U double precision array (NX,NY)
c The solution in the rectangular domain, with
c boundary conditions
c U(1,J)=BA(J) for J=2,NY-1
c U(NX,J)=BB(J) for J=2,NY-1
c U(I,1)=BC(I) for I=1,NX
c U(I,NY)=BD(I) for I=1,NX.
c
c NITT integer constant
c The number of iterations required to get the
c solution within the prescribed margin of error
c given by CCON.
c
c TRIGS double precision array (2*(NX-1))
c An array of coefficients required by the FFT
c routine employed by the fast solver. This array
c is initialized on the first call to DRECCN and
c must not be altered between successive calls to
c the routine with different source or diffusion
c functions.
c
c IERROR integer constant
c An error flag on input variables. On return IERROR
c has the following meaning:
c
c IERROR=0 - normal status, no error detected.
c =1 - value of NX-1 or NY-1 is less than 8.
c =2 - illegal value for NX.
c =3 - an illegal value of G detected.
c =4 - maximum number of iterations reached.
c =5 - nonconvergence of the iteration detected.
c =6 - CCON is less than or equal to zero.
c =7 - illegal value for DX or DY.
c
c Note: Error checking is normally only done on the first call
c to the routine. The user may force the routine to do
c error checking on subsequent calls by setting
c IERROR=-1 on entry.
c
c *********************************************************************
c
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER(MMAX=6,MMAX1=MMAX-1)
DOUBLE PRECISION F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
DOUBLE PRECISION P(NX,NY),G(NX-1,NY-1),RHO
DOUBLE PRECISION BA(NY),BB(NY),BC(NX),BD(NX),AN1(MMAX)
INTEGER IFAX(10)
SAVE LFIRST,NX1,NY1,NXSAVE,NYSAVE,IFAX
DATA LFIRST/1/
c initialize arrays
IF ((LFIRST.EQ.1).OR.(NXSAVE.NE.NX).OR.(NYSAVE.NE.NY)) THEN
IERROR=0
CALL DCHKV3(G,NX,NY,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
NX1=NX-1
NY1=NY-1
CALL DFAX(IFAX,NX1,4)
CALL DFFTRG(TRIGS,NX1,4)
NXSAVE=NX
NYSAVE=NY
LFIRST=0
END IF
c If required, execute the error checking routine
IF(IERROR.EQ.-1) THEN
IERROR=0
CALL DCHKV3(G,NX,NY,DX,DY,IERROR)
IF (IERROR.NE.0) RETURN
END IF
IF(CCON.LE.0.) THEN
IERROR=6
RETURN
END IF
NITT=1
DO 5 M=1,MMAX
AN1(M)=0.D0
5 CONTINUE
RHO=DY/DX
DO 10 I=1,NX
P(I,1)=BC(I)
P(I,NY)=BD(I)
10 CONTINUE
DO 20 J=2,NY1
P(1,J)=BA(J)
P(NX,J)=BB(J)
20 CONTINUE
1000 CONTINUE
c form the right hand side
DO 35 J=2,NY1
DO 30 I=2,NX1
U(I,J)=(4.D0*DY*DY*F(I,J)-RHO*RHO*(P(I+1,J)-P(I-1,J))*
# (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
# - (P(I,J+1)-P(I,J-1))*
# (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
# (G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
30 CONTINUE
35 CONTINUE
c modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the boundary
DO 40 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
40 CONTINUE
DO 50 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
50 CONTINUE
c get the solution in the domain
CALL DPOISS(U,0.D0,RHO,WORK,TRIGS,IFAX,NX,NY)
c set the boundary conditions
DO 120 I=1,NX
U(I,1)=BC(I)
U(I,NY)=BD(I)
120 CONTINUE
DO 130 J=2,NY1
U(1,J)=BA(J)
U(NX,J)=BB(J)
130 CONTINUE
c check for convergence of the solution
ADEN=0.D0
ADIFF=0.D0
DL2N=0.D0
DENM=0.D0
CVD$ ASSOC (D)
DO 145 J=2,NY1
DO 140 I=2,NX1
ADEN=ADEN+U(I,J)*U(I,J)
ADIFF=ADIFF+(U(I,J)-P(I,J))*(U(I,J)-P(I,J))
FD=0.5D0*((G(I,J)+G(I,J-1))*(U(I+1,J)-U(I,J))
# - (G(I-1,J)+G(I-1,J-1))*(U(I,J)-U(I-1,J)))/(DX*DX)
# +0.5D0*((G(I,J)+G(I-1,J))*(U(I,J+1)-U(I,J))
# - (G(I,J-1)+G(I-1,J-1))*(U(I,J)-U(I,J-1)))/(DY*DY)
DL2N=DL2N+(FD-F(I,J))**2
DENM=DENM+F(I,J)**2
140 CONTINUE
145 CONTINUE
IF(ADEN.NE.0.D0) THEN
ANORM=SQRT(ADIFF/ADEN)
ELSE
ANORM=SQRT(ADIFF)
END IF
IF(DENM.NE.0.D0) THEN
RRESD=SQRT(DL2N/DENM)
ELSE IF(ADEN.NE.0.D0) THEN
RRESD=SQRT(DL2N/ADEN)
ELSE
RRESD=SQRT(DL2N)
END IF
c If required, print out the iteration by iteration convergence
c of the solver.
IF(IWRITE.NE.0) THEN
WRITE(IWRITE,101) NITT,ANORM,RRESD
END IF
101 FORMAT(' ','ITERATION ',I4,' ANORM= ',1PG11.5,' RESIDUAL= ',
# 1PG11.5)
IF((ANORM.LE.CCON).AND.(RRESD.LE.CCON)) RETURN
IF(ITMAX.GT.0) THEN
DO 160 M=2,MMAX
MM=MMAX-M+2
AN1(MM)=AN1(MM-1)
160 CONTINUE
AN1(1)=SQRT(ADIFF)
RTSUM1=0.D0
DO 170 M=1,MMAX1
RAT1=0.D0
IF(AN1(M+1).NE.0.D0) THEN
RAT1=AN1(M)/AN1(M+1)
IF(RAT1.GT.1.D0) RTSUM1=RTSUM1+1.D0
END IF
170 CONTINUE
IF((RTSUM1.GT.2.99D0)) THEN
IERROR=5
RETURN
END IF
END IF
IF(NITT.GE.ABS(ITMAX)) THEN
IERROR=4
RETURN
END IF
DO 600 J=1,NY
DO 500 I=1,NX
P(I,J)=U(I,J)
500 CONTINUE
600 CONTINUE
NITT=NITT+1
GO TO 1000
END
c
c Routine to check input variables and the G array
c
SUBROUTINE DCHKV3 (G,NX,NY,DX,DY,IERROR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION G(NX-1,NY-1),DX,DY
DATA MAXCK/40/
NX1=NX-1
NY1=NY-1
c Check if the values of NX and NY are too small
IF (NX1.LT.8.OR.NY1.LT.8) THEN
IERROR=1
RETURN
END IF
c Check if the value of NX is admissible
LR=0
LP=0
LQ=0
N1=NX-1
DO 100 IT=1,MAXCK
MX=MOD(N1,2)
IF(MX.NE.0) GOTO 110
LP=LP+1
N1=N1/2
IF(N1.EQ.1) GOTO 401
100 CONTINUE
110 CONTINUE
IF(LP.EQ.0) THEN
IERROR=2
RETURN
END IF
DO 200 IT=1,MAXCK
MX=MOD(N1,3)
IF(MX.NE.0) GOTO 220
LQ=LQ+1
N1=N1/3
IF(N1.EQ.1) GOTO 401
200 CONTINUE
220 CONTINUE
DO 300 IT=1,MAXCK
MX=MOD(N1,5)
IF(MX.NE.0) GOTO 401
LR=LR+1
N1=N1/5
IF(N1.EQ.1) GOTO 401
300 CONTINUE
401 CONTINUE
NXX=(2**LP)*(3**LQ)*(5**LR) + 1
IF(NXX.NE.NX) THEN
IERROR=2
RETURN
END IF
IF(IERROR.NE.0) RETURN
c Check if the diffusion function is everywhere positive over
c the interior of the rectangular domain
DO 750 J=1,NY1
DO 700 I=1,NX1
IF(G(I,J).LE.0.D0) IERROR=3
700 CONTINUE
750 CONTINUE
c Check values of DX and DY
IF((DX.LE.0.D0).OR.(DY.LE.0.D0)) IERROR=7
RETURN
END
c***DECK : capcn.f
SUBROUTINE CAPCN(F,G,U,P,WORK,MASK,TRIGS,NX,NY,DX,DY,
# C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
# BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
c Programmer: P. Cummins
c Version 1.0
c
c *********************************************************************
c
c CAPCN is a routine to solve a second-order accurate finite
c difference approximation to the self-adjoint elliptic
c equation
c
c (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = F(X,Y)
c
c with Dirichlet boundary conditions in irregular polygonal
c two-dimensional domains in Cartesian coordinates using the
c capacitance matrix method.
c
c The diffusion function, G, is staggered with respect
c to the solution and the source function on the grid.
c
c | | |
c --U(I-1,J+1)-------U(I,J+1)----U(I+1,J+1)--
c | | |
c | | |
c | G(I-1,J) | G(I,J) |
c | | |
c | | |
c --U(I-1,J)---------U(I,J)------U(I+1,J)--
c | | |
c | | |
c | G(I-1,J-1) | G(I,J-1) |
c | | |
c | | |
c --U(I-1,J-1)-------U(I,J-1)----U(I-1,J+1)--
c | | |
c
c
c ******************* ON ENTRY ********************
c
c NX integer constant
c The number of mesh points in the X direction
c Note that the value of NX is restricted such
c that (NX-1) is of the form (2**I * 3**J * 5**K),
c where I is an integer greater than or equal to
c one, and J and K are integers greater than or
c equal to zero. In addition, NX-1 must be greater
c than or equal to 8.
c
c NY integer constant
c The number of mesh points in the Y direction. NY-1
c must be greater than or equal to 8.
c
c NP integer constant
c The number of irregular boundary points. NP must be
c greater than or equal to one.
c
c F real array (NX,NY)
c The source function on the right hand side. The value
c of F is arbitrary at points in the rectangle which are
c outside the irregular solution domain. This array is
c unchanged on exit.
c
c P real array (NX,NY)
c The initial estimate to the solution; if the user does
c not have an initial estimate, this array should be set
c to zero. This array is overwritten on exit.
c
c G real array (NX-1,NY-1)
c The array containing values of the diffusion function.
c Note that this array is staggered with respect to the
c mesh on which the source function and the solution are
c defined. To ensure a well posed problem, the condition
c G(I,J) > 0 must hold for all points within the interior
c of the irregular domain. The value of G is arbitrary at
c points in the rectangle which are outside the irregular
c solution domain. This array is unchanged on exit.
c
c MASK integer array (NX,NY)
c This array is used by the routine to identify the mesh points
c within the embedding rectangle that lie inside the irregular
c (polygonal) domain. The user must set MASK(I,J)=1 for all
c (I,J) values which are inside the irregular domain and
c MASK(I,J)=0 for all other points, including the boundary
c points. The mesh of grid points defined by the MASK array
c coincides with that of the source function and the solution.
c This array is unchanged on exit.
c
c Additional notes on defining the MASK array are given below.
c
c WORK real array (NX,NY)
c A work space array.
c
c CHARGE real array (NP)
c A work space array.
c
c DX real constant
c The discretization interval in the X direction.
c The condition DX > 0 must hold.
c
c DY real constant
c The discretization interval in the Y direction.
c The condition DY > 0 must hold.
c
c BA real array (NY)
c The boundary values on the regular boundary
c at (I=1, J=2,NY-1).
c
c BB real array (NY)
c The boundary values on the regular boundary
c at (I=NX, J=2,NY-1).
c
c BC real array (NX)
c The boundary values on the regular boundary
c at (I=1,NX, J=1).
c
c BD real array (NX)
c The boundary values on the regular boundary
c at (I=1,NX, J=NY).
c
c VP real array (NP)
c The boundary values on the irregular boundary.
c
c IP integer array (NP)
c The coordinates in the X direction of the
c irregular boundary points. The condition,
c 1 < IP(N) < NX, must hold for all N.
c
c JP integer array (NP)
c The coordinates in the Y direction of the
c irregular boundary points. The condition,
c 1 < JP(N) < NY, must hold for all N.
c
c Note: The order in which the irregular boundary grid points, defined
c by the (IP(N),JP(N)) pairs, are specified is immaterial.
c However, these grid points must all be distinct; there must
c not be a repeated (IP(N),JP(N)) pair. Failure to ensure this
c will result in the error condition IERROR=7 (see below).
c
c Additional notes on defining the IP and JP arrays are
c given below.
c
c ICFLAG integer constant
c A flag: ICFLAG=0 - do preprocessing only and compute
c the matrix C and indices IPS.
c
c ICFLAG=1 - do preprocessing to compute the
c matrix C and indices IPS, and solve
c the elliptic problem.
c
c ICFLAG=2 - solve the elliptic problem only; this
c assumes that preprocessing has been
c done and that C and IPS are defined.
c
c Any other value for ICFLAG gives IERROR=5 and no attempt
c is made to do preprocessing or to solve. On return,
c ICFLAG=2.
c
c CCON real constant
c The convergence criterion. The routine will proceed with
c the iteration to reach a solution until the following
c two conditions are satisfied
c
c (1) || U(k)-U(k-1) || / ||U(k)|| <= CCON
c
c (2) || L(U(k))-F || / ||F|| <= CCON,
c
c or until the number of iterations has reached the maximum
c allowable, given by parameter ITMAX. Note that ||U(k)||
c refers to the L2 norm of the kth iterate of U and that L
c refers to the discrete elliptic operator. If ||F||=0, then
c convergence is determined from condition (1) alone. CCON
c must be greater than zero or the routine returns without
c attempting a solution and IERROR=10.
c
c The minimum possible value of CCON is limited by
c the accumulation of round-off error and depends
c mainly on the machine precision and the size of
c the grid. The following table of approximate minimum
c values of CCON was constructed from solution of the
c elliptic equation, in a unit square domain with the
c source and diffusion functions given by a uniformly
c distributed random number, on an Alliant FX/40. This table
c is intended as a guide to the user; the actual minimum
c attainable value of CCON depends on the particular problem
c and on the computer used.
c
c (NX,NY) Single Precision (32 bit) Double Precision (64 bit)
c
c (17,17) 4.E-6 2.E-14
c
c (33,33) 3.E-5 5.E-14
c
c (65,65) 2.E-4 4.E-13
c
c (129,129) 2.E-3 4.E-12
c
c (257,257) 2.E-2 4.E-11
c
c (513,513) ----- 2.E-10
c
c (1025,1025) ----- 2.E-09
c
c Under normal operation, the solver will iterate until
c the convergence conditions are satisfied and then
c return. The convergence of the routine is monitored
c at each iteration and if a lack of convergence is
c detected, then the routine will return with
c IERROR=9. Specifically the ratio D(k+1)/D(k), where
c D(k)=|| U(k+1)-U(k) ||, is monitored for cases when
c it is greater than one. Should this occur repeatedly
c (three times over six iterations), then lack of
c convergence is assumed and the routine returns to
c the calling program with IERROR=9. This feature is
c invoked automatically, but may be disabled by
c setting ITMAX (see below) to a negative value.
c
c One situation which may lead to a lack of convergence
c arises if CCON is set to a value which is too small
c for a given machine precision. In such cases, the routine
c initially converges to the solution for a number of
c iterations and then fails to converge further with
c additional iterations. This may be monitored in the
c output to unit IWRITE (see below). The remedy is to
c increase the value of CCON or to use the double precision
c routine DCAPCN.
c
c IWRITE integer constant
c Logical unit number to which the routine will write out
c the iteration number, k, and the convergence conditions,
c
c ANORM = || U(k)-U(k-1)|| / ||U(k)||
c
c RESIDUAL = || L(U(k))-F || / ||F||.
c
c This is useful for examining the iteration by iteration
c convergence of the solver. IWRITE=0 suppresses this output.
c
c ITMAX integer constant
c The absolute value, |ITMAX|, gives the maximum number of
c iteration permitted. If the number of iterations reaches
c this maximum without obtaining a solution to the required
c accuracy, then IERROR=4 on output. The routine accepts
c positive or negative values for ITMAX, with negative
c values suppressing the iterative check on convergence
c described above under parameter CCON.
c
c ******************* ON RETURN ********************
c
c U real array (NX,NY)
c The solution in the irregular domain, with
c U(IP(N),JP(N)) = VP(N) for N=1,NP on the irregular
c portion of the boundary and also with
c U(1,J)=BA(J) for J=2,NY-1
c U(NX,J)=BB(J) for J=2,NY-1
c U(I,1)=BC(I) for I=1,NX
c U(I,NY)=BD(I) for I=1,NX.
c
c NITT integer constant
c The number of iterations required to get the
c solution within the prescribed margin of error
c given by CCON.
c
c C real array (NP,NP)
c The inverse capacitance matrix. This array must not be
c be altered on successive calls to the routine with
c different source functions when ICFLAG=2.
c
c IPS integer array (NP)
c An integer vector of pivot indices. This array must not
c not be alter on successive calls to the routine with
c different source functions when ICFLAG=2.
c
c TRIGS real array (2*(NX-1))
c An array of coefficients required by the FFT
c routine employed by the fast solver. This array
c is initialized on the first call to CAPCN,
c or whenever ICFLAG=0 or 1. The array must not
c be altered between successive calls to the routine with
c different source functions when ICFLAG=2.
c
c IERROR integer constant
c An error flag on input variables. On return IERROR
c has the following meaning:
c
c IERROR=0 - normal status, no error detected.
c =1 - value of NX-1 or NY-1 is less than 8.
c =2 - illegal value for NX.
c =3 - illegal value for IP or JP.
c =4 - the boundary points in IP and JP are
c inconsistent with the array MASK.
c =5 - illegal value for ICFLAG.
c =6 - an illegal value of G detected.
c =7 - an error occurred in factoring the inverse
c capacitance matrix, i.e. the matrix is singular.
c This condition will occur if there the user has
c defined duplicate irregular boundary points.
c =8 - maximum number of iterations reached.
c =9 - nonconvergence of the iteration detected.
c =10 - CCON is less than or equal to zero.
c =11 - illegal value for DX or DY.
c
c Note: Error checking is normally only done on the first call
c to the routine. The user may force the routine to do error
c checking at any time by setting IERROR=-1 on entry.
c
c
c ********* NOTES ON CALLING ROUTINE CAPCN ********
c
c In applying the capacitance matrix method to the solution
c of a self-adjoint elliptic equation, routine CAPCN works in
c two stages. The first stage is a preprocessing one in which
c the arrays C and IPS are determined. The second stage is
c the actual solution of the problem for a given source
c function. The flag ICFLAG determines the operation of the
c routine:
c
c If the routine is called with ICFLAG=0, only
c preprocessing is done, the arrays C and IPC
c are determined, and CAPCN immediately returns to
c the calling program without solving. On return
c ICFLAG=2
c
c If the routine is called with ICFLAG=1,
c preprocessing is done, the arrays C and IPC
c are determined, and the solution for the
c given source function is determined. On return
c ICFLAG=2.
c
c If the routine is called with ICFLAG=2,
c the solution for the given source function
c is determined. This assumes that preprocessing
c has been done and that the C and IPS arrays
c are defined. On return ICFLAG=2.
c
c If the user wishes to solve the elliptic equation in a
c given domain repeatedly with different source functions,
c then the routine could be called with ICFLAG=1 initially.
c Afterwards this flag is set to 2 and need not be altered.
c If the user has stored the arrays C and IPS from a previous
c run, then these may be supplied as input data and the
c routine may be called initially with ICFLAG=2 to obtain
c the solution directly.
c
c Preprocessing need only be done once for a given domain
c geometry. The user-specified variables listed below must not
c change between successive calls to CAPCN with ICFLAG=2. If
c the user changes any one of the following variables:
c
c NX,
c NY,
c DX,
c DY,
c IP,
c JP,
c NP,
c
c then the preprocessing must be redone (ICFLAG=0 or 1). The
c arrays TRIGS, C, and IPS should also remain unchanged on
c successive calls with ICFLAG=2.
c
c
c ********* NOTES ON DEFINING THE SOLUTION DOMAIN *********
c
c The routine solves an elliptic equation in an irregular,
c polygonal domain. The domain is polygonal because its boundary
c consists of straight line segments parallel or diagonal to
c the discrete mesh. This domain is embedded in a rectangle
c dimensioned NX by NY. The user is free to position the irregular
c domain as desired within the embedding rectangle and the sides of
c the rectangle may be part of the boundary of the solution domain.
c
c The boundary of the irregular domain is described by a set of
c grid points which are classified as 'regular' boundary points
c and 'irregular' boundary points. A regular boundary point is one
c which lies on one of the sides of the embedding rectangle. An
c irregular boundary point is one which is not on a side of the
c embedding rectangle. The routine requires that the user provide
c the coordinates of the irregular boundary points through the
c one-dimensional arrays IP and JP, both dimensioned NP. The user
c must ensure that the specified coordinates of the irregular boundary
c points do not lie on a side of the embedding rectangle, i.e. the
c conditions 1 < IP(N) < NX and 1 < JP(N) < NY must hold for all N.
c Failure to ensure this will result in an error condition. In
c addition the routine requires that there be at least one irregular
c grid point ( NP greater than or equal to 1 ).
c
c Boundary values at the irregular boundary points are supplied
c through the one-dimensional array VP, while boundary values on
c the edge of the embedding rectangle are supplied through the BA,
c BB, BC, and BD arrays. For grid points on the edge of the embedding
c rectangle which are not regular boundary points (i.e. not part of
c the boundary of the solution domain), the corresponding boundary values
c specified through the BA, BB, BC, and BD arrays are arbitrary. For
c definiteness they should be set to zero.
c
c Routine CAPCN also requires that a MASK array, dimensioned
c NX by NY, be specified to allow the routine to identify interior
c points where the solution is to be obtained. The array elements of
c MASK have either a value of 1, for interior points, or 0, for all
c other points.
c
c Three examples follow. In each case the grid points over
c the rectangle are labelled according to the following
c convention: Interior Points - labelled 1,
c Regular Boundary Points - labelled 2,
c Irregular Boundary Points - labelled 3,
c Exterior Points - labelled 0.
c
c **** Example 1 - rectangle with a slot removed ****
c
c NX=17, NY=9, NP=10
c IP=(8,8,8,8,9,10,11,11,11,11)
c JP=(2,3,4,5,5, 5, 5, 4, 3, 2)
c
c 22222222222222222
c 21111111111111112
c 21111111111111112
c 21111111111111112
c 21111113333111112
c 21111113003111112
c 21111113003111112
c 21111113003111112
c 22222222002222222
c
c Boundary values for the problem are specified at the irregular
c boundary points through the VP array, and at the regular boundary
c points through the BA, BB, BC and BD arrays. The
c BC(I), I=1,NX array specifies boundary values on the lower edge
c (J=1) on the embedding rectangle. These values are arbitrary for
c I=9,10 and should be set to zero.
c
c The corresponding MASK array for Example 1: interior points set
c equal to one, all other points equal to zero.
c
c 00000000000000000
c 01111111111111110
c 01111111111111110
c 01111111111111110
c 01111110000111110
c 01111110000111110
c 01111110000111110
c 01111110000111110
c 00000000000000000
c
c **** Example 2 - rectangle with a triangular hole ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 22222222222222222
c 21111111111111112
c 21133333333333112
c 21113000000031112
c 21111300000311112
c 21111130003111112
c 21111113031111112
c 21111111311111112
c 22222222222222222
c
c The corresponding MASK array for Example 2
c
c 00000000000000000
c 01111111111111110
c 01100000000000110
c 01110000000001110
c 01111000000011110
c 01111100000111110
c 01111110001111110
c 01111111011111110
c 00000000000000000
c
c
c **** Example 3 - triangular region ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 00000000000000000
c 00000000000000000
c 00033333333333000
c 00003111111130000
c 00000311111300000
c 00000031113000000
c 00000003130000000
c 00000000300000000
c 00000000000000000
c
c There are no regular boundary points in this case, and so the
c values of the BA, BB, BC and BD arrays are arbitrary. They
c should be set to zero for definiteness.
c
c The corresponding MASK array for Example 3:
c
c 00000000000000000
c 00000000000000000
c 00000000000000000
c 00000111111100000
c 00000011111000000
c 00000001110000000
c 00000000100000000
c 00000000000000000
c 00000000000000000
c
c Note that the irregular boundary points are the same for
c examples 2 and 3. The MASK array allows the routine to
c identify the interior region in each case.
c
c Note also that, in example 3, repositioning the triangular region
c so that one of its sides coincides with one of the sides of the
c embedding rectangle would reduce the number of irregular grid
c points and reduce time for preprocessing and solving.
c
c *************************************************************
c
PARAMETER(MMAX=6,MMAX1=MMAX-1)
REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
REAL P(NX,NY),G(NX-1,NY-1),RHO
REAL C(NP,NP),CHARGE(NP),VP(NP)
REAL BA(NY),BB(NY),BC(NX),BD(NX),AN1(MMAX)
INTEGER IPS(NP),IP(NP),JP(NP),MASK(NX,NY)
INTEGER IFAX(10)
SAVE LFIRST,NX1,NY1,IFAX
DATA LFIRST/1/
c Initialize arrays, error checking is done on the first call
IF((LFIRST.EQ.1).OR.(ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
IERROR=0
CALL CHKV2(G,MASK,NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
NX1=NX-1
NY1=NY-1
CALL FAX(IFAX,NX1,4)
CALL FFTRIG(TRIGS,NX1,4)
LFIRST=0
END IF
c If required, execute the error checking routine
IF(IERROR.EQ.-1) THEN
IERROR=0
CALL CHKV2(G,MASK,NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
END IF
RHO=DY/DX
IF(ICFLAG.EQ.2) GOTO 5000
c preprocessing: generate and factor the Green's function matrix
IF((ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
ICFG=ICFLAG
CALL CPGN2(C,U,WORK,TRIGS,IFAX,IPS,0.,
# RHO,NX,NY,IP,JP,NP,IERROR)
ICFLAG=2
ELSE
IERROR=5
RETURN
END IF
IF(ICFG.EQ.0) RETURN
5000 CONTINUE
IF(CCON.LE.0.) THEN
IERROR=10
RETURN
END IF
NITT=1
DO 5 M=1,MMAX
AN1(M)=0.
5 CONTINUE
DO 10 I=1,NX
P(I,1)=BC(I)
P(I,NY)=BD(I)
10 CONTINUE
DO 20 J=2,NY1
P(1,J)=BA(J)
P(NX,J)=BB(J)
20 CONTINUE
1000 CONTINUE
c Form the right hand side
DO 35 J=2,NY1
DO 30 I=2,NX1
U(I,J)=(4.*DY*DY*F(I,J)-RHO*RHO*MASK(I,J)*(P(I+1,J)-P(I-1,J))*
# (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
# - MASK(I,J)*(P(I,J+1)-P(I,J-1))*
# (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
# (MASK(I,J)*(G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
# -(MASK(I,J)-1.))
30 CONTINUE
35 CONTINUE
c Modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the regular boundary
DO 40 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
40 CONTINUE
DO 50 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
50 CONTINUE
c Get the solution in the regular domain
CALL POISS(U,0.,RHO,WORK,TRIGS,IFAX,NX,NY)
c Obtain the correction to the rhs at the irregular boundary points
CDIR$ IVDEP
CVD$L NOSYNC
DO 60 N1=1,NP
CHARGE(N1)=VP(N1)-U(IP(N1),JP(N1))
60 CONTINUE
CALL SGESL(C,NP,NP,IPS,CHARGE,0)
c Reconstruct the rhs
DO 75 J=2,NY1
DO 70 I=2,NX1
U(I,J)=(4.*DY*DY*F(I,J)-RHO*RHO*MASK(I,J)*(P(I+1,J)-P(I-1,J))*
# (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
# - MASK(I,J)*(P(I,J+1)-P(I,J-1))*
# (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
# (MASK(I,J)*(G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
# -(MASK(I,J)-1.))
70 CONTINUE
75 CONTINUE
c Modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the regular boundary
DO 80 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
80 CONTINUE
DO 90 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
90 CONTINUE
c Modify the rhs at the irregular points
CDIR$ IVDEP
CVD$L NOSYNC
DO 110 N1=1,NP
U(IP(N1),JP(N1))=U(IP(N1),JP(N1))+CHARGE(N1)
110 CONTINUE
c Get the solution in the irregular domain
CALL POISS(U,0.,RHO,WORK,TRIGS,IFAX,NX,NY)
c Set boundary conditions on the regular and irregular boundaries
DO 120 I=1,NX
U(I,1)=BC(I)
U(I,NY)=BD(I)
120 CONTINUE
DO 130 J=2,NY1
U(1,J)=BA(J)
U(NX,J)=BB(J)
130 CONTINUE
CDIR$ IVDEP
CVD$L NOSYNC
DO 140 N1=1,NP
U(IP(N1),JP(N1))=VP(N1)
140 CONTINUE
c Check for convergence of the solution
ADEN=0.
ADIFF=0.
DL2N=0.
DENM=0.
CVD$ ASSOC (D)
DO 150 J=2,NY1
DO 160 I=2,NX1
ADEN=ADEN+MASK(I,J)*U(I,J)*U(I,J)
ADIFF=ADIFF+MASK(I,J)*(U(I,J)-P(I,J))*(U(I,J)-P(I,J))
FD=MASK(I,J)*0.5*((G(I,J)+G(I,J-1))*(U(I+1,J)-U(I,J))
# - (G(I-1,J)+G(I-1,J-1))*(U(I,J)-U(I-1,J)))/(DX*DX)
# +MASK(I,J)*0.5*((G(I,J)+G(I-1,J))*(U(I,J+1)-U(I,J))
# - (G(I,J-1)+G(I-1,J-1))*(U(I,J)-U(I,J-1)))/(DY*DY)
DL2N=DL2N+(FD-MASK(I,J)*F(I,J))*(FD-MASK(I,J)*F(I,J))
DENM=DENM+MASK(I,J)*F(I,J)*F(I,J)
160 CONTINUE
150 CONTINUE
IF(ADEN.NE.0.) THEN
ANORM=SQRT(ADIFF/ADEN)
ELSE
ANORM=SQRT(ADIFF)
END IF
IF(DENM.NE.0.) THEN
RRESD=SQRT(DL2N/DENM)
ELSE IF(ADEN.NE.0.) THEN
RRESD=SQRT(DL2N/ADEN)
ELSE
RRESD=SQRT(DL2N)
END IF
c If required, print out the iteration by iteration convergence
c of the solver.
IF(IWRITE.NE.0) THEN
WRITE(IWRITE,101) NITT,ANORM,RRESD
END IF
101 FORMAT(' ','ITERATION ',I4,' ANORM= ',1PG11.5,' RESIDUAL= ',
# 1PG11.5)
IF((ANORM.LE.CCON).AND.(RRESD.LE.CCON)) RETURN
IF(ITMAX.GT.0) THEN
DO 180 M=2,MMAX
MM=MMAX-M+2
AN1(MM)=AN1(MM-1)
180 CONTINUE
AN1(1)=SQRT(ADIFF)
RTSUM1=0.
DO 190 M=1,MMAX1
RAT1=0.
IF(AN1(M+1).NE.0.) THEN
RAT1=AN1(M)/AN1(M+1)
IF(RAT1.GT.1.) RTSUM1=RTSUM1+1.
END IF
190 CONTINUE
IF((RTSUM1.GT.2.99)) THEN
IERROR=9
RETURN
END IF
END IF
IF(NITT.GE.ABS(ITMAX)) THEN
IERROR=8
RETURN
END IF
DO 500 J=1,NY
DO 510 I=1,NX
P(I,J)=U(I,J)
510 CONTINUE
500 CONTINUE
NITT=NITT+1
GO TO 1000
END
c Routine to generate and factor the inverse
c capacitance matrix
SUBROUTINE CPGN2(C,R,WORK,TRIGS,IFAX,IPS,HELM,
# RHO,NX,NY,IP,JP,NP,IERROR)
REAL C(NP,NP),R(NX,NY),WORK(NX,NY),TRIGS(*)
REAL HELM,RHO
INTEGER IP(NP),JP(NP),IPS(NP),IFAX(10)
CVD$ NOSELECT
DO 10 N1=1,NP
DO 30 J=1,NY
DO 20 I=1,NX
R(I,J)=0.
20 CONTINUE
30 CONTINUE
R(IP(N1),JP(N1))=1.
CALL POISS(R,HELM,RHO,WORK,TRIGS,IFAX,NX,NY)
DO 40 N2=1,NP
C(N2,N1)=R(IP(N2),JP(N2))
40 CONTINUE
10 CONTINUE
c Linpack routine to factor the matrix
CALL SGEFA(C,NP,NP,IPS,IER)
IF(IER.NE.0) THEN
IERROR=7
END IF
RETURN
END
c
c Routine to check input variables, MASK array and G array
c
SUBROUTINE CHKV2 (G,MASK,NX,NY,IP,JP,NP,DX,DY,IERROR)
REAL G(NX-1,NY-1),DX,DY
INTEGER MASK(NX,NY),IP(NP), JP(NP)
DATA MAXCK/40/
NX1=NX-1
NY1=NY-1
c Check if the values of NX and NY are too small
IF (NX1.LT.8.OR.NY1.LT.8) THEN
IERROR=1
RETURN
END IF
c Check if the value of NX is admissible
LR=0
LP=0
LQ=0
N1=NX-1
DO 100 IT=1,MAXCK
MX=MOD(N1,2)
IF(MX.NE.0) GOTO 110
LP=LP+1
N1=N1/2
IF(N1.EQ.1) GOTO 401
100 CONTINUE
110 CONTINUE
IF(LP.EQ.0) THEN
IERROR=2
RETURN
END IF
DO 200 IT=1,MAXCK
MX=MOD(N1,3)
IF(MX.NE.0) GOTO 220
LQ=LQ+1
N1=N1/3
IF(N1.EQ.1) GOTO 401
200 CONTINUE
220 CONTINUE
DO 300 IT=1,MAXCK
MX=MOD(N1,5)
IF(MX.NE.0) GOTO 401
LR=LR+1
N1=N1/5
IF(N1.EQ.1) GOTO 401
300 CONTINUE
401 CONTINUE
NXX=(2**LP)*(3**LQ)*(5**LR) + 1
IF(NXX.NE.NX) THEN
IERROR=2
RETURN
END IF
c Check values of the irregular boundary points
DO 500 IN=1,NP
IF(IP(IN).LT.2.OR.IP(IN).GT.NX-1) IERROR=3
IF(JP(IN).LT.2.OR.JP(IN).GT.NY-1) IERROR=3
500 CONTINUE
IF(IERROR.NE.0) RETURN
c Check if the mask array is compatible with the irregular
c boundary points
CDIR$ IVDEP
CVD$L NOSYNC
DO 600 IN=1,NP
IF(MASK(IP(IN),JP(IN)).NE.0) IERROR=4
600 CONTINUE
IF(IERROR.NE.0) RETURN
c Check if the diffusion function is everywhere positive over
c the interior of the irregular domain
CVD$ SELECT
DO 800 J=2,NY1
DO 700 I=2,NX1
IF(MASK(I,J).EQ.1) THEN
IF((G(I,J).LE.0.).OR.(G(I-1,J).LE.0.).OR.
# (G(I,J-1).LE.0.).OR.(G(I-1,J-1).LE.0.)) THEN
IERROR=6
END IF
END IF
700 CONTINUE
800 CONTINUE
IF((DX.LE.0.).OR.(DY.LE.0.)) IERROR=11
RETURN
END
c***DECK : dcapcn.f
SUBROUTINE DCAPCN(F,G,U,P,WORK,MASK,TRIGS,NX,NY,DX,DY,
# C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
# BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
c Programmer: P. Cummins
c Version 1.0
c
c *********************************************************************
c
c DCAPCN is a double precision routine to solve a second-order
c accurate finite difference approximation to the self-adjoint
c elliptic equation
c
c (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = F(X,Y)
c
c with Dirichlet boundary conditions in irregular polygonal
c two-dimensional domains in Cartesian coordinates using the
c capacitance matrix method.
c
c The diffusion function, G, is staggered with respect
c to the solution and the source function on the grid.
c
c | | |
c --U(I-1,J+1)-------U(I,J+1)----U(I+1,J+1)--
c | | |
c | | |
c | G(I-1,J) | G(I,J) |
c | | |
c | | |
c --U(I-1,J)---------U(I,J)------U(I+1,J)--
c | | |
c | | |
c | G(I-1,J-1) | G(I,J-1) |
c | | |
c | | |
c --U(I-1,J-1)-------U(I,J-1)----U(I-1,J+1)--
c | | |
c
c ******************* ON ENTRY ********************
c
c NX integer constant
c The number of mesh points in the X direction
c Note that the value of NX is restricted such
c that (NX-1) is of the form (2**I * 3**J * 5**K),
c where I is an integer greater than or equal to
c one, and J and K are integers greater than or
c equal to zero. In addition, NX-1 must be greater
c than or equal to 8.
c
c NY integer constant
c The number of mesh points in the Y direction. NY-1
c must be greater than or equal to 8.
c
c NP integer constant
c The number of irregular boundary points. NP must be
c greater than or equal to one.
c
c F double precision array (NX,NY)
c The source function on the right hand side. The value
c of F is arbitrary at points in the rectangle which are
c outside the irregular solution domain. This array is
c unchanged on exit.
c
c P double precision array (NX,NY)
c The initial estimate to the solution; if the user does
c not have an initial estimate, this array should be set
c to zero. This array is overwritten on exit.
c
c G double precision array (NX-1,NY-1)
c The array containing values of the diffusion function.
c Note that this array is staggered with respect to the
c mesh on which the source function and the solution are
c defined. To ensure a well posed problem, the condition
c G(I,J) > 0 must hold for all points within the interior
c of the irregular domain. The value of G is arbitrary at
c points in the rectangle which are outside the irregular
c solution domain. This array is unchanged on exit.
c
c MASK integer array (NX,NY)
c This array is used by the routine to identify the mesh points
c within the embedding rectangle that lie inside the irregular
c (polygonal) domain. The user must set MASK(I,J)=1 for all
c (I,J) values which are inside the irregular domain and
c MASK(I,J)=0 for all other points, including the boundary
c points. The mesh of grid points defined by the MASK array
c coincides with that of the source function and the solution.
c This array is unchanged on exit.
c
c Additional notes on defining the MASK array are given below.
c
c WORK double precision array (NX,NY)
c A work space array.
c
c CHARGE double precision array (NP)
c A work space array.
c
c DX double precision constant
c The discretization interval in the X direction.
c The condition DX > 0 must hold.
c
c DY double precision constant
c The discretization interval in the Y direction.
c The condition DY > 0 must hold.
c
c BA double precision array (NY)
c The boundary values on the regular boundary
c at (I=1, J=2,NY-1).
c
c BB double precision array (NY)
c The boundary values on the regular boundary
c at (I=NX, J=2,NY-1).
c
c BC double precision array (NX)
c The boundary values on the regular boundary
c at (I=1,NX, J=1).
c
c BD double precision array (NX)
c The boundary values on the regular boundary
c at (I=1,NX, J=NY).
c
c VP double precision array (NP)
c The boundary values on the irregular boundary.
c
c IP integer array (NP)
c The coordinates in the X direction of the
c irregular boundary points. The condition,
c 1 < IP(N) < NX, must hold for all N.
c
c JP integer array (NP)
c The coordinates in the Y direction of the
c irregular boundary points. The condition,
c 1 < JP(N) < NY, must hold for all N.
c
c Note: The order in which the irregular boundary grid points, defined
c by the (IP(N),JP(N)) pairs, are specified is immaterial.
c However, these grid points must all be distinct; there must
c not be a repeated (IP(N),JP(N)) pair. Failure to ensure this
c will result in the error condition IERROR=7 (see below).
c
c Additional notes on defining the IP and JP arrays are
c given below.
c
c ICFLAG integer constant
c A flag: ICFLAG=0 - do preprocessing only and compute
c the matrix C and indices IPS.
c
c ICFLAG=1 - do preprocessing to compute the
c matrix C and indices IPS, and solve
c the elliptic problem.
c
c ICFLAG=2 - solve the elliptic problem only; this
c assumes that preprocessing has been
c done and that C and IPS are defined.
c
c Any other value for ICFLAG gives IERROR=5 and no attempt
c is made to do preprocessing or to solve. On return,
c ICFLAG=2.
c
c CCON double precision constant
c The convergence criterion. The routine will proceed with
c the iteration to reach a solution until the following
c two conditions are satisfied
c
c (1) || U(k)-U(k-1) || / ||U(k)|| <= CCON
c
c (2) || L(U(k))-F || / ||F|| <= CCON,
c
c or until the number of iterations has reached the maximum
c allowable, given by parameter ITMAX. Note that ||U(k)||
c refers to the L2 norm of the kth iterate of U and that L
c refers to the discrete elliptic operator. If ||F||=0, then
c convergence is determined from condition (1) alone. CCON
c must be greater than zero or the routine returns without
c attempting a solution and IERROR=10.
c
c The minimum possible value of CCON is limited by
c the accumulation of round-off error and depends
c mainly on the machine precision and the size of
c the grid. The following table of approximate minimum
c values of CCON was constructed from solution of the
c elliptic equation, in a unit square domain with the
c source and diffusion functions given by a uniformly
c distributed random number, on an Alliant FX/40. This table
c is intended as a guide to the user; the actual minimum
c attainable value of CCON depends on the particular problem
c and on the computer used.
c
c (NX,NY) Single Precision (32 bit) Double Precision (64 bit)
c
c (17,17) 4.E-6 2.E-14
c
c (33,33) 3.E-5 5.E-14
c
c (65,65) 2.E-4 4.E-13
c
c (129,129) 2.E-3 4.E-12
c
c (257,257) 2.E-2 4.E-11
c
c (513,513) ----- 2.E-10
c
c (1025,1025) ----- 2.E-09
c
c Under normal operation, the solver will iterate until
c the convergence conditions are satisfied and then
c return. The convergence of the routine is monitored
c at each iteration and if a lack of convergence is
c detected, then the routine will return with
c IERROR=9. Specifically the ratio D(k+1)/D(k), where
c D(k)=|| U(k+1)-U(k) ||, is monitored for cases when
c it is greater than one. Should this occur repeatedly
c (three times over six iterations), then lack of
c convergence is assumed and the routine returns to
c the calling program with IERROR=9. This feature is
c invoked automatically, but may be disabled by
c setting ITMAX (see below) to a negative value.
c
c One situation which may lead to a lack of convergence
c arises if CCON is set to a value which is too small
c for a given machine precision. In such cases, the routine
c initially converges to the solution for a number of
c iterations and then fails to converge further with
c additional iterations. This may be monitored in the
c output to unit IWRITE (see below). The remedy is to
c increase the value of CCON or to use the double precision
c routine DCAPCN.
c
c IWRITE integer constant
c Logical unit number to which the routine will write out
c the iteration number, k, and the convergence conditions,
c
c ANORM = || U(k)-U(k-1)|| / ||U(k)||
c
c RESIDUAL = || L(U(k))-F || / ||F||.
c
c This is useful for examining the iteration by iteration
c convergence of the solver. IWRITE=0 suppresses this output.
c
c ITMAX integer constant
c The absolute value, |ITMAX|, gives the maximum number of
c iteration permitted. If the number of iterations reaches
c this maximum without obtaining a solution to the required
c accuracy, then IERROR=4 on output. The routine accepts
c positive or negative values for ITMAX, with negative
c values suppressing the iterative check on convergence
c described above under parameter CCON.
c
c ******************* ON RETURN ********************
c
c U double precision array (NX,NY)
c The solution in the irregular domain, with
c U(IP(N),JP(N)) = VP(N) for N=1,NP on the irregular
c portion of the boundary and also with
c U(1,J)=BA(J) for J=2,NY-1
c U(NX,J)=BB(J) for J=2,NY-1
c U(I,1)=BC(I) for I=1,NX
c U(I,NY)=BD(I) for I=1,NX.
c
c NITT integer constant
c The number of iterations required to get the
c solution within the prescribed margin of error
c given by CCON.
c
c C double precision array (NP,NP)
c The inverse capacitance matrix. This array must not be
c be altered on successive calls to the routine with
c different source functions when ICFLAG=2.
c
c IPS integer array (NP)
c An integer vector of pivot indices. This array must not
c not be alter on successive calls to the routine with
c different source functions when ICFLAG=2.
c
c TRIGS double precision array (2*(NX-1))
c An array of coefficients required by the FFT
c routine employed by the fast solver. This array
c is initialized on the first call to DCAPCN,
c or whenever ICFLAG=0 or 1. The array must not
c be altered between successive calls to the routine with
c different source functions when ICFLAG=2.
c
c IERROR integer constant
c An error flag on input variables. On return IERROR
c has the following meaning:
c
c IERROR=0 - normal status, no error detected.
c =1 - value of NX-1 or NY-1 is less than 8.
c =2 - illegal value for NX.
c =3 - illegal value for IP or JP.
c =4 - the boundary points in IP and JP are
c inconsistent with the array MASK.
c =5 - illegal value for ICFLAG.
c =6 - an illegal value of G detected.
c =7 - an error occurred in factoring the inverse
c capacitance matrix, i.e. the matrix is singular.
c This condition will occur if there the user has
c defined duplicate irregular boundary points.
c =8 - maximum number of iterations reached.
c =9 - nonconvergence of the iteration detected.
c =10 - CCON is less than or equal to zero.
c =11 - illegal value for DX or DY.
c
c Note: Error checking is normally only done on the first call
c to the routine. The user may force the routine to do error
c checking at any time by setting IERROR=-1 on entry.
c
c
c ********* NOTES ON CALLING ROUTINE DCAPCN ********
c
c In applying the capacitance matrix method to the solution
c of a self-adjoint elliptic equation, routine DCAPCN works in
c two stages. The first stage is a preprocessing one in which
c the arrays C and IPS are determined. The second stage is
c the actual solution of the problem for a given source
c function. The flag ICFLAG determines the operation of the
c routine:
c
c If the routine is called with ICFLAG=0, only
c preprocessing is done, the arrays C and IPC
c are determined, and DCAPCN immediately returns to
c the calling program without solving. On return
c ICFLAG=2
c
c If the routine is called with ICFLAG=1,
c preprocessing is done, the arrays C and IPC
c are determined, and the solution for the
c given source function is determined. On return
c ICFLAG=2.
c
c If the routine is called with ICFLAG=2,
c the solution for the given source function
c is determined. This assumes that preprocessing
c has been done and that the C and IPS arrays
c are defined. On return ICFLAG=2.
c
c If the user wishes to solve the elliptic equation in a
c given domain repeatedly with different source functions,
c then the routine could be called with ICFLAG=1 initially.
c Afterwards this flag is set to 2 and need not be altered.
c If the user has stored the arrays C and IPS from a previous
c run, then these may be supplied as input data and the
c routine may be called initially with ICFLAG=2 to obtain
c the solution directly.
c
c Preprocessing need only be done once for a given domain
c geometry. The user-specified variables listed below must not
c change between successive calls to DCAPCN with ICFLAG=2. If
c the user changes any one of the following variables:
c
c NX,
c NY,
c DX,
c DY,
c IP,
c JP,
c NP,
c
c then the preprocessing must be redone (ICFLAG=0 or 1). The
c arrays TRIGS, C, and IPS should also remain unchanged on
c successive calls with ICFLAG=2.
c
c
c ********* NOTES ON DEFINING THE SOLUTION DOMAIN *********
c
c The routine solves an elliptic equation in an irregular,
c polygonal domain. The domain is polygonal because its boundary
c consists of straight line segments parallel or diagonal to
c the discrete mesh. This domain is embedded in a rectangle
c dimensioned NX by NY. The user is free to position the irregular
c domain as desired within the embedding rectangle and the sides of
c the rectangle may be part of the boundary of the solution domain.
c
c The boundary of the irregular domain is described by a set of
c grid points which are classified as 'regular' boundary points
c and 'irregular' boundary points. A regular boundary point is one
c which lies on one of the sides of the embedding rectangle. An
c irregular boundary point is one which is not on a side of the
c embedding rectangle. The routine requires that the user provide
c the coordinates of the irregular boundary points through the
c one-dimensional arrays IP and JP, both dimensioned NP. The user
c must ensure that the specified coordinates of the irregular boundary
c points do not lie on a side of the embedding rectangle, i.e. the
c conditions 1 < IP(N) < NX and 1 < JP(N) < NY must hold for all N.
c Failure to ensure this will result in an error condition. In
c addition the routine requires that there be at least one irregular
c grid point ( NP greater than or equal to 1 ).
c
c Boundary values at the irregular boundary points are supplied
c through the one-dimensional array VP, while boundary values on
c the edge of the embedding rectangle are supplied through the BA,
c BB, BC, and BD arrays. For grid points on the edge of the embedding
c rectangle which are not regular boundary points (i.e. not part of
c the boundary of the solution domain), the corresponding boundary values
c specified through the BA, BB, BC, and BD arrays are arbitrary. For
c definiteness they should be set to zero.
c
c Routine DCAPCN also requires that a MASK array, dimensioned
c NX by NY, be specified to allow the routine to identify interior
c points where the solution is to be obtained. The array elements of
c MASK have either a value of 1, for interior points, or 0, for all
c other points.
c
c Three examples follow. In each case the grid points over
c the rectangle are labelled according to the following
c convention: Interior Points - labelled 1,
c Regular Boundary Points - labelled 2,
c Irregular Boundary Points - labelled 3,
c Exterior Points - labelled 0.
c
c **** Example 1 - rectangle with a slot removed ****
c
c NX=17, NY=9, NP=10
c IP=(8,8,8,8,9,10,11,11,11,11)
c JP=(2,3,4,5,5, 5, 5, 4, 3, 2)
c
c 22222222222222222
c 21111111111111112
c 21111111111111112
c 21111111111111112
c 21111113333111112
c 21111113003111112
c 21111113003111112
c 21111113003111112
c 22222222002222222
c
c Boundary values for the problem are specified at the irregular
c boundary points through the VP array, and at the regular boundary
c points through the BA, BB, BC and BD arrays. The
c BC(I), I=1,NX array specifies boundary values on the lower edge
c (J=1) on the embedding rectangle. These values are arbitrary for
c I=9,10 and should be set to zero.
c
c The corresponding MASK array for Example 1: interior points set
c equal to one, all other points equal to zero.
c
c 00000000000000000
c 01111111111111110
c 01111111111111110
c 01111111111111110
c 01111110000111110
c 01111110000111110
c 01111110000111110
c 01111110000111110
c 00000000000000000
c
c **** Example 2 - rectangle with a triangular hole ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 22222222222222222
c 21111111111111112
c 21133333333333112
c 21113000000031112
c 21111300000311112
c 21111130003111112
c 21111113031111112
c 21111111311111112
c 22222222222222222
c
c The corresponding MASK array for Example 2
c
c 00000000000000000
c 01111111111111110
c 01100000000000110
c 01110000000001110
c 01111000000011110
c 01111100000111110
c 01111110001111110
c 01111111011111110
c 00000000000000000
c
c
c **** Example 3 - triangular region ****
c
c NX=17, NY=9, NP=20
c IP=(9,8,7,6,5,4,5,6,7,8,9,10,11,12,13,14,13,12,11,10)
c JP=(2,3,4,5,6,7,7,7,7,7,7, 7, 7, 7, 7, 7, 6, 5, 4, 3)
c
c 00000000000000000
c 00000000000000000
c 00033333333333000
c 00003111111130000
c 00000311111300000
c 00000031113000000
c 00000003130000000
c 00000000300000000
c 00000000000000000
c
c There are no regular boundary points in this case, and so the
c values of the BA, BB, BC and BD arrays are arbitrary. They
c should be set to zero for definiteness.
c
c The corresponding MASK array for Example 3:
c
c 00000000000000000
c 00000000000000000
c 00000000000000000
c 00000111111100000
c 00000011111000000
c 00000001110000000
c 00000000100000000
c 00000000000000000
c 00000000000000000
c
c Note that the irregular boundary points are the same for
c examples 2 and 3. The MASK array allows the routine to
c identify the interior region in each case.
c
c Note also that, in example 3, repositioning the triangular region
c so that one of its sides coincides with one of the sides of the
c embedding rectangle would reduce the number of irregular grid
c points and reduce time for preprocessing and solving.
c
c *************************************************************
c
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
PARAMETER(MMAX=6,MMAX1=MMAX-1)
DOUBLE PRECISION F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(*)
DOUBLE PRECISION P(NX,NY),G(NX-1,NY-1),RHO
DOUBLE PRECISION C(NP,NP),CHARGE(NP),VP(NP)
DOUBLE PRECISION BA(NY),BB(NY),BC(NX),BD(NX),AN1(MMAX)
INTEGER IPS(NP),IP(NP),JP(NP),MASK(NX,NY)
INTEGER IFAX(10)
SAVE LFIRST,NX1,NY1,IFAX
DATA LFIRST/1/
c Initialize arrays, error checking is done on the first call
IF((LFIRST.EQ.1).OR.(ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
IERROR=0
CALL DCHKV2(G,MASK,NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
NX1=NX-1
NY1=NY-1
CALL DFAX(IFAX,NX1,4)
CALL DFFTRG(TRIGS,NX1,4)
LFIRST=0
END IF
c If required, execute the error checking routine
IF(IERROR.EQ.-1) THEN
IERROR=0
CALL DCHKV2(G,MASK,NX,NY,IP,JP,NP,DX,DY,IERROR)
IF(IERROR.NE.0) RETURN
END IF
RHO=DY/DX
IF(ICFLAG.EQ.2) GOTO 5000
c preprocessing: generate and factor the Green's function matrix
IF((ICFLAG.EQ.0).OR.(ICFLAG.EQ.1)) THEN
ICFG=ICFLAG
CALL DCPGN2(C,U,WORK,TRIGS,IFAX,IPS,0.D0,
# RHO,NX,NY,IP,JP,NP,IERROR)
ICFLAG=2
ELSE
IERROR=5
RETURN
END IF
IF(ICFG.EQ.0) RETURN
5000 CONTINUE
IF(CCON.LE.0.D0) THEN
IERROR=10
RETURN
END IF
NITT=1
DO 5 M=1,MMAX
AN1(M)=0.D0
5 CONTINUE
DO 10 I=1,NX
P(I,1)=BC(I)
P(I,NY)=BD(I)
10 CONTINUE
DO 20 J=2,NY1
P(1,J)=BA(J)
P(NX,J)=BB(J)
20 CONTINUE
1000 CONTINUE
c Form the right hand side
DO 35 J=2,NY1
DO 30 I=2,NX1
U(I,J)=(4.D0*DY*DY*F(I,J)-RHO*RHO*MASK(I,J)*(P(I+1,J)-P(I-1,J))*
# (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
# - MASK(I,J)*(P(I,J+1)-P(I,J-1))*
# (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
# (MASK(I,J)*(G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
# -(MASK(I,J)-1.D0))
30 CONTINUE
35 CONTINUE
c Modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the regular boundary
DO 40 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
40 CONTINUE
DO 50 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
50 CONTINUE
c Get the solution in the regular domain
CALL DPOISS(U,0.D0,RHO,WORK,TRIGS,IFAX,NX,NY)
c Obtain the correction to the rhs at the irregular boundary points
CDIR$ IVDEP
CVD$L NOSYNC
DO 60 N1=1,NP
CHARGE(N1)=VP(N1)-U(IP(N1),JP(N1))
60 CONTINUE
CALL DGESL(C,NP,NP,IPS,CHARGE,0)
c Reconstruct the rhs
DO 75 J=2,NY1
DO 70 I=2,NX1
U(I,J)=(4.D0*DY*DY*F(I,J)-RHO*RHO*MASK(I,J)*(P(I+1,J)-P(I-1,J))*
# (G(I,J)+G(I,J-1)-G(I-1,J)-G(I-1,J-1))
# - MASK(I,J)*(P(I,J+1)-P(I,J-1))*
# (G(I,J)+G(I-1,J)-G(I,J-1)-G(I-1,J-1))) /
# (MASK(I,J)*(G(I,J)+G(I-1,J)+G(I,J-1)+G(I-1,J-1))
# -(MASK(I,J)-1.D0))
70 CONTINUE
75 CONTINUE
c Modify the first interior pt of the rhs to allow for
c nonzero boundary conditions along the regular boundary
DO 80 I=2,NX1
U(I,2)=U(I,2)-BC(I)
U(I,NY1)=U(I,NY1)-BD(I)
80 CONTINUE
DO 90 J=2,NY1
U(2,J)=U(2,J)-RHO*RHO*BA(J)
U(NX1,J)=U(NX1,J)-RHO*RHO*BB(J)
90 CONTINUE
c Modify the rhs at the irregular points
CDIR$ IVDEP
CVD$L NOSYNC
DO 110 N1=1,NP
U(IP(N1),JP(N1))=U(IP(N1),JP(N1))+CHARGE(N1)
110 CONTINUE
c Get the solution in the irregular domain
CALL DPOISS(U,0.D0,RHO,WORK,TRIGS,IFAX,NX,NY)
c Set boundary conditions on the regular and irregular boundaries
DO 120 I=1,NX
U(I,1)=BC(I)
U(I,NY)=BD(I)
120 CONTINUE
DO 130 J=2,NY1
U(1,J)=BA(J)
U(NX,J)=BB(J)
130 CONTINUE
CDIR$ IVDEP
CVD$L NOSYNC
DO 140 N1=1,NP
U(IP(N1),JP(N1))=VP(N1)
140 CONTINUE
c Check for convergence of the solution
ADEN=0.D0
ADIFF=0.D0
DL2N=0.D0
DENM=0.D0
CVD$ ASSOC (D)
DO 150 J=2,NY1
DO 160 I=2,NX1
ADEN=ADEN+MASK(I,J)*U(I,J)*U(I,J)
ADIFF=ADIFF+MASK(I,J)*(U(I,J)-P(I,J))*(U(I,J)-P(I,J))
FD=MASK(I,J)*0.5D0*((G(I,J)+G(I,J-1))*(U(I+1,J)-U(I,J))
# - (G(I-1,J)+G(I-1,J-1))*(U(I,J)-U(I-1,J)))/(DX*DX)
# +MASK(I,J)*0.5D0*((G(I,J)+G(I-1,J))*(U(I,J+1)-U(I,J))
# - (G(I,J-1)+G(I-1,J-1))*(U(I,J)-U(I,J-1)))/(DY*DY)
DL2N=DL2N+(FD-MASK(I,J)*F(I,J))*(FD-MASK(I,J)*F(I,J))
DENM=DENM+MASK(I,J)*F(I,J)*F(I,J)
160 CONTINUE
150 CONTINUE
IF(ADEN.NE.0.D0) THEN
ANORM=SQRT(ADIFF/ADEN)
ELSE
ANORM=SQRT(ADIFF)
END IF
IF(DENM.NE.0.D0) THEN
RRESD=SQRT(DL2N/DENM)
ELSE IF(ADEN.NE.0.D0) THEN
RRESD=SQRT(DL2N/ADEN)
ELSE
RRESD=SQRT(DL2N)
END IF
c If required, print out the iteration by iteration convergence
c of the solver.
IF(IWRITE.NE.0) THEN
WRITE(IWRITE,101) NITT,ANORM,RRESD
END IF
101 FORMAT(' ','ITERATION ',I4,' ANORM= ',1PG11.5,' RESIDUAL= ',
# 1PG11.5)
IF((ANORM.LE.CCON).AND.(RRESD.LE.CCON)) RETURN
IF(ITMAX.GT.0) THEN
DO 180 M=2,MMAX
MM=MMAX-M+2
AN1(MM)=AN1(MM-1)
180 CONTINUE
AN1(1)=SQRT(ADIFF)
RTSUM1=0.
DO 190 M=1,MMAX1
RAT1=0.
IF(AN1(M+1).NE.0.) THEN
RAT1=AN1(M)/AN1(M+1)
IF(RAT1.GT.1.D0) RTSUM1=RTSUM1+1.D0
END IF
190 CONTINUE
IF((RTSUM1.GT.2.99D0)) THEN
IERROR=9
RETURN
END IF
END IF
IF(NITT.GE.ABS(ITMAX)) THEN
IERROR=8
RETURN
END IF
DO 500 J=1,NY
DO 510 I=1,NX
P(I,J)=U(I,J)
510 CONTINUE
500 CONTINUE
NITT=NITT+1
GO TO 1000
END
c Routine to generate and factor the inverse
c capacitance matrix
SUBROUTINE DCPGN2(C,R,WORK,TRIGS,IFAX,IPS,HELM,
# RHO,NX,NY,IP,JP,NP,IERROR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION C(NP,NP),R(NX,NY),WORK(NX,NY),TRIGS(*)
DOUBLE PRECISION HELM,RHO
INTEGER IP(NP),JP(NP),IPS(NP),IFAX(10)
CVD$ NOSELECT
DO 10 N1=1,NP
DO 30 J=1,NY
DO 20 I=1,NX
R(I,J)=0.D0
20 CONTINUE
30 CONTINUE
R(IP(N1),JP(N1))=1.D0
CALL DPOISS(R,HELM,RHO,WORK,TRIGS,IFAX,NX,NY)
DO 40 N2=1,NP
C(N2,N1)=R(IP(N2),JP(N2))
40 CONTINUE
10 CONTINUE
c Linpack routine to factor the matrix
CALL DGEFA(C,NP,NP,IPS,IER)
IF(IER.NE.0) THEN
IERROR=7
END IF
RETURN
END
c
c Routine to check input variables, MASK array and G array
c
SUBROUTINE DCHKV2(G,MASK,NX,NY,IP,JP,NP,DX,DY,IERROR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION G(NX-1,NY-1),DX,DY
INTEGER MASK(NX,NY),IP(NP), JP(NP)
DATA MAXCK/40/
NX1=NX-1
NY1=NY-1
c Check if the values of NX and NY are too small
IF (NX1.LT.8.OR.NY1.LT.8) THEN
IERROR=1
RETURN
END IF
c Check if the value of NX is admissible
LR=0
LP=0
LQ=0
N1=NX-1
DO 100 IT=1,MAXCK
MX=MOD(N1,2)
IF(MX.NE.0) GOTO 110
LP=LP+1
N1=N1/2
IF(N1.EQ.1) GOTO 401
100 CONTINUE
110 CONTINUE
IF(LP.EQ.0) THEN
IERROR=2
RETURN
END IF
DO 200 IT=1,MAXCK
MX=MOD(N1,3)
IF(MX.NE.0) GOTO 220
LQ=LQ+1
N1=N1/3
IF(N1.EQ.1) GOTO 401
200 CONTINUE
220 CONTINUE
DO 300 IT=1,MAXCK
MX=MOD(N1,5)
IF(MX.NE.0) GOTO 401
LR=LR+1
N1=N1/5
IF(N1.EQ.1) GOTO 401
300 CONTINUE
401 CONTINUE
NXX=(2**LP)*(3**LQ)*(5**LR) + 1
IF(NXX.NE.NX) THEN
IERROR=2
RETURN
END IF
c Check values of the irregular boundary points
DO 500 IN=1,NP
IF(IP(IN).LT.2.OR.IP(IN).GT.NX-1) IERROR=3
IF(JP(IN).LT.2.OR.JP(IN).GT.NY-1) IERROR=3
500 CONTINUE
IF(IERROR.NE.0) RETURN
c Check if the mask array is compatible with the irregular
c boundary points
CDIR$ IVDEP
CVD$L NOSYNC
DO 600 IN=1,NP
IF(MASK(IP(IN),JP(IN)).NE.0) IERROR=4
600 CONTINUE
IF(IERROR.NE.0) RETURN
c Check if the diffusion function is everywhere positive over
c the interior of the irregular domain
CVD$ SELECT
DO 800 J=2,NY1
DO 700 I=2,NX1
IF(MASK(I,J).EQ.1) THEN
IF((G(I,J).LE.0.D0).OR.(G(I-1,J).LE.0.D0).OR.
# (G(I,J-1).LE.0.D0).OR.(G(I-1,J-1).LE.0.D0)) THEN
IERROR=6
END IF
END IF
700 CONTINUE
800 CONTINUE
IF((DX.LE.0.D0).OR.(DY.LE.0.D0)) IERROR=11
RETURN
END
C-END-OF-FILE
cat > testcapc.f << C-END-OF-FILE
c --------------------------------------------------------------
c
c Program TESTCAPC - Example of CAPC usage
c
c --------------------------------------------------------------
c
c This program illustrates the use of routine CAPC, a solver
c for the Helmholtz equation:
c (d/dX)(dU/dX) + (d/dY)(dU/dY) -H*U = Q(X,Y),
c in irregular (polygonal) two-dimensional domains with
c Dirichlet boundary conditions.
c
c The solution is obtain in the region defined by the equation
c |X-0.5| + |Y-0.5| <= 0.375, which is embedded in the unit
c square. The source function is Q(X,Y)=COS(PI*X)*SIN(2*PI*Y)
c and the Helmholtz coefficient H=1. The boundary
c conditions are U(X,Y)=0 on |X-0.5| + |Y-0.5| = 0.375. The
c coordinates of the boundary are specified in the BLOCK DATA
c routine at the end.
c
c The solution produced by solver is tested by applying the discrete
c Helmholtz operator to the computed solution and examining the
c difference (the residual) between this and the specified source
c function over the domain.
c
c --------------
c Declarations
c --------------
c
PARAMETER(NP=24,NX=17,NY=17,NX1=NX-1,NY1=NY-1,NTRG=2*(NX-1))
REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(NTRG)
REAL C(NP,NP),CHARGE(NP),VP(NP)
REAL BA(NY),BB(NY),BC(NX),BD(NX)
INTEGER IPS(NP)
COMMON /BNDY/ IP(NP),JP(NP)
c
c --------------------------------------
c Set the grid intervals and constants
c --------------------------------------
c
DX=1./REAL(NX1)
DY=1./REAL(NY1)
PI=4.*ATAN(1.)
HELM=1.
c
c -------------------------------------------
c Define the right-hand side source function
c -------------------------------------------
c
DO 60 J=2,NY1
DO 50 I=2,NX1
F(I,J)=COS(PI*REAL(I-1)*DX)*SIN(2.*PI*REAL(J-1)*DY)
50 CONTINUE
60 CONTINUE
c
c ------------------------------------------------------
c For definiteness, set the arrays for boundary values
c on the sides of the embedding rectangle to zero.
c ------------------------------------------------------
c
DO 70 I=1,NX
BC(I)=0.
BD(I)=0.
70 CONTINUE
c
DO 80 J=1,NY
BA(J)=0.
BB(J)=0.
80 CONTINUE
c
c --------------------------------------------------------
c Set the boundary data on the irregular boundary points
c --------------------------------------------------------
c
DO 90 N=1,NP
VP(N)=0.
90 CONTINUE
c
c ---------------------------------------------------------------
c Call CAPC with ICFLAG=1 to do the preprocessing and
c to obtain the solution in the irregular domain
c ---------------------------------------------------------------
c
ICFLAG=1
CALL CAPC(F,U,WORK,TRIGS,NX,NY,HELM,DX,DY,
# C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
# BA,BB,BC,BD,IERROR)
c
IF(IERROR.NE.0) THEN
WRITE(6,*)'*** ERROR CONDITION *** IERROR=',IERROR
STOP
END IF
c
c -------------------------------------------------------------------
c Apply the discrete Helmholtz operator to the computed solution
c over the interior of the solution domain and find the normalized
c L2 norm of the residual between the specified source function
c and the computed right hand side. Also find the largest residual
c and the number of correct digits, in an L2 sense, of the computed
c right hand side.
c -------------------------------------------------------------------
c
DINF=0.
DL2N=0.
DENM=0.
DO 110 J=2,NY1
DO 100 I=2,NX1
RR=ABS(REAL(I-1)*DX-0.5) + ABS(REAL(J-1)*DY-0.5)
IF(RR.LT.0.375) THEN
c
FD = (U(I+1,J)+U(I-1,J)-2.*U(I,J))/DX**2
# + (U(I,J+1)+U(I,J-1)-2.*U(I,J))/DY**2
# - HELM*U(I,J)
c
DL2N=DL2N+(FD-F(I,J))**2
DENM=DENM+F(I,J)**2
IF(ABS(FD-F(I,J)).GT.DINF) THEN
DINF=ABS(FD-F(I,J))
ID=I
JD=J
END IF
c
END IF
100 CONTINUE
110 CONTINUE
c
DL2N=SQRT(DL2N/DENM)
DIGT=-LOG10(DL2N)
WRITE(6,101) DIGT
WRITE(6,102) DL2N
WRITE(6,103) DINF, ID, JD
c
101 FORMAT(' NUMBER OF CORRECT DIGITS IN RESIDUAL =',F6.2)
102 FORMAT(' NORMALIZED L2 NORM OF RESIDUAL = ',1PE11.5)
103 FORMAT(' LARGEST ABSOLUTE RESIDUAL = ',1PE11.5,
# ' AT (I,J)= ',2(1X,I3))
STOP
END
c
c -----------------------------------------------------------
c Initialize the indices of the irregular boundary points:
c the coordinates of the curve |X-0.5| + |Y-0.5| = 0.375
c embedded in the unit square with discretization interval
c 1/16 in the X and Y directions.
c -----------------------------------------------------------
c
BLOCK DATA
PARAMETER(NP=24)
COMMON /BNDY/IP(NP),JP(NP)
DATA IP/9,10,11,12,13,14,15,14,13,12,11,10,9,8,7,
# 6,5,4,3,4,5,6,7,8/
DATA JP/3,4,5,6,7,8,9,10,11,12,13,14,15,14,13,12,11,
# 10,9,8,7,6,5,4/
END
C-END-OF-FILE
cat > testcapcn.f << C-END-OF-FILE
c -----------------------------------------------------------
c
c Program TESTCAPCN - Example of CAPCN usage
c
c -----------------------------------------------------------
c
c This program illustrates the use routine CAPCN, a solver for
c the nonseparable self-adjoint elliptic equation:
c (d/dX)(G dU/dX) + (d/dY)(G dU/dY) = Q(X,Y),
c in irregular (polygonal) two-dimensional domains with
c Dirichlet boundary conditions.
c
c The solution is obtain in the region defined by the equation
c |X-0.5| + |Y-0.5| <= 0.375, which is embedded in the unit
c square. The diffusion function is G(X,Y)=EXP(X+Y), and the
c source function is Q(X,Y)=COS(PI*X)*SIN(2*PI*Y). The
c boundary conditions are U(X,Y)=0 on |X-0.5| + |Y-0.5| = 0.375.
c The coordinates of the boundary are specified in the BLOCK DATA
c routine at the end.
c
c The solution produced by the solver is tested by applying the
c discrete elliptic operator to the computed solution and examining
c the difference (the residual) between this and the specified source
c function over the domain.
c
c --------------
c Declarations
c --------------
c
PARAMETER(NP=24,NX=17,NY=17,NX1=NX-1,NY1=NY-1,NTRG=2*(NX-1))
REAL F(NX,NY),U(NX,NY),WORK(NX,NY),TRIGS(NTRG)
REAL G(NX1,NY1),P(NX,NY)
REAL C(NP,NP),CHARGE(NP),VP(NP)
REAL BA(NY),BB(NY),BC(NX),BD(NX)
INTEGER IPS(NP),MASK(NX,NY)
COMMON /BNDY/ IP(NP),JP(NP)
c
c --------------------------------------
c Set the grid intervals and constants
c --------------------------------------
c
DX=1./REAL(NX1)
DY=1./REAL(NY1)
PI=4.*ATAN(1.)
IWRITE=6
c
c -------------------------------------------------------
c Set up the mask array to identify the points that are
c inside the irregular domain
c -------------------------------------------------------
c
c MASK(I,J) = 1 for interior points
c = 0 elsewhere, including all boundary points
c
DO 20 J=1,NY
DO 10 I=1,NX
XX=REAL(I-1)*DX
YY=REAL(J-1)*DY
RR=ABS(XX-0.5) + ABS(YY-0.5)
IF(RR.LT.0.375) THEN
MASK(I,J)=1
ELSE
MASK(I,J)=0
END IF
10 CONTINUE
20 CONTINUE
c
c -------------------------------
c Define the diffusion function
c -------------------------------
c
C1=1.
DO 50 J=1,NY1
DO 40 I=1,NX1
G(I,J)=EXP((REAL(I-1)+0.5)*DX+(REAL(J-1)+0.5)*DY)-C1
40 CONTINUE
50 CONTINUE
c
c -------------------------------------------
c Define the right-hand side source function
c -------------------------------------------
c
DO 70 J=2,NY1
DO 60 I=2,NX1
F(I,J)=COS(PI*REAL(I-1)*DX)*SIN(2.*PI*REAL(J-1)*DY)
60 CONTINUE
70 CONTINUE
c
c -----------------------------------------------------------
c Set the initial guess to the solution to zero. Also, for
c definiteness, set the arrays for boundary values
c on the sides of the embedding rectangle to zero.
c -----------------------------------------------------------
c
DO 90 J=1,NY
DO 80 I=1,NX
P(I,J)=0.
80 CONTINUE
90 CONTINUE
c
DO 100 I=1,NX
BC(I)=0.
BD(I)=0.
100 CONTINUE
c
DO 110 J=1,NY
BA(J)=0.
BB(J)=0.
110 CONTINUE
c
c --------------------------------------------------------
c Set the boundary data on the irregular boundary points
c --------------------------------------------------------
c
DO 120 N=1,NP
VP(N)=0.
120 CONTINUE
c
c --------------------------------------
c Set the maximum number of iterations
c --------------------------------------
c
ITMAX=25
c
c -------------------------------
c Set the convergence criterion
c -------------------------------
c
CCON=2.0E-4
c
c ------------------------------------------------------
c Call CAPCN with ICFLAG=1 to do the preprocessing and
c to obtain the solution in the irregular domain
c ------------------------------------------------------
c
ICFLAG=1
CALL CAPCN(F,G,U,P,WORK,MASK,TRIGS,NX,NY,DX,DY,
# C,CHARGE,IPS,VP,IP,JP,NP,ICFLAG,
# BA,BB,BC,BD,ITMAX,CCON,NITT,IWRITE,IERROR)
c
IF(IERROR.NE.0) THEN
WRITE(6,*)'*** ERROR CONDITION *** IERROR=',IERROR
STOP
END IF
c
WRITE(6,*) 'NUMBER OF ITERATIONS REQUIRED = ',NITT
c
c -------------------------------------------------------------------
c Apply the discrete Helmholtz operator to the computed solution
c over the interior of the solution domain and find the normalized
c L2 norm of the residual between the specified source function
c and the computed right hand side. Also find the largest residual
c and the number of correct digits, in an L2 sense, of the computed
c right hand side.
c -------------------------------------------------------------------
c
DINF=0.
DL2N=0.
DENM=0.
DO 140 J=2,NY1
DO 130 I=2,NX1
IF(MASK(I,J).EQ.1) THEN
c
FD=0.5*((G(I,J)+G(I,J-1))*(U(I+1,J)-U(I,J))
# - (G(I-1,J)+G(I-1,J-1))*(U(I,J)-U(I-1,J)))/DX**2
# + 0.5*((G(I,J)+G(I-1,J))*(U(I,J+1)-U(I,J))
# - (G(I,J-1)+G(I-1,J-1))*(U(I,J)-U(I,J-1)))/DY**2
DL2N=DL2N+(FD-F(I,J))**2
DENM=DENM+F(I,J)**2
IF(ABS(FD-F(I,J)).GT.DINF) THEN
DINF=ABS(FD-F(I,J))
ID=I
JD=J
END IF
c
END IF
130 CONTINUE
140 CONTINUE
c
DL2N=SQRT(DL2N/DENM)
DIGT=-LOG10(DL2N)
WRITE(6,101) DIGT
WRITE(6,102) DL2N
WRITE(6,103) DINF, ID, JD
c
101 FORMAT(' NUMBER OF CORRECT DIGITS IN RESIDUAL =',F6.2)
102 FORMAT(' NORMALIZED L2 NORM OF RESIDUAL = ',1PE11.5)
103 FORMAT(' LARGEST ABSOLUTE RESIDUAL = ',1PE11.5,
# ' AT (I,J)= ',2(1X,I3))
STOP
END
c
c -----------------------------------------------------------
c Initialize the indices of the irregular boundary points:
c the coordinates of the curve |X-0.5| + |Y-0.5| = 0.375
c embedded in the unit square with discretization interval
c 1/16 in the X and Y directions.
c -----------------------------------------------------------
c
BLOCK DATA
PARAMETER(NP=24)
COMMON /BNDY/IP(NP),JP(NP)
DATA IP/9,10,11,12,13,14,15,14,13,12,11,10,9,8,7,
# 6,5,4,3,4,5,6,7,8/
DATA JP/3,4,5,6,7,8,9,10,11,12,13,14,15,14,13,12,11,
# 10,9,8,7,6,5,4/
END
C-END-OF-FILE
cat > testdcapc.f <