C ALGORITHM 593, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 1, C MAR., 1983, P. 117-124. C MAN 10 C MAN 20 C DRIVER PROGRAM TO ILLUSTRATE THE USE OF CMMEXP MAN 30 C MAN 40 REAL C(1200), F(20,20), W(320), HX(40), HY(40), TX(40), TY(40), MAN 50 * BX(40), BY(40) MAN 60 INTEGER IW(400), IC(40), JC(40) MAN 70 LOGICAL DIR, LAPL, TEST MAN 80 C MAN 90 C NDIM IS THE FIRST DIMENSION OF THE ARRAY F. MAN 100 C THIS PARAMETER IS USED TO SPECIFY THE VARIABLE MAN 110 C DIMENSION OF F. NDIM MUST NOT BE LESS THAN N. MAN 120 C MAN 130 C LEW IS THE TOTAL DIMENSION OF THE WORKSPACE W. MAN 140 C LEW MUST NOT BE LESS THAN 8*IP MAN 150 C MAN 160 C LEIW IS THE TOTAL DIMENSION OF THE WORKSPACE IW. MAN 170 C LEIW MUST NOT BE LESS THAN 10*IP MAN 180 C MAN 190 C THE TOTAL DIMENSION OF THE ARRAY C MUST MAN 200 C NOT BE LESS THAN IP*IP. MAN 210 C MAN 220 C IPMAX IS THE DIMENSION OF VECTORS IC,JC,HX,HY,TX,TY,BX,BY MAN 230 C AND IS AN UPPER ESTIMATE OF IP, THE NUMBER OF IRREGULAR MAN 240 C MESH POINTS. MAN 250 C MAN 260 C MAN 270 C ******************** MAN 280 C * INITIALIZE * MAN 290 C ******************** MAN 300 C MAN 310 NDIM = 20 MAN 320 LEW = 320 MAN 330 LEIW = 400 MAN 340 IPMAX = 40 MAN 350 CON = 0.0 MAN 360 DIR = .TRUE. MAN 370 LAPL = .FALSE. MAN 380 TEST = .TRUE. MAN 390 N = 16 MAN 400 M = N MAN 410 INTL = 0 MAN 420 IOUT = 6 MAN 430 C MAN 440 C THE TEST IS PERFORMED IN THE CIRCULAR REGION MAN 450 C MAN 460 C ( X - 0.5 ) ** 2 + ( Y - 0.5 ) ** 2 = R ** 2, R=0.375, MAN 470 C MAN 480 C WITH THE EXACT SOLUTION BEING U(X,Y) = X**5 + Y**5. MAN 490 C MAN 500 C ******************** MAN 510 C * INPUT DATA * MAN 520 C ******************** MAN 530 C MAN 540 C PROCESS THE INFORMATION ABOUT THE REGION AND ITS BOUNDARY MAN 550 C MAN 560 C INPUT DATA ABOUT THE REGION: DEFINE IP,IPP,H, MAN 570 C IC,JC,HX,HY AND TX,TY(FOR THE NEUMANN PROBLEM). MAN 580 C MAN 590 CALL DOMAIN(N, M, IC, JC, HX, HY, TX, TY, IPMAX, IP, IPP, H, DIR) MAN 600 C MAN 610 C ******************** MAN 620 C * CALL CMMEXP * MAN 630 C ******************** MAN 640 C MAN 650 C PREPROCESSING STAGE MAN 660 C MAN 670 CALL CMMEXP(INTL, IC, JC, HX, HY, TX, TY, BX, BY, F, NDIM, N, M, MAN 680 * H, CON, DIR, LAPL, TEST, C, IPMAX, IP, IPP, W, LEW, IW, LEIW, MAN 690 * IOUT, IERROR) MAN 700 C MAN 710 IF (IERROR.GT.0) WRITE (IOUT,99999) IERROR MAN 720 IF (IERROR.GT.0) STOP MAN 730 C MAN 740 C ******************** MAN 750 C * INPUT DATA * MAN 760 C ******************** MAN 770 C MAN 780 C INPUT BOUNDARY DATA: DEFINE BX AND BY. MAN 790 C MAN 800 CALL BNDRY(IC, JC, HX, HY, BX, BY, IP, IPP, H, DIR) MAN 810 C MAN 820 C INPUT RIGHT SIDE: DEFINE F MAN 830 C MAN 840 CALL CHARGE(F, NDIM, N, M, H, CON) MAN 850 C MAN 860 C ******************** MAN 870 C * CALL CMMEXP * MAN 880 C ******************** MAN 890 C MAN 900 C SOLVING STAGE MAN 910 C MAN 920 INTL = INTL + 1 MAN 930 C MAN 940 CALL CMMEXP(INTL, IC, JC, HX, HY, TX, TY, BX, BY, F, NDIM, N, M, MAN 950 * H, CON, DIR, LAPL, TEST, C, IPMAX, IP, IPP, W, LEW, IW, LEIW, MAN 960 * IOUT, IERROR) MAN 970 C MAN 980 C ******************** MAN 990 C * PRINTOUT * MAN 1000 C ******************** MAN 1010 C MAN 1020 C FOR THE TEST EXAMPLE MAN 1030 C PRINT OUT THE NUMBER OF CORRECT DIGITS MAN 1040 C FOR ALL POINTS INSIDE THE IRREGULAR REGION. MAN 1050 C MAN 1060 IF (TEST .AND. IERROR.EQ.0) CALL OUTPUT(F, NDIM, N, M, H, DIR, MAN 1070 * IOUT) MAN 1080 C MAN 1090 C NOTE: SUBROUTINES DOMAIN, BNDRY, CHARGE AND OUTPUT ARE USED MAN 1100 C ONLY IN THIS DRIVER AND ARE NOT REQUIRED IN CMMEXP. MAN 1110 C MAN 1120 99999 FORMAT (//1X, 19H DETECTED ERROR NO., I3//) MAN 1130 STOP MAN 1140 END MAN 1150 C MAN 10 C DRIVER PROGRAM TO ILLUSTRATE THE USE OF CMMSIX MAN 20 C MAN 30 REAL C(5000), F(40,40), W(3500), TX(1), TY(1), HX(80), HY(80), MAN 40 * BX(80), BY(80) MAN 50 DIMENSION IC(80), JC(80), IW(200) MAN 60 LOGICAL TEST MAN 70 C MAN 80 C NDIM IS THE FIRST DIMENSION OF THE ARRAY F. MAN 90 C THIS PARAMETER IS USED TO SPECIFY THE VARIABLE MAN 100 C DIMENSION OF F. NDIM MUST NOT BE LESS THAN N. MAN 110 C MAN 120 C LEW IS THE TOTAL DIMENSION OF THE WORKSPACE W. MAN 130 C LEW MUST NOT BE LESS THAN 2*N*M + (8+4*NCOR)*IP MAN 140 C MAN 150 C LEIW IS THE TOTAL DIMENSION OF THE WORKSPACE IW. MAN 160 C LEIW MUST NOT BE LESS THAN 2*IP MAN 170 C MAN 180 C THE TOTAL DIMENSION OF THE ARRAY C MUST MAN 190 C NOT BE LESS THAN IP*IP. MAN 200 C MAN 210 C IPMAX IS THE DIMENSION OF VECTORS IC,JC,HX,HY,BX,BY MAN 220 C AND IS AN UPPER ESTIMATE OF IP, THE NUMBER OF IRREGULAR MAN 230 C MESH POINTS. MAN 240 C MAN 250 C MAN 260 C ******************** MAN 270 C * INITIALIZE * MAN 280 C ******************** MAN 290 C MAN 300 NDIM = 40 MAN 310 LEW = 3500 MAN 320 LEIW = 200 MAN 330 IPMAX = 80 MAN 340 TEST = .TRUE. MAN 350 N = 32 MAN 360 M = N MAN 370 IOUT = 6 MAN 380 DO 10 NC=1,3 MAN 390 NCOR = NC - 1 MAN 400 INTL = 0 MAN 410 C MAN 420 C THE TEST IS PERFORMED IN THE CIRCULAR REGION MAN 430 C MAN 440 C ( X - 0.5 ) ** 2 + ( Y - 0.5 ) ** 2 = R ** 2, R=0.375, MAN 450 C MAN 460 C WITH THE EXACT SOLUTION BEING U(X,Y) = X**5 + Y**5. MAN 470 C MAN 480 C ******************** MAN 490 C * INPUT DATA * MAN 500 C ******************** MAN 510 C MAN 520 C PROCESS THE INFORMATION ABOUT THE REGION AND ITS BOUNDARY MAN 530 C MAN 540 C INPUT DATA ABOUT THE REGION: DEFINE IP,H,IC,JC,HX,HY. MAN 550 C MAN 560 CALL DOMAIN(N, M, IC, JC, HX, HY, TX, TY, IPMAX, IP, IPP, H, MAN 570 * .TRUE.) MAN 580 C MAN 590 C ******************** MAN 600 C * CALL CMMSIX * MAN 610 C ******************** MAN 620 C MAN 630 C PREPROCESSING STAGE MAN 640 C MAN 650 CALL CMMSIX(INTL, IC, JC, HX, HY, BX, BY, F, NDIM, N, M, H, MAN 660 * TEST, C, IPMAX, IP, W, LEW, IW, LEIW, NCOR, IOUT, IERROR) MAN 670 C MAN 680 IF (IERROR.GT.0) WRITE (IOUT,99999) IERROR MAN 690 IF (IERROR.GT.0) STOP MAN 700 C MAN 710 C ******************** MAN 720 C * INPUT DATA * MAN 730 C ******************** MAN 740 C MAN 750 C INPUT BOUNDARY DATA: DEFINE BX AND BY. MAN 760 C MAN 770 CALL BNDRY(IC, JC, HX, HY, BX, BY, IP, IPP, H, .TRUE.) MAN 780 C MAN 790 C INPUT RIGHT SIDE: DEFINE F MAN 800 C MAN 810 CALL CHARGE(F, NDIM, N, M, H, 0.0) MAN 820 C MAN 830 C ******************** MAN 840 C * CALL CMMSIX * MAN 850 C ******************** MAN 860 C MAN 870 C SOLVING STAGE MAN 880 C MAN 890 INTL = INTL + 1 MAN 900 C MAN 910 CALL CMMSIX(INTL, IC, JC, HX, HY, BX, BY, F, NDIM, N, M, H, MAN 920 * TEST, C, IPMAX, IP, W, LEW, IW, LEIW, NCOR, IOUT, IERROR) MAN 930 C MAN 940 C ******************** MAN 950 C * PRINTOUT * MAN 960 C ******************** MAN 970 C MAN 980 C FOR THE TEST EXAMPLE MAN 990 C PRINT OUT THE NUMBER OF CORRECT DIGITS MAN 1000 C FOR ALL POINTS INSIDE THE IRREGULAR REGION. MAN 1010 C MAN 1020 IF (TEST .AND. IERROR.EQ.0) CALL OUTPUT(F, NDIM, N, M, H, MAN 1030 * .TRUE., IOUT) MAN 1040 C MAN 1050 C NOTE: SUBROUTINES DOMAIN, BNDRY, CHARGE AND OUTPUT ARE USED MAN 1060 C ONLY IN THIS DRIVER AND ARE NOT REQUIRED IN CMMSIX. MAN 1070 C MAN 1080 10 CONTINUE MAN 1090 99999 FORMAT (//1X, 19H DETECTED ERROR NO., I3//) MAN 1100 STOP MAN 1110 END MAN 1120 SUBROUTINE CMMIMP(IC, JC, HX, HY, TX, TY, BX, BY, F, NDIM, N, M, CMM 10 * H, CON, DIR, LAPL, TEST, ACCY, ITMAX, IPMAX, IP, IPP, W, LEW, * IW, LEIW, IOUT, IERROR) REAL F(NDIM,M), W(LEW), HX(IPP), HY(IPP), TX(IPP), TY(IPP), * BX(IPP), BY(IPP) INTEGER IW(LEIW), IC(IP), JC(IP) LOGICAL DIR, LAPL, TEST C C THIS PROGRAM SOLVES THE DIRICHLET AND NEUMANN PROBLEMS C FOR THE HELMHOLTZ EQUATION OVER AN ARBITRARY BOUNDED PLANAR REGION C IMBEDDED IN A RECTANGLE A,B;C,D : C C ( -L + C )U = F IN THE REGION, C C U = G C OR ON THE BOUNDARY, C DU/DN = G C C WHERE L IS THE LAPLACIAN, DU/DN IS THE NORMAL DERIVATIVE, C IS C A REAL CONSTANT, AND F AND G ARE SPECIFIED FUNCTIONS OF X AND Y. C C THE ALGORITHM IS AN ITERATIVE VARIANT OF THE CAPACITANCE C MATRIX METHOD. FOR THE DESCRIPTION OF THE METHOD SEE C W.PROSKUROWSKI - ACM TRANS.MATH.SOFTWARE, VOL.5, 1979, P.36-49. C C THE SUBROUTINE CMMIMP IS DESIGNED FOR USE WITH LARGE MESHES C ( APPROXIMATELY UP TO N*M = 500*500 AND IP = 1500 ASSUMING C THE CORE STORAGE OF 65K WORDS) IF ONE OR VERY FEW PROBLEMS C ARE SOLVED FOR A GIVEN REGION. C FOR PROBLEMS WITH MULTIPLE RIGHT HAND SIDES F ONE SHOULD USE C THE SUBROUTINE CMMEXP INSTEAD, IF THE MEMORY SPACE REQUIRED C DOES NOT EXCEED THE FAST CORE MEMORY. C C THE EMBEDDING RECTANGLE IS COVERED WITH A SQUARE MESH OF LENGTH H. C IT SHOULD BE CHOSEN TO FOLLOW CLOSELY THE SHAPE OF THE REGION C IN SUCH A WAY THAT AT LEAST ONE MESH LENGTH IS LEFT C BETWEEN THE BOUNDARIES OF THE REGION AND THE RECTANGLE. C C FOR MEMORY SPACE REQUIREMENT SEE PARAMETERS F, LEW, LEIW AND IPMAX. C C THIS PROGRAM WAS DEVELOPED BY W.PROSKUROWSKI. C THIS IS AN AUGUST 1981 VERSION. C C THE ARGUMENTS ARE DEFINED AS : C C C * * * * * * * * * * ON INPUT * * * * * * * * * * * * * * * * C C C IC,JC C C INTEGER VECTORS OF LENGTH IP CONTAINING COORDINATES (I,J) C OF ALL IRREGULAR POINTS ENUMERATED IN AN ARBITRARY ORDER. C THE VALUES OF I MUST BE LARGER THAN 1 AND LESS THAN N. C THE VALUES OF J MUST BE LARGER THAN 1 AND LESS THAN M. C A CALL OF CMMIMP DOES NOT DESTROY THE VALUES OF IC AND JC. C C HX,HY C C REAL VECTORS OF LENGTH IPP WITH SIGNED SHORTEST C DISTANCES IN BOTH COORDINATE DIRECTIONS FROM THE BOUNDARY C TO AN IRREGULAR POINT (I,J) IN UNITS OF STEP SIZE H. C IF (I,J) HAS A NEIGHBOR OUTSIDE THE REGION THEN THE VALUES OF C HX OR HY MUST BELONG TO THE INTERVAL -1.0 TO +1.0 (EXCLUDING 0.0 C OTHERWISE THE PROPER VALUE OF MAGNITUDE LARGER THAN 1.0 C MUST BE PUT IN. C A CALL OF CMMIMP DOES NOT DESTROY THE VALUES OF HX AND HY. C C TX,TY C C REAL VECTORS OF LENGTH IPP WITH ABSOLUTE VALUES OF C TANGENTS OF THE ANGLE BETWEEN THE NORMAL AND THE X-AXIS C AT (I+HX,J) ( AND THE Y-AXIS AT (I,J+HY) ). C TX,TY ARE USED ONLY FOR A PROBLEM WITH THE NEUMANN BOUNDARY C CONDITIONS. THE VALUES OF TX,TY MUST BE POSITIVE. C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0, A DUMMY C VALUE OF ARBITRARY MAGNITUDE CAN BE PUT IN. C A CALL OF CMMIMP DOES NOT DESTROY THE VALUES OF TX AND TY. C C BX,BY C C REAL VECTORS OF LENGTH IPP WITH BOUNDARY VALUES C AT (IC(I)+HX,JC(I)) AND (IC(I),JC(I)+HY). C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0 THEN C A DUMMY VALUE OF BX OR BY SHOULD BE PUT IN. C A CALL OF CMMIMP DESTROYS THE VALUES OF BX AND BY. C C IP C C THE NUMBER OF IRREGULAR MESH POINTS, I.E., POINTS C STRICTLY INSIDE THE REGION WHICH HAVE AT LEAST ONE C NEIGHBOR OUTSIDE OR ON THE BOUNDARY. C C IPP C C LENGTH OF VECTORS HX,HY,TX,TY,BX,BY. C C THE SET OF IRREGULAR POINTS CONSISTS OF THE SET (1) OF IP1 C POINTS WITH AT MOST 2 NEIGHBORS OUTSIDE THE REGION (OR ON C THE BOUNDARY) AND THE REMAINING SET (2) OF IP2 POINTS WITH C EXACTLY 3 NEIGHBORS OUTSIDE THE REGION (OR ON THE BOUNDARY). C THUS IP = IP1 + IP2. C IF A POINT BELONGS TO THE SET (2) THEN THE DISTANCES FROM C THE BOUNDARY MUST BE STORED IN A PAIR OF LOCATIONS: C (HX(I),HX(I+1)) AND (HY(I),HY(I+1)). NOTE THAT THE SECOND ENTRY C MUST BE LARGER IN MAGNITUDE THAN THE FIRST. A SIMILAR SITUATION C OCCURS FOR THE BOUNDARY VALUES BX,BY AND THE TANGENTS TX,TY. C THUS THE LENGTH OF THESE VECTORS MUST BE IPP = IP1 + IP2*2. C C IPMAX C C IS AN A PRIORI ESTIMATE OF IP AND THE DIMENSION OF VECTORS C IC,JC,HX,HY,TX,TY,BX,BY AS IT APPEARS IN THE PROGRAM CALLING C CMMIMP. IPMAX MUS BE AT LEAST IPP. C C F C C REAL ARRAY OF ORDER NDIM BY M CONTAINING THE VALUES OF THE C RIGHT HAND SIDE OF THE HELMHOLTZ EQUATION FOR THE GRID POINTS C INSIDE THE IRREGULAR REGION SUPPLEMENTED BY SOME ARBITRARY C VALUES, FOR INSTANCE BY ZERO VALUES, FOR THE GRID POINTS C OUTSIDE THE IRREGULAR REGION. C F DOES NOT NEED TO BE INITIALIZED IF LAPL=.TRUE. C C NDIM C C THE FIRST DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING CMMIMP. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. NDIM MUST BE AT LEAST N. C C N C C THE NUMBER OF MESH LENGTHS IN THE X DIRECTION. C N MUST BE A POWER OF 2 AND NOT LESS THAN 8. C C M C C THE NUMBER OF MESH LENGTHS IN THE Y DIRECTION. C M SHOULD APPROXIMATELY BE M = N * (D-C)/(B-A), C WHERE A,B IS THE RANGE OF X IN THE RECTANGULAR REGION, C I.E., A =< X =< B, AND C,D IS THE RANGE OF Y C IN THE RECTANGULAR REGION, I.E., C =< Y =< D. C M MUST BE GREATER THAN OR EQUAL TO N. C C H C C THE MESH SIZE IN BOTH COORDINATE DIRECTIONS. C THE VALUE USED IN THE PACKAGE IS H=(B-A)/N. C C CON C C REAL VARIABLE, THE CONSTANT IN THE HELMHOLTZ EQUATION. C C DIR C C LOGICAL VARIABLE, EQUAL TO .TRUE. FOR A DIRICHLET PROBLEM, C EQUAL TO .FALSE. FOR A NEUMANN PROBLEM. C NOTE THAT THE NEUMANN PROBLEM FOR THE LAPLACE EQUATION, I.E., C WITH DIR=.FALSE. AND CON=0.0, IS SINGULAR. IN SUCH CASES THE C LEAST SQUARE SOLUTION IS COMPUTED AND IT MAY NEED NORMALIZATION C BY SUBTRACTING A PROPER CONSTANT. C C LAPL C C LOGICAL VARIABLE, EQUAL TO .TRUE. IF THE LAPLACE EQUATION C IS TO BE SOLVED, .FALSE. OTHERWISE, I.E., IF THE RIGHT SIDE C F IS NONZERO. A CALL WITH LAPL=.TRUE. IS LESS EXPENSIVE C BY APPROXIMATELY 5 TO 10 %. C C TEST C C LOGICAL VARIABLE, EQUAL TO .TRUE. FOR TESTING OR DEBUGGING C PURPOSES, .FALSE. OTHERWISE,I.E., IN THE PRODUCTION MODE. C C ACCY C C REAL VARIABLE, THE REQUIRED TOLERANCE IN THE CONJUGATE GRADIENT C ITERATION. THE ITERATION IS TERMINATED IF THE L2-NORM C OF THE RESIDUAL DROPS BELOW ACCY. THE RECOMMENDED VALUE IS 1E-3. C CORRELATION BETWEEN ACCY AND THE SOLUTION ERROR DEPENDS C HEAVILY ON THE DISCRETIZATION (TRUNCATION) ERROR IN THE PROBLEM. C C ITMAX C C INTEGER VARIABLE, THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. C THE RECOMMENDED VALUE IS 20. C C W C C A REAL ARRAY THAT MUST BE PROVIDED FOR WORKSPACE. C C LEW C C THE TOTAL DIMENSION OF THE WORKSPACE W. C FOR THE DIRICHLET PROBLEM C LEW MUST BE NOT LESS THAN 6*IP C FOR THE NEUMANN PROBLEM C LEW MUST BE NOT LESS THAN 7*IP C C IW C C AN INTEGER ARRAY THAT MUST BE PROVIDED FOR WORKSPACE. C C LEIW C C THE TOTAL DIMENSION OF THE WORKSPACE IW. C FOR THE DIRICHLET PROBLEM C LEIW MUST BE NOT LESS THAN 19*IP C FOR THE NEUMANN PROBLEM C LEIW MUST BE NOT LESS THAN 13*IP C C IOUT C C INTEGER VARIABLE, THE NUMBER OF THE OUTPUT DEVICE. C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * * * * * * * * * C C F C C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINTS INSIDE THE IRREGULAR REGION, C AND LARGELY ARBITRARY VALUES FOR POINTS OUTSIDE THE REGION. C C IERROR C C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. C IF IERROR IS GREATER THAN 0 AND LESS THAN 7 NO SOLUTION IS SOUGH C C = 0: NO ERROR C = 1: N IS NOT A POWER OF 2 OR LESS THAN 16, OR M IS LESS THAN N C = 2: NDIM (THE FIRST DIMENSION OF ARRAY F) IS LESS THAN N. C = 3: IPMAX (THE ESTIMATE OF IP) IS TOO SMALL. C = 4: LEW OR LEIW (THE DIMENSIONS OF WORKSPACE) ARE TOO SMALL. C = 5: REGION TOO CLOSE TO THE EMBEDDING RECTANGLE. C = 6: INTERNAL ERROR IN THE PACKAGE. C = 7: NUMBER OF ITERATIONS EXCEEDED THE PRESCRIBED MAXIMUM C OR ITERATION DIVERGES ( BAD DATA ). C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO CMMIMP, THE USER SHOULD TEST IERROR AFTER A CALL. C IF IERROR > 0, ONE CAN RUN CMMIMP AGAIN WITH TEST=.TRUE. C (DEBUGGING MODE) TO OBTAIN SOME ADDITIONAL INFORMATION. C C ITMAX C C CONTAINS THE VALUE 999 IF THE PRESCRIBED MAXIMUM NUMBER OF C ITERATIONS HAS BEEN EXCEEDED OR ITERATION DIVERGES ( BAD DATA ). C C C ********************************************************************** C C IERROR = 0 NN = 4 10 NN = NN*2 IF (N.GT.NN) GO TO 10 IF (N.LT.8 .OR. N.NE.NN .OR. M.LT.N) IERROR = 1 IF (IERROR.EQ.1) RETURN IF (TEST) WRITE (IOUT,99994) N, M, NDIM IF (N.LE.NDIM) GO TO 20 IERROR = 2 RETURN 20 CONTINUE IF (IPP.GT.IPMAX) IERROR = 3 IF (TEST) WRITE (IOUT,99999) IP, IPP, IPMAX IP1 = IP - (IPP-IP) IF (DIR) LEN = MIN0(LEW/6,LEIW/19) IF (.NOT.DIR) LEN = MIN0(LEW/7,LEIW/13) IF (LEN.LT.IPP) IERROR = 4 FACTOR = FLOAT(IPP)/FLOAT(LEN) IF (TEST .AND. IERROR.EQ.4) WRITE (IOUT,99995) FACTOR IF (IERROR.GT.0) RETURN DO 30 I=1,IP K = IC(I) L = JC(I) IF (K.LE.1 .OR. K.GE.N .OR. L.LE.1 .OR. L.GE.N) IERROR = 5 30 CONTINUE C IF (.NOT.TEST) GO TO 50 C C CHECK THE INPUT VALUES BY PRINTING THEM OUT C WRITE (IOUT,99993) IF (DIR) WRITE (IOUT,99992) IF (.NOT.DIR) WRITE (IOUT,99991) IP1 = IP - (IPP-IP) DO 40 I=1,IP J = I IF (J.GT.IP1) J = IP1 + (I-IP1)*2 - 1 IF (DIR) WRITE (IOUT,99990) I, IC(I), JC(I), HX(J), HY(J) IF (DIR .AND. I.GT.IP1) WRITE (IOUT,99989) HX(J+1), HY(J+1) IF (.NOT.DIR) WRITE (IOUT,99990) I, IC(I), JC(I), HX(J), HY(J), * TX(J), TY(J) IF (.NOT.DIR .AND. I.GT.IP1) WRITE (IOUT,99989) HX(J+1), * HY(J+1), TX(J+1), TY(J+1) IF (IC(I).LE.1 .OR. IC(I).GE.N .OR. JC(I).LE.1 .OR. JC(I).GE.N) * WRITE (IOUT,99988) I, IC(I), JC(I) 40 CONTINUE 50 CONTINUE IF (IERROR.GT.0) RETURN C C GUESS THE APPROXIMATE VALUES OF NR,IL,NO AND IN. C NR = IPMAX*2 NROLD = NR IF (DIR) Z = 2.5 IF (.NOT.DIR) Z = 1. IL = FLOAT(IPMAX)*Z ILOLD = IL NO = IPMAX*3 NOOLD = NO I2 = NR + 1 I3 = I2 + NR I4 = I3 + IL I5 = I4 + IL I6 = I5 + NO IN = NR INOLD = IN IPP1 = IPP + 1 C DO 60 I=1,IP IW(I) = IC(I) NRI = NR + I IW(NRI) = JC(I) 60 CONTINUE C C PREPROCESS THE BOUNDARY INFORMATION INTO A COMPACT FORM C AND COMPUTE THE VALUES OF NR,IL,NO AND IN. C CALL COMPRS(IW(1), IW(I2), NR, IW(I3), IW(I4), IL, IW(I5), NO, * IW(I6), IN, HX, HY, IPP, IP, DIR) C IF (NR.GT.NROLD .OR. IL.GT.ILOLD .OR. NO.GT.NOOLD .OR. * IN.GT.INOLD) IERROR = 6 IF (TEST .AND. IERROR.EQ.6) WRITE (IOUT,99996) NR, NROLD, IL, * ILOLD, NO, NOOLD, IN, INOLD IF (TEST) WRITE (IOUT,99997) NR, IL, NO, IN IF (IERROR.EQ.6) RETURN C C COMPRESS THE VECTOR IW. C L = IPMAX*2 DO 70 I=1,NR INR = I + NR LI = I + L IW(INR) = IW(LI) 70 CONTINUE K = NR*2 L = IPMAX*4 DO 90 J=1,2 DO 80 I=1,IL LI = I + L KI = I + K IW(KI) = IW(LI) 80 CONTINUE K = K + IL IZ = FLOAT(IPMAX)*Z L = L + IZ 90 CONTINUE DO 100 I=1,NO LI = I + L KI = I + K IW(KI) = IW(LI) 100 CONTINUE K = K + NO L = L + IPMAX*3 DO 110 I=1,IN LI = I + L KI = I + K IW(KI) = IW(LI) 110 CONTINUE MX = MAX0(NR,IL) I2 = NR + 1 I3 = I2 + NR I4 = I3 + IL I5 = I4 + IL I6 = I5 + NO IADD = IN IF (.NOT.DIR) IADD = 0 I7 = I6 + IADD I8 = I7 + MX I9 = I8 + MX LI = I9 + N - 1 J2 = 1 + IP J3 = J2 + MX J4 = J3 + MX J5 = J4 + M J6 = J5 + M J7 = J6 + N J8 = J7 + N LW = J8 + N/4 - 1 IF (TEST) WRITE (IOUT,99998) LW, LI, LEW, LEIW C C CALL THE MAIN SUBROUTINE C CALL HLMHLZ(F, NDIM, N, M, H, CON, DIR, LAPL, TEST, ACCY, ITMAX, * IW(1), IW(I2), IW(I3), IW(I4), IW(I5), IW(I6), IW(I7), IW(I8), * IW(I9), HX, HY, TX, TY, BX, BY, W(1), W(J2), W(J3), W(J4), * W(J5), W(J6), W(J7), W(J8), IPP, IP, NR, IL, NO, IN, IOUT) C IF (ITMAX.EQ.999) IERROR = 7 C RETURN 99999 FORMAT (/1X, 4H IP=, I5, 5H IPP=, I5, 7H IPMAX=, I5/) 99998 FORMAT (/48H THE ACTUAL NEED FOR THE WORKSPACE W AND IW WAS, * 2I5/48H VS THE ESTIMATED SPACE RESERVED AS LEW AND LEIW, 2I5/) 99997 FORMAT (/1X, 35H THE LENGTH OF VECTORS IC,JC IS NR=, I5/1X, * 35H THE LENGTH OF VECTORS KC,LC IS IL=, I5/1X, 14H THE LENGTH OF, * 21H VECTOR INR IS NO=, I5/1X, 28H THE LENGTH OF VECTOR IDL I, * 7HS IN=, I5) 99996 FORMAT (1X, 50H WRONG ESTIMATE OF THE INTERMEDIATE VALUES NR,IL,N, * 4HO,IN/43H CONSULT W.PROSKUROWSKI AT USC, LOS ANGELES) 99995 FORMAT (1X, 50H INSUFFICIENT SPACE RESERVED FOR THE WORKSPACE W O, * 4HR IW/49H INCREASE THE DIMENSIONS LEW AND LEIW BY A FACTOR, * 12H OF AT LEAST, F5.2) 99994 FORMAT (1X, 3H N=, I3, 4H M=, I3, 11H AND NDIM=, I5) 99993 FORMAT (/1X, 24H THIS IS A TEST PRINTOUT/1X, 17H THE COORDINATES , * 24HOF THE IRREGULAR POINTS,/1X, 30H THE DISTANCES TO THE BOUNDARY * /1X, 40H AND TANGENTS (FOR THE NEUMANN PROBLEM) /1X, 9H FOR ALL , * 32HIRREGULAR POINTS ARE PRINTED OUT/) 99992 FORMAT (/52H NO I J HX HY /) 99991 FORMAT (/52H NO I J HX HY TX TY/) 99990 FORMAT (3I5, 4F10.5) 99989 FORMAT (15X, 4F10.5) 99988 FORMAT (/3I5, 21H BOUNDARIES TOO CLOSE) END SUBROUTINE HLMHLZ(F, NDIM, N, M, H, CON, DIR, LAPL, TEST, ACCY, HLM 10 * ITMAX, IC, JC, KC, LC, INR, IDL, IIS, IIT, IB, HXX, HYY, TNX, * TNY, P, RHS, SOL, AP, W, RX, AIX, SAJN, COSAJN, SI, IPP, IP, NR, * IL, NO, IN, IOUT) DIMENSION F(NDIM,M) LOGICAL DIR, LAPL, TEST DIMENSION IC(NR), JC(NR), KC(IL), LC(IL), INR(NO), IDL(IN) DIMENSION HXX(IPP), HYY(IPP), SOL(IP), RHS(IP), P(IP), AP(IL), * W(IL) DIMENSION TNX(IPP), TNY(IPP), RX(M), AIX(M), IIS(IL), IIT(IL) DIMENSION SAJN(N), COSAJN(N), IB(N), SI(1) C C C PARAMETERS C ARE DESCRIBED BELOW SO FAR AS THEY ARE NOT INCLUDED C OR DIFFER FROM THOSE OF THE SUBROUTINE CNNIMP. C C NR C LENGTH OF VECTORS IC,JC. APPROXIMATELY EQUAL TO 2*IP. C C IL C LENGTH OF VECTORS KC,LC,W,AP. APPROXIMATELY EQUAL TO 2.5*IP. C C NO C LENGTH OF VECTOR INR. APPROXIMATELY EQUAL TO 3*IP. C C IN C LENGTH OF VECTOR IDL. APPROXIMATELY EQUAL TO 2*IP. C C IC,JC C INTEGER VECTORS OF ORDER NR. CONTAIN COORDINATES OF THE C IRREGULAR POINTS AND THEIR NEIGHBORS INSIDE THE REGION. C C KC,LC C INTEGER VECTORS OF ORDER IL. CONTAIN COORDINATES OF THE C DISCRETE DIPOLES. C C INR,IDL C INTEGER VECTORS, CONTAIN CERTAIN INFORMATION ABOUT THE C COORDINATES OF POINTS IN IC,JC,KC,LC PREPROCESSED IN COMPRS. C C SOL,RHS,W C REAL VECTORS OF LENGTH IP,IP & IL. C SERVE AS A WORKSPACE IN CONJUGATE GRADIENT ITERATIONS C ( TOGETHER WITH VECTORS P AND AP ). C C IIS,IIT C INTEGER VECTORS OF LENGTH MAX(NR,IL). C USED FOR THE SPARSE SOLVER. C C RX,AIX C REAL VECTORS OF LENGTH M. C USED FOR THE SPARSE SOLVER. C C SAJN,COSAJN C REAL VECTORS OF LENGTH N. C USED FOR THE SPARSE SOLVER. C C IB C INTEGER VECTOR OF LENGTH N. C USED FOR THE FFT SOLVER. C C SI C REAL VECTOR OF LENGTH N/4. C USED FOR THE FFT SOLVER. C C C APPROXIMATE SPACE REQUIREMENTS FOR VECTORS C IN UNITS OF LENGTH EQUAL TO IP. C C DIRICHLET NEUMANN C IC,JC 2 2 C KC,LC 2.5 1 C INR 3 3 C IDL 2 0 C HXX,HYY,SOL,RHS,P 1 1 C W,AP 2.5 2 C TNX,TNY 0 1 C C ADDITIONAL SPACE USED IN THE SUBROUTINE STRIP C C IIT,IIS 2.5 2 C RX,AIX ABOUT 0.5 ABOUT 0.5 C SAJN,COSAJN,IB ABOUT 0.5 ABOUT 0.5 C C TOTAL ABOUT 31 ABOUT 26 C C ********************* C LIST OF SUBROUTINES : C ********************* C C COMPRS - PROCESSING THE COORDINATES OF THE IRREGULAR POINTS C INTO A COMPACT FORM. C C RIGHT - COMPUTES THE RIGHT HAND SIDE OF THE CAPACITANCE SYSTEM. C C STEN5 - MAPPING THE BOUNDARY STENCIL INTO THE IRREGULAR POINTS. C C NETS5 - A TRANSPOSE OF STEN5. C C DIPOLE - MAPPING THE DIPOLES INTO THE IRREGULAR POINTS. C C ELOPID - A TRANSPOSE OF DIPOLE. C C STRIP - FAST HELMHOLTZ SOLVER. C C RFORT,FORT - FAST FOURIER TRANSFORMATION ROUTINES. C HSQ = H*H CH2 = CON*HSQ IP1 = IP - (IPP-IP) R1 = 1.0 C C T T T C COMPUTE RHS = V G(U A) SOL C DO 20 J=1,M DO 10 I=1,N F(I,J) = F(I,J)*HSQ 10 CONTINUE 20 CONTINUE C CALL RIGHT(F, NDIM, N, M, P, RHS, SOL, IC, JC, IPP, IP, HXX, HYY, * TNX, TNY, H, CH2, DIR, LAPL) IF (LAPL) GO TO 40 C C CALCULATE THE SPACE POTENTIAL. C CALL STRIP(P, KC, LC, IP, W, IC, JC, NR, .TRUE., F, NDIM, N, M, * CH2, 2, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) CALL STEN5(W, NR, RHS, IPP, IP, INR, HXX, HYY, TNX, TNY, CH2, DIR) DO 30 I=1,IP SOL(I) = SOL(I) - RHS(I) 30 CONTINUE 40 CONTINUE CALL NETS5(SOL, IPP, IP, W, NR, INR, HXX, HYY, TNX, TNY, CH2, DIR) CALL STRIP(W, IC, JC, NR, AP, KC, LC, IL, LAPL, F, NDIM, N, M, * CH2, 1, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) CALL ELOPID(AP, IL, RHS, IPP, IP, IDL, HXX, HYY, DIR) C C CONJUGATE GRADIENT ITERATION C FIP = FLOAT(IP) RR = 0. DO 50 L=1,IP RR = RR + RHS(L)*RHS(L) P(L) = RHS(L) SOL(L) = 0. 50 CONTINUE C C ************************ C * MAIN ITERATION LOOP * C ************************ C IF (TEST) WRITE (IOUT,99998) C DO 90 KIT=1,ITMAX C T T T T C AP = V .G.(U A) .(U A).G.V.P C CALL DIPOLE(P, IPP, IP, W, IL, IDL, HXX, HYY, DIR) CALL STRIP(W, KC, LC, IL, AP, IC, JC, NR, .FALSE., F, NDIM, N, * M, CH2, 1, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) CALL STEN5(AP, NR, W, IPP, IP, INR, HXX, HYY, TNX, TNY, CH2, * DIR) CALL NETS5(W, IPP, IP, AP, NR, INR, HXX, HYY, TNX, TNY, CH2, * DIR) CALL STRIP(AP, IC, JC, NR, W, KC, LC, IL, .FALSE., F, NDIM, N, * M, CH2, 1, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) CALL ELOPID(W, IL, AP, IPP, IP, IDL, HXX, HYY, DIR) C C CALCULATE STEP LENGTH C PAP = 0. DO 60 L=1,IP PAP = PAP + P(L)*AP(L) 60 CONTINUE ALPHA = RR/PAP C C CALCULATE NEW ITERATE AND RESIDUAL C RROLD = RR RR = 0. DO 70 L=1,IP SOL(L) = SOL(L) + ALPHA*P(L) RHS(L) = RHS(L) - ALPHA*AP(L) RR = RR + RHS(L)*RHS(L) 70 CONTINUE BETA = RR/RROLD C C TERMINATE IF ANSWER SUFFICIENTLY ACCURATE C R = SQRT(RR/FIP) IF (ACCY.LT.0.0 .AND. KIT.EQ.1) R1 = R IF (TEST) WRITE (IOUT,99999) KIT, ALPHA, BETA, PAP, R IF (R.LT.ABS(ACCY)*R1) GO TO 110 IF (R.GT.1E5) GO TO 100 C C CALCULATE NEW STEP DIRECTION C DO 80 L=1,IP P(L) = RHS(L) + BETA*P(L) 80 CONTINUE 90 CONTINUE 100 ITMAX = 999 110 CONTINUE C C ********************** C COMPUTE FINAL SOLUTION C ********************** C CALL DIPOLE(SOL, IPP, IP, AP, IL, IDL, HXX, HYY, DIR) MODE = 3 IF (.NOT.LAPL) GO TO 140 DO 130 J=1,M DO 120 I=1,N F(I,J) = 0.0 120 CONTINUE 130 CONTINUE 140 CONTINUE CALL STRIP(AP, KC, LC, IL, W, IC, JC, IP, .FALSE., F, NDIM, N, M, * CH2, MODE, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) C IF (TEST) WRITE (IOUT,99997) RETURN C 99999 FORMAT (I5, 4E15.3) 99998 FORMAT (/1X, 49H ITER NO ALPHA BETA PAP , * 9H RR/) 99997 FORMAT (/1X, 34H THE SOLUTION IS STORED IN ARRAY F//) END SUBROUTINE COMPRS(IC, JC, NR, KC, LC, IL, INR, NO, IDL, IN, HXX, COM 10 * HYY, IPP, IP, DIR) DIMENSION IC(1), JC(1), KC(1), LC(1), INR(1), IDL(1), HXX(IPP), * HYY(IPP) LOGICAL SANT(4) LOGICAL DIR C C THIS SUBROUTINE PROCESSES THE INFORMATION ABOUT COORDINATES OF THE C IRREGULAR POINTS INTO A MORE COMPACT FORM. IC AND JC ARE EXTENDED C TO ALSO CONTAIN COORDINATES OF THE INTERNAL NEIGHBORS OF THE C IRREGULAR POINTS. VECTORS KC AND LC ARE CREATED. THE FIRST IP C ELEMENTS ARE IDENTICAL WITH THOSE OF IC AND JC. THE REST C CONTAIN COORDINATES OF THE DISCRETE DIPOLES. VECTORS INR C AND IDL ALSO ARE GENERATED. C IP1 = IP - (IPP-IP) NR = IP IL = IP NO = 0 IN = 0 DO 120 I=1,IP DO 10 J=1,4 SANT(J) = .FALSE. 10 CONTINUE K = IC(I) L = JC(I) J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) II = SIGN(1.0,HXX(J)) JJ = SIGN(1.0,HYY(J)) KK = -II LL = -JJ C C X INCREMENTS C IF (HX.GT.1.0) GO TO 30 IF (I.LE.IP1) GO TO 20 IF (ABS(HXX(J+1)).GT.1.0) GO TO 20 SANT(KK+2) = .TRUE. 20 CONTINUE SANT(II+2) = .TRUE. 30 CONTINUE C C Y INCREMENTS C IF (HY.GT.1.0) GO TO 50 IF (I.LE.IP1) GO TO 40 IF (ABS(HYY(J+1)).GT.1.0) GO TO 40 SANT(LL+3) = .TRUE. 40 CONTINUE SANT(JJ+3) = .TRUE. 50 CONTINUE C C FIND AN INTERNAL NEIGHBOR C DO 70 J=1,4 KX = (K+J) - 2 LY = (L+J) - 3 IF (J.EQ.4) KX = K IF (J.EQ.1) LY = L IF (SANT(J)) GO TO 70 NO = NO + 1 INR(NO) = 0 C C CHECK WHETHER THE POINT IS ALREADY ON THE LIST OF COORDINATES. C DO 60 M=1,NR IF (KX.NE.IC(M)) GO TO 60 IF (LY.NE.JC(M)) GO TO 60 INR(NO) = M GO TO 70 60 CONTINUE NR = NR + 1 IC(NR) = KX JC(NR) = LY 70 CONTINUE C C FIND THE COORDINATES OF THE DISCRETE MONOPOLES AND DIPOLES. C KC(I) = K LC(I) = L IF (.NOT.DIR) GO TO 120 KX = K + II LY = L + JJ C C FIND AN EXTERNAL POINT OF A DIPOLE C DO 110 J=1,2 IN = IN + 1 IDL(IN) = 0 IF (I.EQ.1) GO TO 90 IPP1 = IP + 1 C C CHECK WHETHER THE POINT IS ALREADY ON THE LIST OF COORDINATES. C DO 80 M=IPP1,IL IF (KX.NE.KC(M)) GO TO 80 IF (LY.NE.LC(M)) GO TO 80 IDL(IN) = M GO TO 100 80 CONTINUE 90 IL = IL + 1 KC(IL) = KX LC(IL) = LY 100 CONTINUE IF (HX.LE.HY) LY = L IF (HX.GT.HY) KX = K 110 CONTINUE C 120 CONTINUE RETURN END SUBROUTINE RIGHT(F, NDIM, N, M, RX, RY, RHS, IC, JC, IPP, IP, RIG 10 * HXX, HYY, TNX, TNY, H, CH2, DIR, LAPL) LOGICAL DIR, LAPL DIMENSION F(NDIM,M), IC(IP), JC(IP), TNX(IPP), TNY(IPP) DIMENSION RX(IPP), RY(IPP), RHS(IP), HXX(IPP), HYY(IPP) C C THIS SUBROUTINE COMPUTES RIGHT HAND SIDE TO THE CAPACITANCE SYSTEM C T C RHS = U F C USING BOUNDARY DATA STORED IN RX AND RY AND APPLYING SHORTLEY- C WELLER 5-POINT STENCIL AT THE IRREGULAR POINT. C IP1 = IP - (IPP-IP) SPOT = 0.0 DO 110 I=1,IP J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) C C X INCREMENTS C S = 0.0 T = 0.0 IF (HX.GT.1.0) GO TO 40 IF (.NOT.DIR) TX = TNX(J) IF (I.LE.IP1) GO TO 20 HX2 = ABS(HXX(J+1)) IF (HX2.GT.1.0) GO TO 20 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE X DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 10 TX2 = TNX(J+1) XX = 2.0/(HX+HX2) DIAG = XX*(TX+TX2) S = -(H*XX)*SQRT(1.0+TX**2) T = -(H*XX)*SQRT(1.0+TX2**2)*RX(J+1) GO TO 50 10 CONTINUE DIAG = 2.0/(HX*HX2) S = 2.0/(HX*(HX+HX2)) T = RX(J+1)*2.0/(HX2*(HX+HX2)) GO TO 50 20 CONTINUE C C ONLY ONE NEIGHBOR IN THIS DIRECTION C IF (DIR) GO TO 30 XX = 2.0/(1.0+HX) DIAG = XX*(1.0+TX) S = -(H*XX)*SQRT(1.0+TX**2) GO TO 50 30 CONTINUE DIAG = 2.0/HX S = 2.0/(HX*(HX+1.0)) GO TO 50 40 DIAG = 2.0 50 R = S*RX(J) + T C C Y INCREMENTS C S = 0.0 T = 0.0 IF (HY.GT.1.0) GO TO 90 IF (.NOT.DIR) TY = TNY(J) IF (I.LE.IP1) GO TO 70 HY2 = ABS(HYY(J+1)) IF (HY2.GT.1.0) GO TO 70 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE Y DIRECTION, AT C AT A DISTANCE LESS THAN OR EQUAL TO H. C IF (DIR) GO TO 60 TY2 = TNY(J+1) YY = 2.0/(HY+HY2) DIAG = YY*(TY+TY2) + DIAG S = -(H*YY)*SQRT(1.0+TY**2) T = -(H*YY)*SQRT(1.0+TY2**2)*RY(J+1) GO TO 100 60 CONTINUE DIAG = 2.0/(HY*HY2) + DIAG S = 2.0/(HY*(HY+HY2)) T = RY(J+1)*2.0/(HY2*(HY+HY2)) GO TO 100 70 CONTINUE C C ONLY ONE NEIGHBOR IN THIS DIRECTION C IF (DIR) GO TO 80 YY = 2.0/(1.0+HY) DIAG = YY*(1.0+TY) + DIAG S = -(H*YY)*SQRT(1.0+TY**2) GO TO 100 80 CONTINUE DIAG = 2.0/HY + DIAG S = 2.0/(HY*(HY+1.0)) GO TO 100 90 DIAG = 2.0 + DIAG 100 R = S*RY(J) + T + R C C FORM THE ELEMENT OF RHS C SC = 1.0/(DIAG+CH2) K = IC(I) L = JC(I) IF (.NOT.LAPL) SPOT = F(K,L) RHS(I) = (R+SPOT)*SC 110 CONTINUE RETURN END SUBROUTINE STEN5(ENT, NR, OUT, IPP, IP, INR, HXX, HYY, TNX, TNY, STE 10 * CH2, DIR) LOGICAL DIR DIMENSION ENT(NR), OUT(IP), INR(1), HXX(IP), HYY(IP), TNX(IP), * TNY(IP) DIMENSION S(4) C T C THIS SUBROUTINE COMPUTES OUT = (U A) ENT C T C U A REPRESENTS THE RESTRICTION OF THE DISCRETE LAPLACIAN C TO THE IRREGULAR MESH POINTS USING SHORTLEY-WELLER 5-POINT C STENCIL AT THE BOUNDARY OF THE REGION, I.E. THE INFORMATION C AT THE IRREGULAR POINT AND ITS INTERNAL NEIGHBORS. C THE EQUATIONS ARE SCALED SO THAT THE MAIN DIAGONAL C ELEMENTS ARE 1.0. C IP1 = IP - (IPP-IP) IN = IP NI = IP NA = 0 DO 140 I=1,IP DO 10 J=1,4 S(J) = 0.0 10 CONTINUE J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 II = SIGN(1.0,HXX(J)) JJ = SIGN(1.0,HYY(J)) KK = -II LL = -JJ HX = ABS(HXX(J)) HY = ABS(HYY(J)) C C X INCREMENTS C IF (HX.GT.1.0) GO TO 50 IF (.NOT.DIR) TX = TNX(J) IF (I.LE.IP1) GO TO 30 HX2 = ABS(HXX(J+1)) IF (HX2.GT.1.0) GO TO 30 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE X DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 20 TX2 = TNX(J+1) XX = 2.0/(HX+HX2) DIAG = XX*(TX+TX2) S(LL+3) = -XX*(TX+TX2) GO TO 60 20 DIAG = 2.0/(HX*HX2) GO TO 60 C C ONLY ONE NEIGHBOR IN THIS DIRECTION C 30 XX = 2.0/(1.0+HX) S(KK+2) = -XX IF (DIR) GO TO 40 DIAG = XX*(1.0+TX) S(LL+3) = -TX*XX GO TO 60 40 DIAG = 2.0/HX GO TO 60 50 DIAG = 2.0 S(1) = -1.0 S(3) = -1.0 60 CONTINUE C C Y INCREMENTS C IF (HY.GT.1.0) GO TO 100 IF (.NOT.DIR) TY = TNY(J) IF (I.LE.IP1) GO TO 80 HY2 = ABS(HYY(J+1)) IF (HY2.GT.1.0) GO TO 80 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE Y DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 70 TY2 = TNY(J+1) YY = 2.0/(HY+HY2) DIAG = YY*(TY+TY2) + DIAG S(KK+2) = S(KK+2) - YY*(TY+TY2) GO TO 110 70 DIAG = 2.0/(HY*HY2) + DIAG GO TO 110 C C ONLY ONE NEIGHBOR IN THIS DIRECTION C 80 YY = 2.0/(1.0+HY) S(LL+3) = S(LL+3) - YY IF (DIR) GO TO 90 DIAG = YY*(1.0+TY) + DIAG S(KK+2) = S(KK+2) - TY*YY GO TO 110 90 DIAG = 2.0/HY + DIAG GO TO 110 100 DIAG = 2.0 + DIAG S(LL+3) = S(LL+3) - 1.0 S(JJ+3) = -1.0 110 CONTINUE C C FIND AN INTERNAL NEIGHBOR C SC = 1.0/(DIAG+CH2) R = 0.0 DO 130 J=1,4 IF (S(J).EQ.0.0) GO TO 130 NA = NA + 1 IN = NI + 1 NI = IN IF (INR(NA).EQ.0) GO TO 120 IN = INR(NA) NI = NI - 1 120 CONTINUE R = ENT(IN)*S(J) + R 130 CONTINUE C C FORM THE ELEMENT OF OUT C OUT(I) = ENT(I) + R*SC 140 CONTINUE RETURN END SUBROUTINE NETS5(ENT, IPP, IP, OUT, NR, INR, HXX, HYY, TNX, TNY, NET 10 * CH2, DIR) LOGICAL DIR DIMENSION ENT(IP), OUT(NR), INR(1), HXX(IP), HYY(IP), TNX(IP), * TNY(IP) DIMENSION S(4) C T C THIS SUBROUTINE COMPUTES OUT = (U A) ENT C T C U A REPRESENTS THE EXTENSION OF THE DISCRETE LAPLACIAN C TO THE IRREGULAR MESH POINTS USING SHORTLEY-WELLER 5-POINT C STENCIL AT THE BOUNDARY OF THE REGION, I.E. THE INFORMATION C AT THE IRREGULAR POINT AND ITS INTERNAL NEIGHBORS. C THE EQUATIONS ARE SCALED SO THAT THE MAIN DIAGONAL C ELEMENTS ARE 1.0. C NETS5 IS A TRANSPOSE OF STEN5. C IP1 = IP - (IPP-IP) DO 10 I=1,NR OUT(I) = 0.0 10 CONTINUE IN = IP NI = IP NA = 0 DO 150 I=1,IP DO 20 J=1,4 S(J) = 0.0 20 CONTINUE J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) II = SIGN(1.0,HXX(J)) JJ = SIGN(1.0,HYY(J)) KK = -II LL = -JJ C C X INCREMENTS C IF (HX.GT.1.0) GO TO 60 IF (.NOT.DIR) TX = TNX(J) IF (I.LE.IP1) GO TO 40 HX2 = ABS(HXX(J+1)) IF (HX2.GT.1.0) GO TO 40 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE X DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 30 TX2 = TNX(J+1) XX = 2.0/(HX+HX2) DIAG = XX*(TX+TX2) S(LL+3) = -XX*(TX+TX2) GO TO 70 30 DIAG = 2.0/(HX*HX2) GO TO 70 C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C 40 XX = 2.0/(1.0+HX) S(KK+2) = -XX IF (DIR) GO TO 50 DIAG = XX*(1.0+TX) S(LL+3) = -TX*XX GO TO 70 50 DIAG = 2.0/HX GO TO 70 60 DIAG = 2.0 S(1) = -1.0 S(3) = -1.0 70 CONTINUE C C Y INCREMENTS C IF (HY.GT.1.0) GO TO 110 IF (.NOT.DIR) TY = TNY(J) IF (I.LE.IP1) GO TO 90 HY2 = ABS(HYY(J+1)) IF (HY2.GT.1.0) GO TO 90 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE Y DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 80 TY2 = TNY(J+1) YY = 2.0/(HY+HY2) DIAG = YY*(TY+TY2) + DIAG S(KK+2) = S(KK+2) - YY*(TY+TY2) GO TO 120 80 DIAG = 2.0/(HY*HY2) + DIAG GO TO 120 C C ONLY ONE NEIGHBOR IN THIS DIRECTION C 90 YY = 2.0/(1.0+HY) S(LL+3) = S(LL+3) - YY IF (DIR) GO TO 100 DIAG = YY*(1.0+TY) + DIAG S(KK+2) = S(KK+2) - TY*YY GO TO 120 100 DIAG = 2.0/HY + DIAG GO TO 120 110 DIAG = 2.0 + DIAG S(LL+3) = S(LL+3) - 1.0 S(JJ+3) = -1.0 120 CONTINUE C C FIND AN INTERNAL NEIGHBOR C R = ENT(I)/(DIAG+CH2) DO 140 J=1,4 IF (S(J).EQ.0.0) GO TO 140 NA = NA + 1 IN = NI + 1 NI = IN IF (INR(NA).EQ.0) GO TO 130 IN = INR(NA) NI = NI - 1 130 CONTINUE C C FORM THE ELEMENT OF OUT C OUT(IN) = OUT(IN) + R*S(J) 140 CONTINUE C C FORM THE ELEMENT OF OUT C OUT(I) = OUT(I) + ENT(I) 150 CONTINUE RETURN END SUBROUTINE DIPOLE(ENT, IPP, IP, OUT, IL, IDL, HXX, HYY, DIR) DIP 10 LOGICAL DIR DIMENSION ENT(IP), OUT(IL), IDL(1), HXX(IP), HYY(IP) DIMENSION S(2) C C THIS SUBROUTINE COMPUTES OUT = V *ENT C IL>=IP, I.E. THE LENGTH OF ENT IS NOT LESS THAN THE ONE OF OUT C V REPRESENTS THE DIPOLES (MONOPOLES) AT THE BOUNDARY OF THE REGION C IP1 = IP - (IPP-IP) IF (DIR) GO TO 20 C C MONOPOLES C DO 10 I=1,IP OUT(I) = ENT(I) 10 CONTINUE RETURN 20 CONTINUE C C DIPOLES C IN = IP NI = IP NA = 0 DO 30 I=1,IL OUT(I) = 0.0 30 CONTINUE DO 60 I=1,IP J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) IF (HX.LE.HY) TG = HX/HY IF (HX.GT.HY) TG = HY/HX S(1) = -TG S(2) = TG - 1.0 R = ENT(I) OUT(I) = OUT(I) + R C C FIND AN EXTERNAL NEIGHBOR C DO 50 J=1,2 NA = NA + 1 IN = NI + 1 NI = IN IF (IDL(NA).EQ.0) GO TO 40 IN = IDL(NA) NI = NI - 1 40 CONTINUE OUT(IN) = OUT(IN) + R*S(J) 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE ELOPID(ENT, IL, OUT, IPP, IP, IDL, HXX, HYY, DIR) ELO 10 LOGICAL DIR DIMENSION ENT(IL), OUT(IP), IDL(1), HXX(IP), HYY(IP) DIMENSION S(2) C T C THIS SUBROUTINE COMPUTES OUT = V *ENT C V REPRESENTS THE DIPOLES (MONOPOLES) AT THE BOUNDARY OF THE REGION C ELOPID IS A TRANSPOSE OF DIPOLE. C IL>=IP, I.E. THE LENGTH OF OUT IS NOT LESS THAN THE ONE OF ENT C IP1 = IP - (IPP-IP) IF (DIR) GO TO 20 C C MONOPOLES C DO 10 I=1,IP OUT(I) = ENT(I) 10 CONTINUE RETURN 20 CONTINUE C C DIPOLES C IN = IP NI = IP NA = 0 DO 50 I=1,IP J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) IF (HX.LE.HY) TG = HX/HY IF (HX.GT.HY) TG = HY/HX S(1) = -TG S(2) = TG - 1.0 R = ENT(I) C C FIND AN EXTERNAL NEIGHBOR C DO 40 J=1,2 NA = NA + 1 IN = NI + 1 NI = IN IF (IDL(NA).EQ.0) GO TO 30 IN = IDL(NA) NI = NI - 1 30 CONTINUE R = R + S(J)*ENT(IN) 40 CONTINUE C C FORM THE ELEMENT OF OUT C OUT(I) = R 50 CONTINUE RETURN END SUBROUTINE CMMEXP(INTL, IC, JC, HX, HY, TX, TY, BX, BY, F, NDIM, CMM 10 * N, M, H, CON, DIR, LAPL, TEST, C, IPMAX, IP, IPP, W, LEW, IW, * LEIW, IOUT, IERROR) REAL C(IP,IP), F(NDIM,M), W(LEW), HX(IPP), HY(IPP), TX(IPP), * TY(IPP), BX(IPP), BY(IPP) INTEGER IW(LEIW), IC(IP), JC(IP) LOGICAL DIR, LAPL, TEST C C C THIS PROGRAM SOLVES THE DIRICHLET AND NEUMANN PROBLEMS C FOR THE HELMHOLTZ EQUATION OVER AN ARBITRARY BOUNDED PLANAR REGION C IMBEDDED IN A RECTANGLE A,B;C,D : C C ( -L + C )U = F IN THE REGION, C C U = G C OR ON THE BOUNDARY, C DU/DN = G C C WHERE L IS THE LAPLACIAN, DU/DN IS THE NORMAL DERIVATIVE, C C IS A REAL C C THE ALGORITHM IS A VARIANT OF THE CAPACITANCE MATRIX METHOD. C FOR THE DETAILED DESCRIPTION OF THE METHOD SEE C W.PROSKUROWSKI AND O.WIDLUND - MATH.COMP., V.30(1976), PP.434-468. C C SUBROUTINE CMMEXP IS DESIGNED FOR REPETITIVE USE WITH MULTIPLE C RIGHT SIDES F AND FOR RELATIVELY SMALL PROBLEMS ( APPROXIMATELY UP C TO N*M = 64*64 AND IP = 150 ASSUMING THE CORE STORAGE OF 65K WORDS), C IF THE MEMORY SPACE REQUIRED DOES NOT EXCEED THE FAST CORE MEMORY. C AS A ONE SHOT SOLVER AND LARGE MESHES ONE SHOULD USE THE SUBROUTINE C CMMIMP INSTEAD. C C FOR PROBLEMS WITH MULTIPLE RIGHT SIDES ONLY THE SOLVING STAGE C (WITH INTL>0) SHOULD BE REPEATED AND THE USER MUST RECOMPUTE C THE VALUES OF F AND BX,BY ACCORDINGLY. C C THE EMBEDDING RECTANGLE IS COVERED WITH A SQUARE MESH OF LENGTH H. C IT SHOULD BE CHOSEN TO FOLLOW CLOSELY THE SHAPE OF THE REGION C IN SUCH A WAY THAT AT LEAST ONE MESH LENGTH IS LEFT C BETWEEN THE BOUNDARIES OF THE REGION AND THE RECTANGLE. C C THE TOTAL STORAGE FOR ARRAYS IS ABOUT IP*(IP+28) + N*(M+5). C FOR MORE DETAILED MEMORY SPACE REQUIREMENT C SEE PARAMETERS C,F,LEW,LEIW AND IPMAX. C C THIS PROGRAM WAS DEVELOPED BY W.PROSKUROWSKI. C THIS IS A JULY 1981 VERSION. C C THE ARGUMENTS ARE DEFINED AS : C C C * * * * * * * * * * ON INPUT * * * * * * * * * * * * * * * * C C C INTL C C = 0: ON INITIAL ENTRY TO CMMEXP; PREPROCESSING STAGE ONLY. C > 0: SOLVING STAGE; CAN BE USED FOR REPETITIVE USE OF C CMMEXP WITH MULTIPLE RIGHT SIDES, IF THE CONSTANT CON, C THE GEOMETRY OF THE PROBLEM AND THE MESH ARE UNCHANGED. C C IC,JC C C INTEGER VECTORS OF LENGTH IP CONTAINING COORDINATES (I,J) C OF ALL IRREGULAR POINTS ENUMERATED IN AN ARBITRARY ORDER. C THE VALUES OF I MUST BE LARGER THAN 1 AND LESS THAN N. C THE VALUES OF J MUST BE LARGER THAN 1 AND LESS THAN M. C A CALL OF CMMEXP DOES NOT DESTROY THE VALUES OF IC AND JC. C C HX,HY C C REAL VECTORS OF LENGTH IPP WITH SIGNED SHORTEST C DISTANCES IN BOTH COORDINATE DIRECTIONS FROM THE BOUNDARY C TO AN IRREGULAR POINT (I,J) IN UNITS OF STEP SIZE H. C IF (I,J) HAS A NEIGHBOR OUTSIDE THE REGION THEN THE VALUES OF C HX OR HY MUST BELONG TO THE INTERVAL -1.0 TO +1.0 (EXCLUDING 0.0 C OTHERWISE THE PROPER VALUE OF MAGNITUDE LARGER THAN 1.0 C MUST BE PUT IN. C A CALL OF CMMEXP DOES NOT DESTROY THE VALUES OF HX AND HY. C C TX,TY C C REAL VECTORS OF LENGTH IPP WITH ABSOLUTE VALUES OF C TANGENTS OF THE ANGLE BETWEEN THE NORMAL AND THE X-AXIS C AT (I+HX,J) ( AND THE Y-AXIS AT (I,J+HY) ). C TX,TY ARE USED ONLY FOR A PROBLEM WITH THE NEUMANN BOUNDARY C CONDITIONS. THE VALUES OF TX,TY MUST BE POSITIVE. C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0, A DUMMY C VALUE OF ARBITRARY MAGNITUDE CAN BE PUT IN. C A CALL OF CMMEXP DOES NOT DESTROY THE VALUES OF TX AND TY. C C BX,BY C C REAL VECTORS OF LENGTH IPP WITH BOUNDARY VALUES C AT (IC(I)+HX,JC(I)) AND (IC(I),JC(I)+HY). C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0 THEN C A DUMMY VALUE OF BX OR BY SHOULD BE PUT IN. C A CALL OF CMMEXP DESTROYS THE VALUES OF BX AND BY. C C IP C C THE NUMBER OF IRREGULAR MESH POINTS, I.E., POINTS C STRICTLY INSIDE THE REGION WHICH HAVE AT LEAST ONE C NEIGHBOR OUTSIDE OR ON THE BOUNDARY. C C IPP C C LENGTH OF VECTORS HX,HY,TX,TY,BX,BY. C C THE SET OF IRREGULAR POINTS CONSISTS OF THE SET (1) OF IP1 C POINTS WITH AT MOST 2 NEIGHBORS OUTSIDE THE REGION (OR ON C THE BOUNDARY) AND THE REMAINING SET (2) OF IP2 POINTS WITH C EXACTLY 3 NEIGHBORS OUTSIDE THE REGION (OR ON THE BOUNDARY). C THUS IP = IP1 + IP2. C IF A POINT BELONGS TO THE SET (2) THEN THE DISTANCES FROM C THE BOUNDARY MUST BE STORED IN A PAIR OF LOCATIONS: C (HX(I),HX(I+1)) AND (HY(I),HY(I+1)). NOTE THAT THE SECOND ENTRY C MUST BE LARGER IN MAGNITUDE THAN THE FIRST. A SIMILAR SITUATION C OCCURS FOR THE BOUNDARY VALUES BX,BY AND THE TANGENTS TX,TY. C THUS THE LENGTH OF THESE VECTORS MUST BE IPP = IP1 + IP2*2. C C IPMAX C C IS AN A PRIORI ESTIMATE OF IP AND THE DIMENSION OF VECTORS C IC,JC,HX,HY,TX,TY,BX,BY AS IT APPEARS IN THE PROGRAM CALLING C CMMEXP. IPMAX MUS BE AT LEAST IPP. C C C C C REAL ARRAY OF DIMENSION IP BY IP CONTAINING THE CAPACITANCE C MATRIX ( LATER OVERWRITTEN BY ITS FACTORED FORM ). C IN THE PROGRAM CALLING CMMEXP C MUST BE DIMENSIONED AS C A VECTOR OF THE LENGTH IPMAX*IPMAX. C C F C C REAL ARRAY OF ORDER N BY M CONTAING THE VALUES OF THE RIGHT C HAND SIDE OF THE HELMHOLTZ EQUATION FOR THE GRID POINTS C INSIDE THE IRREGULAR REGION SUPPLEMENTED BY SOME ARBITRARY C VALUES, FOR INSTANCE BY ZERO VALUES, FOR THE GRID POINTS C OUTSIDE THE IRREGULAR REGION. C F DOES NOT NEED TO BE INITIALIZED IF LAPL=.TRUE. C C NDIM C C THE FIRST DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING CMMEXP. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. NDIM MUST BE AT LEAST N. C C N C C THE NUMBER OF MESH LENGTHS IN THE X DIRECTION. C N MUST BE A POWER OF 2 AND NOT LESS THAN 8. C C M C C THE NUMBER OF MESH LENGTHS IN THE Y DIRECTION. C M SHOULD APPROXIMATELY BE M = N * (D-C)/(B-A). C WHERE A,B IS THE RANGE OF X IN THE RECTANGULAR REGION, C I.E., A =< X =< B, AND C,D IS THE RANGE OF Y C IN THE RECTANGULAR REGION, I.E., C =< Y =< D. C M MUST BE GREATER THAN OR EQUAL TO N. C C H C C THE MESH SIZE IN BOTH COORDINATE DIRECTIONS. C THE VALUE USED IN THIS PACKAGE IS H=(B-A)/N. C C CON C C REAL VARIABLE, THE CONSTANT IN THE HELMHOLTZ EQUATION. C C DIR C C LOGICAL VARIABLE, EQUAL TO .TRUE. FOR A DIRICHLET PROBLEM, C EQUAL TO .FALSE. FOR A NEUMANN PROBLEM. C NOTE THAT THE NEUMANN PROBLEM FOR THE LAPLACE EQUATION, I.E., C WITH DIR=.FALSE. AND CON=0.0, IS SINGULAR. IN SUCH A CASE THE C SOLUTION IS COMPUTED FOR A SLIGHTLY PERTURBED PROBLEM WITH C CON=8*(MACHINE EPSILON)*N**2. C C LAPL C C LOGICAL VARIABLE, EQUAL TO .TRUE. IF THE LAPLACE EQUATION C IS TO BE SOLVED, .FALSE. OTHERWISE, I.E., IF THE RIGHT SIDE C F IS NONZERO. A CALL WITH LAPL=.TRUE. IS SIGNIFICANTLY LESS C EXPENSIVE. C C TEST C C LOGICAL VARIABLE, EQUAL TO .TRUE. FOR TESTING OR DEBUGGING C PURPOSES, .FALSE. OTHERWISE,I.E., IN THE PRODUCTION MODE. C C W C C A REAL ARRAY THAT MUST BE PROVIDED FOR WORKSPACE. C C LEW C C THE TOTAL DIMENSION OF THE WORKSPACE W. C LEW MUST BE NOT LESS THAN 8*IP C C IW C C AN INTEGER ARRAY THAT MUST BE PROVIDED FOR WORKSPACE. C C LEIW C C THE TOTAL DIMENSION OF THE WORKSPACE IW. C LEIW MUST BE NOT LESS THAN 10*IP C C IOUT C C INTEGER VARIABLE, THE NUMBER OF THE OUTPUT DEVICE. C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * * * * * * * * * C C F C C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINTS INSIDE THE IRREGULAR REGION, C AND LARGELY ARBITRARY VALUES FOR POINTS OUTSIDE THE REGION. C C IERROR C C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. C IF IERROR IS GREATER THAN 0 NO SOLUTION IS SOUGHT. C C = 0: NO ERROR C = 1: N IS NOT A POWER OF 2 OR LESS THAN 16, OR M IS LESS THAN N C = 2: NDIM (THE FIRST DIMENSION OF ARRAY F) IS LESS THAN N. C = 3: IPMAX (THE ESTIMATE OF IP) IS TOO SMALL. C = 4: LEW OR LEIW (THE DIMENSIONS OF WORKSPACE) ARE TOO SMALL. C = 5: REGION TOO CLOSE TO THE EMBEDDING RECTANGLE. C = 6: INTERNAL ERROR IN THE PACKAGE. C = 7: CAPACITANCE MATRIX IS SINGULAR. C C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO CMMEXP, THE USER SHOULD TEST IERROR AFTER A CALL. C IF IERROR > 0, ONE CAN RUN CMMEXP AGAIN WITH TEST=.TRUE. C (DEBUGGING MODE) TO OBTAIN SOME ADDITIONAL INFORMATION. C C C ********************************************************************** C C IERROR = 0 IF (INTL.GT.0) GO TO 90 NN = 4 10 NN = NN*2 IF (N.GT.NN) GO TO 10 IF (N.LT.8 .OR. N.NE.NN .OR. M.LT.N) IERROR = 1 IF (IERROR.EQ.1) RETURN IF (TEST) WRITE (IOUT,99993) N, M, NDIM IF (N.LE.NDIM) GO TO 20 IERROR = 2 RETURN 20 CONTINUE IF (IPP.GT.IPMAX) IERROR = 3 IF (TEST) WRITE (IOUT,99998) IP, IPP, IPMAX LEN = MIN0(LEW/8,LEIW/10) IF (LEN.LT.IPP) IERROR = 4 FACTOR = FLOAT(IPP)/FLOAT(LEN) IF (TEST .AND. IERROR.EQ.4) WRITE (IOUT,99994) FACTOR IF (IERROR.GT.0) RETURN DO 30 I=1,IP K = IC(I) L = JC(I) IF (K.LE.1 .OR. K.GE.N .OR. L.LE.1 .OR. L.GE.N) IERROR = 5 30 CONTINUE C C CHECK THE INPUT VALUES BY PRINTING THEM OUT C IF (.NOT.TEST) GO TO 50 WRITE (IOUT,99992) IF (DIR) WRITE (IOUT,99991) IF (.NOT.DIR) WRITE (IOUT,99990) IP1 = IP - (IPP-IP) DO 40 I=1,IP J = I IF (J.GT.IP1) J = IP1 + (I-IP1)*2 - 1 IF (DIR) WRITE (IOUT,99989) I, IC(I), JC(I), HX(J), HY(J) IF (DIR .AND. I.GT.IP1) WRITE (IOUT,99988) HX(J+1), HY(J+1) IF (.NOT.DIR) WRITE (IOUT,99989) I, IC(I), JC(I), HX(J), HY(J), * TX(J), TY(J) IF (.NOT.DIR .AND. I.GT.IP1) WRITE (IOUT,99988) HX(J+1), * HY(J+1), TX(J+1), TY(J+1) IF (IC(I).LE.1 .OR. IC(I).GE.N .OR. JC(I).LE.1 .OR. JC(I).GE.N) * WRITE (IOUT,99987) I, IC(I), JC(I) 40 CONTINUE 50 CONTINUE IF (IERROR.GT.0) RETURN C C GUESS THE APPROXIMATE VALUES OF NR AND NO. C NR = IPMAX*2 NROLD = NR NO = IPMAX*3 NOOLD = NO I2 = NR + 1 I3 = I2 + NR C DO 60 I=1,IP IW(I) = IC(I) NRI = NR + I IW(NRI) = JC(I) 60 CONTINUE C C PREPROCESS THE BOUNDARY INFORMATION INTO A COMPACT FORM C AND COMPUTE THE VALUES OF NR AND NO. C CALL COMPCT(IW(1), IW(I2), NR, IW(I3), NO, HX, HY, IPP, IP) C IF (NR.GT.NROLD .OR. NO.GT.NOOLD) IERROR = 6 IF (TEST .AND. IERROR.EQ.6) WRITE (IOUT,99995) NR, NROLD, NO, * NOOLD IF (TEST) WRITE (IOUT,99996) NR, NO IF (IERROR.EQ.6) RETURN C C COMPRESS THE VECTOR IW. C L = IPMAX*2 DO 70 I=1,NR NRI = NR + I LI = L + I IW(NRI) = IW(LI) 70 CONTINUE K = NR*2 L = IPMAX*4 DO 80 I=1,NO KI = K + I LI = L + I IW(KI) = IW(LI) 80 CONTINUE 90 CONTINUE I2 = NR + 1 I3 = I2 + NR I4 = I3 + NO I5 = I4 + IP I6 = I5 + NR LI = I6 + N - 1 J2 = 1 + NR J3 = J2 + IP J4 = J3 + IP J5 = J4 + IP J6 = J5 + IP J7 = J6 + M J8 = J7 + M J9 = J8 + N J10 = J9 + N LW = J10 + N/4 - 1 IF (TEST .AND. INTL.EQ.0) WRITE (IOUT,99997) LW, LI, LEW, LEIW C C CALL THE MAIN SUBROUTINE C CALL HELMIT(INTL, F, NDIM, N, M, H, CON, DIR, LAPL, TEST, IW(1), * IW(I2), IW(I3), IW(I4), IW(I5), IW(I6), HX, HY, TX, TY, BX, BY, * W(1), W(J2), W(J3), W(J4), W(J5), W(J6), W(J7), W(J8), W(J9), * W(J10), C, IPP, IP, NR, NO, IOUT, IERROR) C IF (IERROR.EQ.7) WRITE (IOUT,99999) IF (IERROR.GT.0) RETURN C RETURN 99999 FORMAT (/1X, 31H CAPACITANCE MATRIX IS SINGULAR/) 99998 FORMAT (/1X, 4H IP=, I5, 5H IPP=, I5, 7H IPMAX=, I5/) 99997 FORMAT (/48H THE ACTUAL NEED FOR THE WORKSPACE W AND IW WAS, * 2I5/48H VS THE ESTIMATED SPACE RESERVED AS LEW AND LEIW, 2I5/) 99996 FORMAT (/1X, 35H THE LENGTH OF VECTORS IC,JC IS NR=, I5/1X, * 35H THE LENGTH OF VECTOR INR IS NO=, I5/) 99995 FORMAT (1X, 48H WRONG ESTIMATE OF THE INTERMEDIATE VALUES NR,NO/ * 43H CONSULT W.PROSKUROWSKI AT USC, LOS ANGELES) 99994 FORMAT (1X, 51H INSUFFICIENT SPACE RESERVED FOR THE WORKSPACE W,IW * /50H INCREASE THE DIMENSIONS LEW AND LEIW BY A FACTOR, 6H OF AT, * 6H LEAST, F5.2) 99993 FORMAT (1X, 3H N=, I3, 4H M=, I3, 11H AND NDIM=, I5) 99992 FORMAT (/1X, 24H THIS IS A TEST PRINTOUT/1X, 17H THE COORDINATES , * 24HOF THE IRREGULAR POINTS,/1X, 30H THE DISTANCES TO THE BOUNDARY * /1X, 40H AND TANGENTS (FOR THE NEUMANN PROBLEM) /1X, 9H FOR ALL , * 32HIRREGULAR POINTS ARE PRINTED OUT/) 99991 FORMAT (/40H NO I J HX HY /) 99990 FORMAT (/52H NO I J HX HY TX TY/) 99989 FORMAT (3I5, 4F10.5) 99988 FORMAT (15X, 4F10.5) 99987 FORMAT (/3I5, 21H BOUNDARIES TOO CLOSE) END C HEL 10 SUBROUTINE HELMIT(INTL, F, NDIM, N, M, H, CON, DIR, LAPL, TEST, HEL 20 * IC, JC, INR, IPS, INTER, IB, HXX, HYY, TNX, TNY, FX, FY, X, Y, * Z, U, W, RX, AIX, SAJN, COSAJN, SI, C, IPP, IP, NR, NO, IOUT, * IERROR) DIMENSION F(NDIM,M), C(IP,IP) DIMENSION IC(NR), JC(NR), INR(NO), IPS(IP), INTER(NR), IB(IP) DIMENSION HXX(IPP), HYY(IPP), TNX(IPP), TNY(IPP), FX(IPP), FY(IPP) DIMENSION X(NR), Y(IP), Z(IP), U(IP), W(IP), RX(M), AIX(M), * SAJN(N), COSAJN(N), SI(N) LOGICAL DIR, LAPL, TEST C DIMENSION KC(1), LC(1), ENT(1), DUM(1), IDUM(1) C C C PARAMETERS C ARE DESCRIBED BELOW SO FAR AS THEY ARE NOT INCLUDED OR DIFFER C FROM THOSE OF THE SUBROUTINE CMMEXP. C C IC,JC C INTEGER VECTORS OF LENGTH NR. CONTAIN COORDINATES C OF THE IRREGULAR POINTS AND THEIR NEIGHBORS INSIDE THE REGION. C C INR C INTEGER VECTOR OF LENGTH NO. CONTAINS A COMPRESSED C INFORMATION ABOUT THE IRREGULAR POINTS AND THEIR NEIGHBORS. C C X C REAL VECTOR OF LENGTH NR. C C Y,Z,U,W C REAL VECTORS OF LENGTH IP. C C IPS C INTEGER VECTOR OF LENGTH IP. C CONTAINS ROW NUMBERS FROM PIVOTING IN GAUSSIAN ELIMINATION. C C IIS C INTEGER VECTOR OF LENGTH IP ( USED IN STRIP ). C C IIT C INTEGER VECTOR OF LENGTH NR ( USED IN STRIP ). C C SAJN,COSAJN C REAL VECTORS OF LENGTH N ( USED IN STRIP ). C C IB C INTEGER VECTOR OF LENGTH N ( USED IN STRIP ). C C RX,AIX C REAL VECTORS OF LENGTH M ( USED IN STRIP ). C C SI C REAL VECTOR OF LENGTH N/4 ( USED IN STRIP ). C C NR C THE LENGTH OF VECTORS X AND IIT. APPROXIMATELY EQUAL TO 2*IP. C C NO C THE LENGTH OF VECTOR INR. APPROXIMATELY EQUAL TO 3*IP. C C C,Z,IPS,IC,JC,INR,HXX,HYY,TNX,TNY C CONTAIN INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF HELMIT C IS TO BE CALLED WITH MULTIPLE RIGHT SIDES (INTL>1). C C C ************************* C THE LIST OF SUBROUTINES : C ************************* C C RHSIDE - COMPUTES RIGHT HAND SIDE TO THE CAPACITANCE SYSTEM. C C STENCL - COMPUTES THE COEFFICIENTS FROM THE BOUNDARY STENCILS. C C CAPACI - COMPUTES THE CAPACITANCE MATRIX. C C COMPCT - PROCESSES THE COORDINATES OF THE IRREGULAR POINTS C INTO A COMPACT FORM. C C STENEX - MAPPING THE BOUNDARY STENCILS INTO THE IRREGULAR POINTS. C C STRIP - FAST HELMHOLTZ SOLVER. C C RFORT,FORT - FAST FOURIER TRANSFORMATION ROUTINES. C C DECOMP,SOLVE - GAUSSIAN ELIMINATION ROUTINES. C C IP1 = IP - (IPP-IP) HSQ = H*H CH2 = CON*HSQ IF (INTL.GT.0) GO TO 40 EPS = 1.0 10 EPS = 0.5*EPS EPSP1 = EPS + 1.0 IF (EPSP1.GT.1.0) GO TO 10 IF (.NOT.DIR .AND. ABS(CH2).LT.EPS*8.0) CH2 = EPS*8.0 C C *********************************************************** C INTL=0. PREPROCESSING PART IN WHICH THE CAPACITANCE C MATRIX IS GENERATED AND FACTORED. C *********************************************************** C C C COMPUTE AND FACTOR THE CAPACITANCE MATRIX C CALL STENCL(F, NDIM, N, M, IC, JC, DUM, IP, IPP, HXX, HYY, TNX, * TNY, X, Y, Z, U, W, CH2, DIR, .TRUE.) ENT(1) = 1.0 KC(1) = 2 LC(1) = 2 DO 30 J=1,M DO 20 I=1,N F(I,J) = 0.0 20 CONTINUE 30 CONTINUE CALL STRIP(ENT, KC, LC, 1, DUM, KC, LC, 1, .TRUE., F, NDIM, N, M, * CH2, 3, RX, AIX, SAJN, COSAJN, IB, SI, IDUM, IDUM, .FALSE.) CALL CAPACI(C, IP, IP, F, NDIM, N, X, Y, Z, U, IC, JC, W) CALL DECOMP(C, IP, IP, Z, IPS, IFLAG) IF (IFLAG.EQ.1) IERROR = 7 RETURN C C END OF PREPROCESSING STAGE C 40 CONTINUE C C *********************************************************** C INTL>0 : THE PREPROCESSING STAGE IS NOT REQUIRED, I.E. THE C CAPACITANCE MATRIX IS ALREADY GENERATED AND DECOMPOSED. C *********************************************************** C C C COMPUTE THE RIGHT HAND SIDE TO THE CAPACITANCE MATRIX SYSTEM. C CALL RHSIDE(F, NDIM, N, M, FX, FY, Y, IC, JC, IP, IPP, HXX, HYY, * TNX, TNY, H, CH2, DIR, LAPL) IF (LAPL) GO TO 80 C C COMPUTE THE SPACE POTENTIAL. C DO 60 J=1,M DO 50 I=1,N F(I,J) = F(I,J)*HSQ 50 CONTINUE 60 CONTINUE CALL STRIP(DUM, KC, LC, 1, X, IC, JC, NR, .FALSE., F, NDIM, N, M, * CH2, 2, RX, AIX, SAJN, COSAJN, IB, SI, IDUM, INTER, .FALSE.) CALL STENEX(X, NR, W, IP, IPP, INR, HXX, HYY, TNX, TNY, CH2, DIR) DO 70 I=1,IP Y(I) = Y(I) - W(I) 70 CONTINUE 80 CONTINUE C C COMPUTE THE CHARGE DISTRIBUTION BY BACKSOLVING THE CAPACITANCE C MATRIX SYSTEM. C CALL SOLVE(C, IP, IP, X, Y, Z, IPS) IF (.NOT.LAPL) GO TO 110 DO 100 J=1,M DO 90 I=1,N F(I,J) = 0.0 90 CONTINUE 100 CONTINUE 110 CONTINUE C C ********************** C COMPUTE FINAL SOLUTION C ********************** C CALL STRIP(X, IC, JC, IP, DUM, KC, LC, 1, .FALSE., F, NDIM, N, M, * CH2, 3, RX, AIX, SAJN, COSAJN, IB, SI, INTER, IDUM, .FALSE.) C RETURN END C RHS 10 SUBROUTINE RHSIDE(F, NDIM, N, M, RX, RY, RHS, IC, JC, IP, IPP, RHS 20 * HXX, HYY, TNX, TNY, H, CH2, DIR, LAPL) LOGICAL LAPL, DIR DIMENSION F(NDIM,M), IC(IP), JC(IP) DIMENSION RX(IPP), RY(IPP), RHS(IP), HXX(IPP), HYY(IPP), * TNX(IPP), TNY(IPP) C C THIS SUBROUTINE COMPUTES RIGHT HAND SIDE TO THE CAPACITANCE SYSTEM C T C RHS = U F C USING BOUNDARY DATA STORED IN RX AND RY AND APPLYING SHORTLEY- C WELLER 5-POINT STENCIL AT THE IRREGULAR POINT. C HSQ = H*H IP1 = IP - (IPP-IP) SPOT = 0.0 DO 110 I=1,IP J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) C C X INCREMENTS C S = 0.0 T = 0.0 IF (HX.GT.1.0) GO TO 40 IF (.NOT.DIR) TX = TNX(J) IF (I.LE.IP1) GO TO 20 HX2 = ABS(HXX(J+1)) IF (HX2.GT.1.0) GO TO 20 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE X DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 10 TX2 = TNX(J+1) XX = 2.0/(HX+HX2) DIAG = XX*(TX+TX2) S = -(H*XX)*SQRT(1.0+TX**2) T = -(H*XX)*SQRT(1.0+TX2**2)*RX(J+1) GO TO 50 10 CONTINUE DIAG = 2.0/(HX*HX2) S = 2.0/(HX*(HX+HX2)) T = RX(J+1)*2.0/(HX2*(HX+HX2)) GO TO 50 20 CONTINUE C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C IF (DIR) GO TO 30 XX = 2.0/(1.0+HX) DIAG = XX*(1.0+TX) S = -H*XX*SQRT(1.0+TX**2) GO TO 50 30 CONTINUE DIAG = 2.0/HX S = 2.0/(HX*(HX+1.0)) GO TO 50 40 DIAG = 2.0 50 R = S*RX(J) + T C C Y INCREMENTS C IF (.NOT.DIR) TY = TNY(J) S = 0.0 T = 0.0 IF (HY.GT.1.0) GO TO 90 IF (.NOT.DIR) TY = TNY(I) IF (I.LE.IP1) GO TO 70 HY2 = ABS(HYY(J+1)) IF (HY2.GT.1.0) GO TO 70 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE Y DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 60 TY2 = TNY(J+1) YY = 2.0/(HY+HY2) DIAG = YY*(TY+TY2) + DIAG S = -(H*YY)*SQRT(1.0+TY**2) T = -(H*YY)*SQRT(1.0+TY2**2)*RY(J+1) GO TO 100 60 CONTINUE DIAG = 2.0/(HY*HY2) + DIAG S = 2.0/(HY*(HY+HY2)) T = RY(J+1)*2.0/(HY2*(HY+HY2)) GO TO 100 70 CONTINUE C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C IF (DIR) GO TO 80 YY = 2.0/(1.0+HY) DIAG = YY*(1.0+TY) + DIAG S = -H*YY*SQRT(1.0+TY**2) GO TO 100 80 CONTINUE DIAG = 2.0/HY + DIAG S = 2.0/(HY*(HY+1.0)) GO TO 100 90 DIAG = 2.0 + DIAG 100 R = S*RY(J) + T + R C C FORM AN ELEMENT OF RHS. C SC = 1.0/(DIAG+CH2) K = IC(I) L = JC(I) IF (.NOT.LAPL) SPOT = F(K,L)*HSQ RHS(I) = (R+SPOT)*SC 110 CONTINUE RETURN END SUBROUTINE STENCL(F, NDIM, N, M, IC, JC, OUT, IP, IPP, HXX, HYY, STE 10 * TNX, TNY, AA, AB, AC, AD, SCAL, CH2, DIR, SAVE) DIMENSION F(NDIM,M), IC(IP), JC(IP), OUT(IP), HXX(IP), HYY(IP) DIMENSION AA(IP), AB(IP), AC(IP), AD(IP), SCAL(IP), TNX(IP), * TNY(IP) DIMENSION S(4) LOGICAL DIR, SAVE C T C THIS SUBROUTINE COMPUTES OUT = (U A) ENT C T C U A REPRESENTS THE RESTRICTION OF THE DISCRETE LAPLACIAN C TO THE IRREGULAR MESH POINTS USING SHORTLEY-WELLER 5-POINT C STENCIL AT THE BOUNDARY OF THE REGION, I.E. THE INFORMATION C AT THE IRREGULAR POINT AND ITS INTERNAL NEIGHBORS. C THE EQUATIONS ARE SCALED SO THAT THE MAIN DIAGONAL C ELEMENTS ARE 1.0. C IF SAVE=.TRUE. THE COEFFICIENTS FROM THE STENCILS AT THE C BOUNDARY ARE STORED AND USED LATER IN GENERATING THE C CAPACITANCE MATRIX. C IP1 = IP - (IPP-IP) DO 130 I=1,IP K = IC(I) L = JC(I) DO 10 J=1,4 S(J) = 0.0 10 CONTINUE J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 II = SIGN(1.0,HXX(J)) JJ = SIGN(1.0,HYY(J)) KK = -II LL = -JJ HX = ABS(HXX(J)) HY = ABS(HYY(J)) C C X INCREMENTS C IF (HX.GT.1.0) GO TO 50 IF (.NOT.DIR) TX = TNX(J) IF (I.LE.IP1) GO TO 30 HX2 = ABS(HXX(J+1)) IF (HX2.GT.1.0) GO TO 30 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE X DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 20 TX2 = TNX(J+1) XX = 2.0/(HX+HX2) DIAG = XX*(TX+TX2) S(LL+3) = -XX*(TX+TX2) GO TO 60 20 DIAG = 2.0/(HX*HX2) GO TO 60 C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C 30 XX = 2.0/(1.0+HX) S(KK+2) = -XX IF (DIR) GO TO 40 DIAG = XX*(1.0+TX) S(LL+3) = -TX*XX GO TO 60 40 CONTINUE DIAG = 2.0/HX GO TO 60 50 DIAG = 2.0 S(1) = -1.0 S(3) = -1.0 60 CONTINUE C C Y INCREMENTS C IF (HY.GT.1.0) GO TO 100 IF (.NOT.DIR) TY = TNY(J) IF (I.LE.IP1) GO TO 80 HY2 = ABS(HYY(J+1)) IF (HY2.GT.1.0) GO TO 80 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE Y DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 70 TY2 = TNY(J+1) YY = 2.0/(HY+HY2) DIAG = YY*(TY+TY2) + DIAG S(KK+2) = S(KK+2) - YY*(TY+TY2) GO TO 110 70 DIAG = 2.0/(HY*HY2) + DIAG GO TO 110 C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C 80 YY = 2.0/(1.0+HY) S(LL+3) = S(LL+3) - YY IF (DIR) GO TO 90 DIAG = YY*(1.0+TY) + DIAG S(KK+2) = S(KK+2) - TY*YY GO TO 110 90 CONTINUE DIAG = 2.0/HY + DIAG GO TO 110 100 DIAG = 2.0 + DIAG S(JJ+3) = -1.0 S(LL+3) = S(LL+3) - 1.0 110 CONTINUE SC = 1.0/(DIAG+CH2) IF (SAVE) GO TO 120 C C FORM THE ELEMENT OF OUT R = S(1)*F(K-1,L) + S(2)*F(K,L-1) + S(3)*F(K+1,L) + * S(4)*F(K,L+1) OUT(I) = F(K,L) + R*SC GO TO 130 120 CONTINUE C SAVE THE COEFFICIENTS FOR GENERATING THE CAPACITANCE MATRIX. SCAL(I) = SC AA(I) = S(1) AB(I) = S(2) AC(I) = S(3) AD(I) = S(4) 130 CONTINUE RETURN END C CAP 10 SUBROUTINE CAPACI(C, NDIMC, IP, F, NDIMF, N, AA, AB, AC, AD, IC, CAP 20 * JC, SCAL) DIMENSION C(NDIMC,IP) DIMENSION F(NDIMF,1) DIMENSION AA(1), AB(1), AC(1), AD(1), IC(1), JC(1), SCAL(1) C C THIS SUBROUTINE COMPUTES THE CAPACITANCE MATRIX C EXPLOITING C THE TRANSLATION INVARIANCY PROPERTY OF THE SOLUTION OF HELMHOLTZ C EQUATION ON AN INFINITE STRIP WITH PERIODIC BOUNDARY CONDITIONS. C C F REPRESENTS A DISCRETE FUNDAMENTAL SOLUTION, I.E. SOLUTION C TO THE PROBLEM WITH A UNIT CHARGE AT THE POINT (I0,J0). C I0 = 2 J0 = 2 C C SINGLE LAYER APPROACH C DO 20 I=1,IP K = IC(I) L = JC(I) SC = SCAL(I) AKML = AA(I) AKLM = AB(I) AKPL = AC(I) AKLP = AD(I) DO 10 J=1,IP KK = IC(J) LL = JC(J) K1 = IABS(K-KK) K2 = IABS(K-1-KK) K3 = IABS(K+1-KK) L1 = IABS(L-LL) + J0 L2 = IABS(L-1-LL) + J0 L3 = IABS(L+1-LL) + J0 M1 = MIN0(K1,N-K1) + I0 M2 = MIN0(K2,N-K2) + I0 M3 = MIN0(K3,N-K3) + I0 S = AKML*F(M2,L1) + AKLM*F(M1,L2) + AKPL*F(M3,L1) + * AKLP*F(M1,L3) C(I,J) = F(M1,L1) + S*SC 10 CONTINUE 20 CONTINUE RETURN END C COM 10 SUBROUTINE COMPCT(IC, JC, NR, INR, NO, HXX, HYY, IPP, IP) COM 20 DIMENSION IC(1), JC(1), INR(1), HXX(IPP), HYY(IPP) LOGICAL SANT(4) C THIS SUBROUTINE PROCESSES THE INFORMATION ABOUT COORDINATES OF THE C IRREGULAR POINTS INTO A MORE COMPACT FORM. C IC AND JC ARE EXTENDED TO ALSO CONTAIN COORDINATES OF THE C INTERNAL NEIGHBORS OF IRREGULAR POINTS. VECTOR INR IS C GENERATED. C IP1 = IP - (IPP-IP) NR = IP NO = 0 DO 80 I=1,IP DO 10 J=1,4 SANT(J) = .FALSE. 10 CONTINUE K = IC(I) L = JC(I) J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 HX = ABS(HXX(J)) HY = ABS(HYY(J)) II = SIGN(1.0,HXX(J)) JJ = SIGN(1.0,HYY(J)) KK = -II LL = -JJ C C X INCREMENTS C IF (HX.GT.1.0) GO TO 30 IF (I.LE.IP1) GO TO 20 IF (ABS(HXX(J+1)).GT.1.0) GO TO 20 SANT(KK+2) = .TRUE. 20 CONTINUE SANT(II+2) = .TRUE. 30 CONTINUE C C Y INCREMENTS C IF (HY.GT.1.0) GO TO 50 IF (I.LE.IP1) GO TO 40 IF (ABS(HYY(J+1)).GT.1.0) GO TO 40 SANT(LL+3) = .TRUE. 40 CONTINUE SANT(JJ+3) = .TRUE. 50 CONTINUE C C FIND AN INTERNAL NEIGHBOR C DO 70 J=1,4 KX = (K+J) - 2 LY = (L+J) - 3 IF (J.EQ.4) KX = K IF (J.EQ.1) LY = L IF (SANT(J)) GO TO 70 NO = NO + 1 INR(NO) = 0 C C CHECK WHETHER THE POINT IS ALREADY ON THE LIST OF COORDINATES. C DO 60 M=1,NR IF (KX.NE.IC(M)) GO TO 60 IF (LY.NE.JC(M)) GO TO 60 INR(NO) = M GO TO 70 60 CONTINUE NR = NR + 1 IC(NR) = KX JC(NR) = LY 70 CONTINUE 80 CONTINUE RETURN END C STE 10 SUBROUTINE STENEX(ENT, NR, OUT, IP, IPP, INR, HXX, HYY, TNX, TNY, STE 20 * CH2, DIR) DIMENSION ENT(NR), OUT(IP), INR(1), HXX(IP), HYY(IP), TNX(IP), * TNY(IP) DIMENSION S(4) LOGICAL DIR C T C THIS SUBROUTINE COMPUTES OUT = (U A) ENT C T C U A REPRESENTS THE RESTRICTION OF THE DISCRETE LAPLACIAN C TO THE IRREGULAR MESH POINTS USING SHORTLEY-WELLER 5-POINT C STENCIL AT THE BOUNDARY OF THE REGION, I.E. THE INFORMATION C AT THE IRREGULAR POINT AND ITS INTERNAL NEIGHBORS. C THE EQUATIONS ARE SCALED SO THAT THE MAIN DIAGONAL C ELEMENTS ARE 1.0. C IN CONTRAST TO THE SUBROUTINE STENCL, HERE THE RESULT IS STORED C IN A VECTOR OUT OF THE LENGTH IP. C IP1 = IP - (IPP-IP) IN = IP NI = IP NA = 0 DO 140 I=1,IP DO 10 J=1,4 S(J) = 0.0 10 CONTINUE J = I IF (I.GT.IP1) J = IP1 + (I-IP1)*2 - 1 II = SIGN(1.0,HXX(J)) JJ = SIGN(1.0,HYY(J)) KK = -II LL = -JJ HX = ABS(HXX(J)) HY = ABS(HYY(J)) C C X INCREMENTS C IF (HX.GT.1.0) GO TO 50 IF (.NOT.DIR) TX = TNX(J) IF (I.LE.IP1) GO TO 30 HX2 = ABS(HXX(J+1)) IF (HX2.GT.1.0) GO TO 30 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE X DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 20 TX2 = TNX(J+1) XX = 2.0/(HX+HX2) DIAG = XX*(TX+TX2) S(LL+3) = -XX*(TX+TX2) GO TO 60 20 DIAG = 2.0/(HX*HX2) GO TO 60 C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C 30 XX = 2.0/(1.0+HX) S(KK+2) = -XX IF (DIR) GO TO 40 DIAG = XX*(1.0+TX) S(LL+3) = -TX*XX GO TO 60 40 DIAG = 2.0/HX GO TO 60 50 DIAG = 2.0 S(1) = -1.0 S(3) = -1.0 60 CONTINUE C C Y INCREMENTS C IF (HY.GT.1.0) GO TO 100 IF (.NOT.DIR) TY = TNY(J) IF (I.LE.IP1) GO TO 80 HY2 = ABS(HYY(J+1)) IF (HY2.GT.1.0) GO TO 80 C C TWO NEIGHBORS, I.E. AT BOTH SIDES IN THE Y DIRECTION, AT C A DISTANCE LESS THAN OR EQUAL H. C IF (DIR) GO TO 70 TY2 = TNY(J+1) YY = 2.0/(HY+HY2) DIAG = YY*(TY+TY2) + DIAG S(KK+2) = S(KK+2) - YY*(TY+TY2) GO TO 110 70 DIAG = 2.0/(HY*HY2) + DIAG GO TO 110 C C ONLY ONE NEIGHBOR IN THIS DIRECTION. C 80 YY = 2.0/(1.0+HY) S(LL+3) = S(LL+3) - YY IF (DIR) GO TO 90 DIAG = YY*(1.0+TY) + DIAG S(KK+2) = S(KK+2) - TY*YY GO TO 110 90 DIAG = 2.0/HY + DIAG GO TO 110 100 DIAG = 2.0 + DIAG S(LL+3) = S(LL+3) - 1.0 S(JJ+3) = -1.0 110 CONTINUE C C FORM THE ELEMENT OF OUT C SC = 1.0/(DIAG+CH2) R = 0.0 DO 130 J=1,4 IF (S(J).EQ.0.0) GO TO 130 NA = NA + 1 IN = NI + 1 NI = IN IF (INR(NA).EQ.0) GO TO 120 IN = INR(NA) NI = NI - 1 120 CONTINUE R = ENT(IN)*S(J) + R 130 CONTINUE OUT(I) = ENT(I) + R*SC 140 CONTINUE RETURN END C CMM 10 SUBROUTINE CMMSIX(INTL, IC, JC, HX, HY, BX, BY, F, NDIM, N, M, H, CMM 20 * TEST, C, IPMAX, IP, W, LEW, IW, LEIW, NCOR, IOUT, IERROR) REAL C(IP,IP), F(NDIM,M), W(LEW), HX(IP), HY(IP), BX(IP), BY(IP) INTEGER IW(LEIW), IC(IP), JC(IP) LOGICAL TEST C C C THIS PROGRAM SOLVES THE DIRICHLET PROBLEM FOR THE POISSON EQUATION C OVER AN ARBITRARY BOUNDED PLANAR REGION IMBEDDED IN A RECTANGLE A,B C C C -LU = F IN THE REGION, C C U = G ON THE BOUNDARY, C C WHERE L IS THE LAPLACIAN, AND F AND G ARE GIVEN FUNCTIONS OF X AND Y C C THE ALGORITHM IS A VARIANT OF THE CAPACITANCE MATRIX METHOD. C FOR THE DETAILED DESCRIPTION OF THE METHOD SEE C V.PEREYRA,W.PROSKUROWSKI AND O.WIDLUND - MATH.COMP., V.31(1977), PP. C C SUBROUTINE CMMEXP IS DESIGNED TO COMPUTE THE SOLUTION WITH HIGH C ORDER ACCURACY. IT CAN BE USED FOR REPETITIVE USE WITH MULTIPLE C RIGHT SIDES F AND FOR RELATIVELY SMALL PROBLEMS ( APPROXIMATELY C UP TO N*M = 64*64 AND IP = 150 ), IF THE MEMORY SPACE REQUIRED C DOES NOT EXCEED THE FAST CORE MEMORY. FOR PROBLEMS WITH THE C NEUMANN BOUNDARY CONDITIONS OR NONZERO CONSTANTS (THE HELMHOLTZ C EQUATION) ONE SHOULD USE THE SUBROUTINE CMMEXP INSTEAD. AS A ONE C SHOT SOLVER AND LARGE MESHES ONE SHOULD USE THE SUBROUTINE C CMMIMP INSTEAD. C C FOR PROBLEMS WITH MULTIPLE RIGHT SIDES ONLY THE SOLVING STAGE C (WITH INTL>0) SHOULD BE REPEATED AND THE USER MUST RECOMPUTE C THE VALUES OF F AND BX,BY ACCORDINGLY. C C THE EMBEDDING RECTANGLE IS COVERED WITH A SQUARE MESH OF LENGTH H. C IT SHOULD BE CHOSEN TO FOLLOW CLOSELY THE SHAPE OF THE REGION C IN SUCH A WAY THAT AT LEAST NCOR MESH LENGTHS ARE LEFT C BETWEEN THE BOUNDARIES OF THE REGION AND THE RECTANGLE. C THE MESH LENGTH MUST BE CHOSEN IN SUCH A WAY THAT INSIDE THE C REGION THERE ARE AT LEAST (2*NCOR+1) POINTS ON EACH MESH LINE. C C THE TOTAL STORAGE FOR ARRAYS IS ABOUT IP*(IP+7+4*NCOR) + N*(3M+5). C FOR MORE DETAILED MEMORY SPACE REQUIREMENT C SEE PARAMETERS C,F,LEW,LEIW AND IPMAX. C C THIS PROGRAM WAS DEVELOPED BY W.PROSKUROWSKI. C THIS IS A JULY 1981 VERSION. C C THE ARGUMENTS ARE DEFINED AS : C C C * * * * * * * * * * ON INPUT * * * * * * * * * * * * * * * * C C C INTL C C = 0: ON INITIAL ENTRY TO CMMSIX; PREPROCESSING STAGE ONLY. C > 0: SOLVING STAGE; CAN BE USED FOR REPETITIVE USE OF C CMMSIX WITH MULTIPLE RIGHT SIDES, IF THE GEOMETRY C OF THE PROBLEM AND THE MESH ARE UNCHANGED. C C IC,JC C C INTEGER VECTORS OF LENGTH IP CONTAINING COORDINATES (I,J) C OF ALL IRREGULAR POINTS ENUMERATED IN AN ARBITRARY ORDER. C THE VALUES OF I MUST BE LARGER THAN 1 AND LESS THAN N. C THE VALUES OF J MUST BE LARGER THAN 1 AND LESS THAN M. C A CALL OF CMMSIX DOES NOT DESTROY THE VALUES OF IC AND JC. C C HX,HY C C REAL VECTORS OF LENGTH IP WITH SIGNED SHORTEST C DISTANCES IN BOTH COORDINATE DIRECTIONS FROM THE BOUNDARY C TO AN IRREGULAR POINT (I,J) IN UNITS OF THE STEP SIZE H. C IF (I,J) HAS A NEIGHBOR OUTSIDE THE REGION THEN THE VALUES OF C HX OR HY MUST BELONG TO THE INTERVAL -1.0 TO +1.0 (EXCLUDING 0.0 C OTHERWISE THE PROPER VALUE OF MAGNITUDE LARGER THAN 1.0 C MUST BE PUT IN. C A CALL OF CMMSIX DOES NOT DESTROY THE VALUES OF HX AND HY. C C BX,BY C C REAL VECTORS OF LENGTH IP WITH BOUNDARY VALUES C AT (IC(I)+HX,JC(I)) AND (IC(I),JC(I)+HY). C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0 THEN C A DUMMY VALUE OF BX OR BY SHOULD BE PUT IN. C A CALL OF CMMSIX DESTROYS THE VALUES OF BX AND BY. C C IP C C THE NUMBER OF IRREGULAR MESH POINTS, I.E., POINTS C STRICTLY INSIDE THE REGION WHICH HAVE AT LEAST ONE C NEIGHBOR OUTSIDE OR ON THE BOUNDARY. C C IPMAX C C IS AN A PRIORI ESTIMATE OF IP AND THE DIMENSION OF VECTORS C IC,JC,HX,HY,BX,BY AS IT APPEARS IN THE PROGRAM CALLING C CMMSIX. IPMAX MUS BE AT LEAST IP. C C C C C REAL ARRAY OF DIMENSION IP BY IP CONTAINING THE CAPACITANCE C MATRIX ( LATER OVERWRITTEN BY ITS FACTORED FORM ). C IN THE PROGRAM CALLING CMMSIX C MUST BE DIMENSIONED AS C A VECTOR OF THE LENGTH NOT LESS THAN IP*IP. C C F C C REAL ARRAY OF ORDER N BY M CONTAING THE VALUES OF THE RIGHT C HAND SIDE OF THE HELMHOLTZ EQUATION FOR THE GRID POINTS C INSIDE THE IRREGULAR REGION SUPPLEMENTED BY SOME ARBITRARY C VALUES, FOR INSTANCE BY ZERO VALUES, FOR THE GRID POINTS C OUTSIDE THE IRREGULAR REGION. C C NDIM C C THE FIRST DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING CMMSIX. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. NDIM MUST BE AT LEAST N. C C N C C THE NUMBER OF MESH LENGTHS IN THE X DIRECTION. C N MUST BE A POWER OF 2 AND NOT LESS THAN 8. C C M C C THE NUMBER OF MESH LENGTHS IN THE Y DIRECTION. C M SHOULD APPROXIMATELY BE M = N * (D-C)/(B-A). C WHERE A,B IS THE RANGE OF X IN THE RECTANGULAR REGION, C I.E., A =< X =< B, AND C,D IS THE RANGE OF Y C IN THE RECTANGULAR REGION, I.E., C =< Y =< D. C M MUST BE GREATER THAN OR EQUAL TO N. C C H C C THE MESH SIZE IN BOTH COORDINATE DIRECTIONS. C THE VALUE USED IN THIS PACKAGE IS H=(B-A)/N. C C TEST C C LOGICAL VARIABLE, EQUAL TO .TRUE. FOR TESTING OR DEBUGGING C PURPOSES, .FALSE. OTHERWISE,I.E., IN THE PRODUCTION MODE. C C W C C A REAL ARRAY THAT MUST BE PROVIDED FOR WORKSPACE. C C LEW C C THE TOTAL DIMENSION OF THE WORKSPACE W. C LEW MUST BE NOT LESS THAN 2*N*M + (8+4*NCOR)*IP. C C IW C C AN INTEGER ARRAY THAT MUST BE PROVIDED FOR WORKSPACE. C C LEIW C C THE TOTAL DIMENSION OF THE WORKSPACE IW. C LEIW MUST BE NOT LESS THAN 2*IP C C NCOR C C THE NUMBER OF DEFERRED CORRECTIONS REQUIRED. C NCOR MUST BE 0,1 OR 2. THE SCHEME IS ACCURATE C UP TO THE (NCOR+1)*2 ORDER. C C IOUT C C INTEGER VARIABLE, THE NUMBER OF THE OUTPUT DEVICE. C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * * * * * * * * * C C F C C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINTS INSIDE THE IRREGULAR REGION, C AND LARGELY ARBITRARY VALUES FOR POINTS OUTSIDE THE REGION. C C IERROR C C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. C IF IERROR IS GREATER THAN 0 NO SOLUTION IS SOUGHT. C C = 0: NO ERROR C = 1: N IS NOT A POWER OF 2 OR LESS THAN 16, OR M IS LESS THAN N C = 2: NDIM (THE FIRST DIMENSION OF ARRAY F) IS LESS THAN N. C = 3: IPMAX (THE ESTIMATE OF IP) IS TOO SMALL. C = 4: LEW OR LEIW (THE DIMENSIONS OF WORKSPACE) ARE TOO SMALL. C = 5: REGION TOO CLOSE TO THE EMBEDDING RECTANGLE. C = 6: INSUFFICIENT NUMBER OF MESH POINTS INSIDE THE REGION. C = 7: CAPACITANCE MATRIX IS SINGULAR. C C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO CMMSIX, THE USER SHOULD TEST IERROR AFTER A CALL. C IF IERROR > 0, ONE CAN RUN CMMSIX AGAIN WITH TEST=.TRUE. C (DEBUGGING MODE) TO OBTAIN SOME ADDITIONAL INFORMATION. C C C ********************************************************************** C C IERROR = 0 IF (INTL.GT.0) GO TO 170 NN = 4 10 NN = NN*2 IF (N.GT.NN) GO TO 10 IF (N.LT.8 .OR. N.NE.NN .OR. M.LT.N) IERROR = 1 IF (IERROR.EQ.1) RETURN IF (TEST) WRITE (IOUT,99990) N, M, NDIM IF (N.LE.NDIM) GO TO 20 IERROR = 2 RETURN 20 CONTINUE IF (IP.GT.IPMAX) IERROR = 3 IF (TEST) WRITE (IOUT,99993) IP, IPMAX IF (IERROR.GT.0) RETURN LEF = LEW - 2*N*M NO = 0 DO 30 I=1,IP SX = 1.0 - ABS(HX(I)) SY = 1.0 - ABS(HY(I)) IF (SX.GE.0.0) NO = NO + NCOR*2 + 1 IF (SY.GE.0.0) NO = NO + NCOR*2 + 1 30 CONTINUE LEN = MIN0(LEF/NO*IP,LEIW/2) IF (LEN.LT.IP) IERROR = 4 FACTOR = FLOAT(IP)/FLOAT(LEN) IF (TEST .AND. IERROR.GT.0) WRITE (IOUT,99991) FACTOR NC = NCOR DO 40 I=1,IP K = IC(I) L = JC(I) IF (K.LE.NC .OR. K.GT.N-NC .OR. L.LE.NC .OR. L.GT.N-NC) IERROR * = 5 40 CONTINUE IF (.NOT.TEST) GO TO 60 C C CHECK THE INPUT VALUES BY PRINTING THEM OUT C WRITE (IOUT,99999) WRITE (IOUT,99998) NC = NCOR DO 50 I=1,IP WRITE (IOUT,99997) I, IC(I), JC(I), HX(I), HY(I) IF (IC(I).LE.NC .OR. IC(I).GT.N-NC .OR. JC(I).LE.NC .OR. * JC(I).GT.N-NC) WRITE (IOUT,99996) I, IC(I), JC(I) 50 CONTINUE C 60 CONTINUE IF (NCOR.EQ.0) GO TO 140 KK = 2*NCOR + 1 C C CHECK WHETHER INSIDE THE REGION ON EACH MESH LINE C THERE ARE AT LEAST 2*NCOR+1 MESH POINTS. C DO 80 J=1,M DO 70 I=1,N F(I,J) = 0.0 70 CONTINUE 80 CONTINUE JMIN = 1000 JMAX = 0 IMIN = 1000 IMAX = 0 DO 90 I=1,IP K = IC(I) L = JC(I) IF (L.LT.JMIN) JMIN = L IF (L.GT.JMAX) JMAX = L IF (K.LT.IMIN) IMIN = K IF (K.GT.IMAX) IMAX = K F(K,L) = 1.0 90 CONTINUE DO 110 I=IMIN,IMAX IFLAG = 0 DO 100 J=JMIN,JMAX IF (F(I,J).EQ.0.0) GO TO 100 IFLAG = IFLAG + 1 IF (IFLAG.EQ.1) L = J IF (IFLAG.GT.1) K = J 100 CONTINUE IF (K-L.LT.KK) IERROR = 6 IF (TEST .AND. K-L.LT.KK) WRITE (6,99995) I, L, K 110 CONTINUE DO 130 J=JMIN,JMAX IFLAG = 0 DO 120 I=IMIN,IMAX IF (F(I,J).EQ.0.0) GO TO 120 IFLAG = IFLAG + 1 IF (IFLAG.EQ.1) L = I IF (IFLAG.GT.1) K = I 120 CONTINUE IF (K-L.LT.KK) IERROR = 6 IF (TEST .AND. K-L.LT.KK) WRITE (6,99995) J, L, K 130 CONTINUE CONTINUE 140 CONTINUE C IF (IERROR.GT.0) RETURN C DO 160 J=1,M DO 150 I=1,N F(I,J) = 0.0 150 CONTINUE 160 CONTINUE 170 CONTINUE I2 = IP + 1 LI = I2 + N - 1 J2 = 1 + NO J3 = J2 + IP J4 = J3 + IP J5 = J4 + IP J6 = J5 + IP J7 = J6 + M J8 = J7 + M J9 = J8 + N J10 = J9 + N J11 = J10 + N/4 J12 = J11 + N*M LW = J12 + N*M - 1 IF (TEST .AND. INTL.EQ.0) WRITE (IOUT,99992) LW, LI, LEW, LEIW IF (IERROR.GT.0) RETURN C C CALL THE MAIN SUBROUTINE C KN = 2*NCOR + 3 CALL HELSIX(INTL, F, NDIM, N, M, H, IC, JC, IW(1), IW(I2), HX, * HY, BX, BY, W(1), W(J2), W(J3), W(J4), W(J5), W(J6), W(J7), * W(J8), W(J9), W(J10), W(J11), W(J12), C, IP, KN, TEST, IERROR) C IF (IERROR.EQ.7) WRITE (IOUT,99994) C RETURN 99999 FORMAT (/1X, 24H THIS IS A TEST PRINTOUT/1X, 17H THE COORDINATES , * 24HOF THE IRREGULAR POINTS,/1X, 30H THE DISTANCES TO THE BOUNDARY * /1X, 41H FOR ALL IRREGULAR POINTS ARE PRINTED OUT/) 99998 FORMAT (/40H NO I J HX HY /) 99997 FORMAT (3I5, 4F10.5) 99996 FORMAT (/3I5, 21H BOUNDARIES TOO CLOSE) 99995 FORMAT (/1X, 38H INSUFFICIENT AMOUNT OF POINTS ON LINE, * I5/8H BETWEEN, I5, 4H AND, I5/27H REFINE OR RESCALE THE MESH/) 99994 FORMAT (/1X, 31H CAPACITANCE MATRIX IS SINGULAR/) 99993 FORMAT (/1X, 4H IP=, I5, 7H IPMAX=, I5/) 99992 FORMAT (/48H THE ACTUAL NEED FOR THE WORKSPACE W AND IW WAS, * 2I5/48H VS THE ESTIMATED SPACE RESERVED AS LEW AND LEIW, 2I5/) 99991 FORMAT (1X, 50H INSUFFICIENT SPACE RESERVED FOR WORKSPACE W OR IW/ * 50H INCREASE THE DIMENSIONS LEW AND LEIW BY A FACTOR, 7H OF AT , * 5HLEAST, F5.2) 99990 FORMAT (1X, 3H N=, I3, 4H M=, I3, 11H AND NDIM=, I5) END C HEL 10 SUBROUTINE HELSIX(INTL, F, NDIM, N, M, H, IC, JC, IPS, IB, HXX, HEL 20 * HYY, FX, FY, AA, SCAL, X, Y, Z, RX, AIX, SAJN, COSAJN, SI, FF, * G, C, IP, KN, TEST, IERROR) DIMENSION C(IP,IP) DIMENSION F(NDIM,M), FF(N,M), G(N,M) DIMENSION IC(IP), JC(IP), HXX(IP), HYY(IP), AA(1), SCAL(IP), * FX(IP), FY(IP), X(IP), Y(IP), Z(IP), IPS(IP), SAJN(N), * COSAJN(N), IB(N), SI(1), RX(M), AIX(M) LOGICAL TEST DIMENSION WE(8,6), CE(10), BB(8), A(10) DIMENSION ENT(1), KC(1), LC(1), IIS(1), IIT(1), DUM(1) C C C PARAMETERS C ARE DESCRIBED BELOW SO FAR AS THEY ARE NOT INCLUDED OR DIFFER C FROM THOSE OF THE SUBROUTINE CMMSIX. C C C FF,G C REAL ARRAYS OF ORDER N BY M. C USED TO STORE THE DEFERRED CORRECTIONS AND SOME INTERMEDIATE C RESULTS IN SUBROUTINE DEFCOR. C C X,Y,Z,SCAL C REAL VECTORS OF LENGTH IP. C Z CONTAINS SCALE FACTORS (PIVOT ELEMENTS) FROM DECOMPOSING C. C SCAL CONTAINS SCALE FACTORS FROM THE DISCRETE LAPLACIAN C MODIFIED AT THE BOUNDARY. C C AA C REAL VECTOR OF LENGTH AT MOST 2*(KN-2)*IP. C CONTAINS COEFFICIENTS FROM APPROXIMATION AT THE BOUNDARY. C C IPS C INTEGER VECTOR OF LENGTH IP. C CONTAINS ROW NUMBERS FROM PIVOTING. C C WE C REAL ARRAY 10 BY 6, CONTAINS WEIGHTS FOR DEFERRED CORRECTIONS. C C ADDITIONAL SPACE USED IN SUBROUTINE STRIP C C SAJN,COSAJN C REAL VECTORS OF LENGTH N. C C IB C INTEGER VECTOR OF LENGTH N. C C RX,AIX C REAL VECTORS OF LENGTH M C C SI C REAL VECTOR OF LENGTH N/4. C C C C,Z,IPS,IC,JC,HXX,HYY,AA,SCAL,WE C CONTAIN INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF HELSIX C WILL BE CALLED WITH INTL>0. C C C THE TOTAL STORAGE FOR ARRAYS IS ABOUT C IP(IP+18) + 2M(N+5) C C ********************* C LIST OF SUBROUTINES : C ********************* C C BDAPRX - MAPPING THE BOUNDARY STENCIL INTO THE IRREGULAR POINTS. C C RHSSIX - COMPUTES RIGHT HAND SIDE TO THE CAPACITANCE SYSTEM. C C CAPSIX - COMPUTES THE CAPACITANCE MATRIX. C C ALFAS - COMPUTES LAGRANGE INTERPOLATION COEFFICIENTS. C C COEGEN - GENERATES COEFFICIENTS FOR THE WEIGHT FACTORS. C C DEFCOR - COMPUTES DEFERRED CORRECTIONS. C C STRIP - FAST HELMHOLTZ SOLVER. C C RFORT,FORT - FAST FOURIER TRANSFORMATION ROUTINES. C C DECOMP,SOLVE - GAUSSIAN ELIMINATION ROUTINES. C C CH2 = 0.0 IF (INTL.GT.0) GO TO 80 C C ************************************************************ C INTL=0. PREPROCESSING PART IN WHICH THE CAPACITANCE MATRIX C IS GENERATED AND FACTORED. C ************************************************************ C C C COMPUTE AND FACTOR THE CAPACITANCE MATRIX C CALL BDAPRX(F, NDIM, N, M, IC, JC, ENT, HXX, HYY, IP, AA, SCAL, * KN, .TRUE.) ENT(1) = 1.0 KC(1) = 2 LC(1) = 2 DO 20 J=1,M DO 10 I=1,N F(I,J) = 0.0 10 CONTINUE 20 CONTINUE CALL STRIP(ENT, KC, LC, 1, DUM, KC, LC, 1, .TRUE., F, NDIM, N, M, * CH2, 3, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) CALL CAPSIX(C, IP, IP, F, NDIM, N, AA, IC, JC, SCAL, HXX, HYY, Z, * KN) CALL DECOMP(C, IP, IP, Z, IPS, IFLAG) IF (IFLAG.EQ.1) IERROR = 7 C C GENERATE THE WEIGHT FACTORS WE. C DO 30 I=1,8 A(I) = 0. BB(I) = 0. 30 CONTINUE BB(1) = 1. A(5) = 2. A(7) = 2. A(9) = 2. DO 50 I=9,12 J = I CALL COEGEN(8, J, CE, BB) J1 = I - 8 DO 40 I1=1,8 WE(I1,J1) = CE(I1) 40 CONTINUE 50 CONTINUE I2 = 3 DO 70 I=3,5,2 I1 = I + 2 CALL COEGEN(I1, I2, CE, A) I2 = I2 + 1 DO 60 J=1,I1 WE(J,I2+1) = CE(J) 60 CONTINUE 70 CONTINUE RETURN C 80 CONTINUE C C ***************************************************** C INTL>0 : THE PREPROCESSING IS NOT REQUIRED, I.E. THE C CAPACITANCE MATRIX IS ALREADY GENERATED AND FACTORED. C ***************************************************** C HSQ = H*H DO 100 J=1,M DO 90 I=1,N G(I,J) = 0.0 FF(I,J) = F(I,J)*HSQ 90 CONTINUE 100 CONTINUE KK = (KN-1)/2 DO 170 KA=1,KK C C SOLVE THE SECOND ORDER SCHEME AND (KA-1) DEFERRED CORRECTIONS. C C C COMPUTE THE RIGHT HAND SIDE TO THE CAPACITANCE MATRIX SYSTEM. C DO 120 J=1,M DO 110 I=1,N F(I,J) = FF(I,J) - G(I,J) 110 CONTINUE 120 CONTINUE CALL RHSSIX(F, NDIM, N, M, FX, FY, Y, IC, JC, HXX, HYY, IP, KN, * .FALSE.) C WRITE(21,11) (Y(I),I=1,IP) C 11 FORMAT (10F8.4) C DO 111 J=1,M C111 WRITE(21,21) (F(I,J),I=1,N) C21 FORMAT (16F7.3) CALL STRIP(ENT, KC, LC, 1, ENT, KC, LC, 1, .FALSE., F, NDIM, N, * M, CH2, 4, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) C DO 211 J=1,M C211 WRITE(21,21) (F(I,J),I=1,N) CALL BDAPRX(F, NDIM, N, M, IC, JC, X, HXX, HYY, IP, AA, SCAL, * KN, .FALSE.) C WRITE(21,11) (X(I),I=1,IP) DO 130 I=1,IP Y(I) = Y(I) - X(I) 130 CONTINUE C WRITE(21,11) (Y(I),I=1,IP) C C COMPUTE THE CHARGE DISTRIBUTION BY BACKSOLVING THE CAPACITANCE C MATRIX SYSTEM. C CALL SOLVE(C, IP, IP, X, Y, Z, IPS) C WRITE(21,11) (X(I),I=1,IP) DO 150 J=1,M DO 140 I=1,N F(I,J) = FF(I,J) - G(I,J) 140 CONTINUE 150 CONTINUE DO 160 I=1,IP K = IC(I) L = JC(I) F(K,L) = F(K,L) + X(I) 160 CONTINUE C C COMPUTE THE SOLUTION AFTER (KA-1) CORRECTIONS AND STORE IN F. C CALL STRIP(ENT, KC, LC, 1, ENT, KC, LC, 1, .FALSE., F, NDIM, N, * M, CH2, 4, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, .FALSE.) IF (KA.EQ.KK) GO TO 170 KB = KA CALL DEFCOR(F, G, NDIM, N, IC, JC, HXX, HYY, WE, KB, IP) 170 CONTINUE C RETURN END C BDA 10 SUBROUTINE BDAPRX(F, NDIM, N, M, IC, JC, OUT, HXX, HYY, IP, A, BDA 20 * SCAL, KN, SAVE) DIMENSION F(NDIM,M), IC(IP), JC(IP), OUT(IP), HXX(IP), HYY(IP), * SCAL(IP), A(1) LOGICAL SAVE DIMENSION ALF(7) C C IF SAVE=.FALSE. C T C THIS SUBROUTINE COMPUTES OUT = (U A) F C T C U A REPRESENTS THE RESTRICTION OF THE DISCRETE LAPLACIAN C TO THE IRREGULAR MESH POINTS USING INTERPOLATION STENCILS C AT THE BOUNDARY OF THE REGION, I.E. INFORMATION AT THE C IRREGULAR POINT AND ITS (KN-2) INTERNAL NEIGHBORS IN C EACH COORDINATE DIRECTION. THE EQUATIONS ARE SCALED SO C THAT THE MAIN DIAGONAL ELEMENTS ARE 1.0. C C IF SAVE=.TRUE. C THE COEFFICIENTS FROM THE STENCILS AT THE BOUNDARY ARE C STORED AND USED LATER IN GENERATING THE CAPACITANCE C MATRIX. C ICOUNT = 0 DO 90 I=1,IP K = IC(I) L = JC(I) II = SIGN(1.0,HXX(I)) JJ = SIGN(1.0,HYY(I)) HX = ABS(HXX(I)) HY = ABS(HYY(I)) SX = 1.0 - HX SY = 1.0 - HY C C X INCREMENTS C R = 0.0 IF (SX.GE.0.0) GO TO 10 DIAG = 2.0 IF (.NOT.SAVE) R = -F(K-1,L) - F(K+1,L) GO TO 40 C C FIND LAGRANGE INTERPOLATION COEFFICIENTS AT THE BOUNDARY. C 10 CALL ALFAS(KN, ALF, SX) DO 20 IK=3,KN ALF(IK) = ALF(IK)/ALF(1) 20 CONTINUE ALF(3) = ALF(3) - 1.0 DO 30 IK=3,KN ICOUNT = ICOUNT + 1 IF (SAVE) A(ICOUNT) = ALF(IK) KK = K - II*(IK-2) IF (.NOT.SAVE) R = F(KK,L)*ALF(IK) + R 30 CONTINUE DIAG = 2.0 + ALF(2)/ALF(1) 40 CONTINUE C C Y INCREMENTS C IF (SY.GE.0.0) GO TO 50 DIAG = 2.0 + DIAG IF (.NOT.SAVE) R = -F(K,L-1) - F(K,L+1) + R GO TO 80 C C FIND LAGRANGE INTERPOLATION COEFFICIENTS AT THE BOUNDARY. C 50 CALL ALFAS(KN, ALF, SY) DO 60 IK=3,KN ALF(IK) = ALF(IK)/ALF(1) 60 CONTINUE ALF(3) = ALF(3) - 1.0 DO 70 IK=3,KN ICOUNT = ICOUNT + 1 IF (SAVE) A(ICOUNT) = ALF(IK) LL = L - JJ*(IK-2) IF (.NOT.SAVE) R = F(K,LL)*ALF(IK) + R 70 CONTINUE DIAG = 2.0 + ALF(2)/ALF(1) + DIAG 80 CONTINUE C C COMPUTE THE ELEMENT OF OUT C SC = 1.0/DIAG IF (.NOT.SAVE) OUT(I) = F(K,L) + R*SC C C SAVE THE COEFFICIENTS FOR GENERATING THE CAPACITANCE MATRIX. C IF (SAVE) SCAL(I) = SC 90 CONTINUE RETURN END C RHS 10 SUBROUTINE RHSSIX(F, NDIM, N, M, RX, RY, RHS, IC, JC, HXX, HYY, RHS 20 * IP, KN, LAPL) DIMENSION F(NDIM,M), IC(IP), JC(IP), HXX(IP), HYY(IP), RX(IP), * RY(IP), RHS(IP) LOGICAL LAPL DIMENSION ALF(7) C C THIS SUBROUTINE COMPUTES THE RIGHT HAND SIDE TO THE C CAPACITANCE MATRIX SYSTEM T C RHS = U F C USING BOUNDARY DATA STORED IN RX AND RY AND APPLYING INTERPOLATION C STENCILS AT THE IRREGULAR POINTS. C SPOT = 0.0 DO 50 I=1,IP K = IC(I) L = JC(I) HX = ABS(HXX(I)) HY = ABS(HYY(I)) SX = 1.0 - HX SY = 1.0 - HY C C X INCREMENTS C S = 0.0 IF (SX.GE.0.0) GO TO 10 DIAG = 2.0 GO TO 20 C C FIND LAGRANGE INTERPOLATION COEFFICIENTS AT THE BOUNDARY. C 10 CALL ALFAS(KN, ALF, SX) DIAG = 2.0 + ALF(2)/ALF(1) S = 1.0/ALF(1) 20 R = S*RX(I) C C Y INCREMENTS C S = 0.0 IF (SY.GE.0.0) GO TO 30 DIAG = 2.0 + DIAG GO TO 40 C C FIND LAGRANGE INTERPOLATION COEFFICIENTS AT THE BOUNDARY. C 30 CALL ALFAS(KN, ALF, SY) DIAG = 2.0 + ALF(2)/ALF(1) + DIAG S = 1.0/ALF(1) 40 R = S*RY(I) + R C C COMPUTE THE ELEMENT OF RHS C SC = 1.0/DIAG IF (.NOT.LAPL) SPOT = F(K,L) RHS(I) = (R+SPOT)*SC 50 CONTINUE RETURN END SUBROUTINE CAPSIX(C, NDIMC, IP, F, NDIMF, N, A, IC, JC, SCAL, CAP 10 * HXX, HYY, P, KN) DIMENSION C(NDIMC,IP), F(NDIMF,1), A(1) DIMENSION IC(IP), JC(IP), SCAL(IP), HXX(IP), HYY(IP), P(IP) C C *** THIS SUBROUTINE COMPUTES THE CAPACITANCE MATRIX C EXPLOITING C *** THE TRANSLATION INVARIANCY PROPERTY OF THE SOLUTION OF HELMHOLTZ C *** EQUATION ON AN INFINITE STRIP WITH PERIODIC BOUNDARY CONDITIONS. C C F REPRESENTS A DISCRETE FUNDAMENTAL SOLUTION, I.E. SOLUTION C TO THE PROBLEM WITH A POINT CHARGE AT (I0,J0). C KN = 3,5 OR 7, I.E. (KN-3)/2 IS THE NUMBER OF DEFERRED CORRECTIONS. C I0 = 2 J0 = 2 ICOUNT = 0 C C SINGLE LAYER APPROACH C DO 130 I=1,IP K = IC(I) L = JC(I) SC = SCAL(I) II = SIGN(1.0,HXX(I)) JJ = SIGN(1.0,HYY(I)) HX = ABS(HXX(I)) HY = ABS(HYY(I)) DO 10 M=1,IP P(M) = 0.0 10 CONTINUE C C X INCREMENTS C IF (HX.LE.1.0) GO TO 30 DO 20 M=1,IP KK = IC(M) LL = JC(M) K2 = IABS(K-1-KK) K3 = IABS(K+1-KK) M2 = MIN0(K2,N-K2) + I0 M3 = MIN0(K3,N-K3) + I0 L1 = IABS(L-LL) + J0 P(M) = -F(M2,L1) - F(M3,L1) 20 CONTINUE GO TO 60 C C BOUNDARY AT A DISTANCE LESS THAN OR EQUAL TO H FROM THE C IRREGULAR POINT IN THIS DIRECTION. C 30 DO 50 IK=3,KN ICOUNT = ICOUNT + 1 AX = A(ICOUNT) KI = K - II*(IK-2) DO 40 M=1,IP KK = IC(M) LL = JC(M) KX = IABS(KI-KK) MX = MIN0(KX,N-KX) + I0 LY = IABS(L-LL) + J0 P(M) = AX*F(MX,LY) + P(M) 40 CONTINUE 50 CONTINUE 60 CONTINUE C C Y INCREMENTS C IF (HY.LE.1.0) GO TO 80 DO 70 M=1,IP KK = IC(M) LL = JC(M) K1 = IABS(K-KK) M1 = MIN0(K1,N-K1) + I0 L2 = IABS(L-1-LL) + J0 L3 = IABS(L+1-LL) + J0 P(M) = -F(M1,L2) - F(M1,L3) + P(M) 70 CONTINUE GO TO 110 C C BOUNDARY AT A DISTANCE LESS THAN OR EQUAL TO H FROM THE C IRREGULAR POINT IN THIS DIRECTION. C 80 DO 100 IK=3,KN ICOUNT = ICOUNT + 1 AY = A(ICOUNT) LJ = L - JJ*(IK-2) DO 90 M=1,IP KK = IC(M) LL = JC(M) LY = IABS(LJ-LL) + J0 KX = IABS(K-KK) MX = MIN0(KX,N-KX) + I0 P(M) = AY*F(MX,LY) + P(M) 90 CONTINUE 100 CONTINUE 110 CONTINUE C C COMPUTE THE ELEMENT OF MATRIX C. C DO 120 M=1,IP KK = IC(M) LL = JC(M) KX = IABS(K-KK) MX = MIN0(KX,N-KX) + I0 LY = IABS(L-LL) + J0 C(I,M) = F(MX,LY) + P(M)*SC 120 CONTINUE 130 CONTINUE RETURN END C ALF 10 SUBROUTINE ALFAS(N, ALF, S) ALF 20 DIMENSION ALF(N) C C *** AFTER V.PEREYRA, STANFORD REPORT CS-348(1973),P.46 C *** ADOPTED FOR THIS PARTICULAR CHOICE OF V & RHS. C *** ALF - SOLUTION OF A SYSTEM OF EQNS. OF ORDER N >=3. C *** WITH A VANDERMONDE MATRIX V(0,1,...,N-1) C *** AND RHS(I)=S**(I-1), I=1,...,N;0<=S<=1. C ALF(1) = 1 DO 10 I=2,N ALF(I) = S*ALF(I-1) 10 CONTINUE NN = N - 1 DO 30 I=2,NN I1 = I - 1 LL = N - I DO 20 J=1,LL K = N - J + 1 ALF(K) = ALF(K) - FLOAT(I1)*ALF(K-1) 20 CONTINUE 30 CONTINUE DO 50 I=1,NN K = N - I DK = K XKIN = 1D0/DK KP1 = K + 1 DO 40 J=KP1,N ALF(J) = ALF(J)*XKIN JM1 = J - 1 ALF(JM1) = ALF(JM1) - ALF(J) 40 CONTINUE 50 CONTINUE RETURN END C COE 10 SUBROUTINE COEGEN(N, NP, CC, BB) COE 20 DIMENSION CC(1), BB(1), ALF(10) C C THIS SUBROUTINE GENERATES COEFFICIENTS C USED IN COMPUTING C THE WEIGHT FACTORS FOR THE DEFERRED CORRECTIONS. DO 10 I=1,N CC(I) = BB(I) 10 CONTINUE DO 20 I=1,N ALF(I) = I - NP 20 CONTINUE NN = N - 1 N1 = N + 1 DO 40 I=1,NN LL = N - I DO 30 J=1,LL K = N1 - J CC(K) = CC(K) - ALF(I)*CC(K-1) 30 CONTINUE 40 CONTINUE DO 60 I=1,NN K = N - I KM1 = K + 1 DO 50 J=KM1,N JK = J - K CC(J) = CC(J)/(ALF(J)-ALF(JK)) JM1 = J - 1 CC(JM1) = CC(JM1) - CC(J) 50 CONTINUE 60 CONTINUE RETURN END C DEF 10 SUBROUTINE DEFCOR(F, G, NDIM, N, IC, JC, HXX, HYY, WE, KA, IP) DEF 20 DIMENSION F(NDIM,1), G(N,1), WE(8,6) DIMENSION IC(IP), JC(IP), HXX(IP), HYY(IP) C C THIS SUBROUTINE COMPUTES DEFERRED CORRECTIONS, STORED IN G, N BY M. C C ON INPUT C C KA = 1 OR 2, THE NUMBER OF THE DEFERRED CORRECTION. C F - N BY M ARRAY CONTAINING THE SOLUTION. C WE - CONTAINS THE WEIGHT FACTORS. C IC,JC - COORDINATES OF IRREGULAR POINTS. C HXX,HYY - SIGNED DISTANCES TO THE BOUNDARY. C IP - NUMBER OF IRREGULAR POINTS C C HORIZONTAL (Y) EXTRAPOLATION AND CORRECTION C N4 = N - 3 K3 = KA + 1 DO 60 I=1,IP IF (ABS(HYY(I)).GT.1.0) GO TO 60 IX = IC(I) JY = JC(I) IF (HYY(I).GT.0.0) GO TO 30 C IRREGULAR ON THE LEFT J1 = JY + 8 DO 20 L1=1,K3 TE = 0. DO 10 I1=1,8 J = J1 - I1 TE = TE + F(IX,J)*WE(I1,L1) 10 CONTINUE L = JY - L1 F(IX,L) = TE 20 CONTINUE GO TO 60 C IRREGULAR ON THE RIGHT 30 J1 = JY - 8 DO 50 L1=1,K3 TE = 0. DO 40 I1=1,8 J = J1 + I1 TE = TE + F(IX,J)*WE(I1,L1) 40 CONTINUE L = JY + L1 F(IX,L) = TE 50 CONTINUE 60 CONTINUE C C HORIZONTAL CORRECTION C K21 = KA + 2 I2 = KA + 4 I21 = 2*KA + 3 DO 90 I=4,N4 DO 80 J=4,N4 I10 = J - K21 TE = 0. DO 70 I1=1,I21 L = I10 + I1 TE = TE + F(I,L)*WE(I1,I2) 70 CONTINUE G(I,J) = TE 80 CONTINUE 90 CONTINUE C C VERTICAL (X) EXTRAPOLATION AND CORRECTION C DO 150 I=1,IP IF (ABS(HXX(I)).GT.1.0) GO TO 150 IX = IC(I) JY = JC(I) IF (HXX(I).LT.0.0) GO TO 120 C IRREGULAR ABOVE J1 = IX - 8 DO 110 L1=1,K3 TE = 0. DO 100 I1=1,8 J = J1 + I1 TE = TE + F(J,JY)*WE(I1,L1) 100 CONTINUE L = IX + L1 F(L,JY) = TE 110 CONTINUE GO TO 150 C IRREGULAR BELOW 120 J1 = IX + 8 DO 140 L1=1,K3 TE = 0. DO 130 I1=1,8 J = J1 - I1 TE = TE + F(J,JY)*WE(I1,L1) 130 CONTINUE L = IX - L1 F(L,JY) = TE 140 CONTINUE 150 CONTINUE C C VERTICAL CORRECTION C DO 180 I=4,N4 DO 170 J=4,N4 I10 = I - K21 TE = 0. DO 160 I1=1,I21 K = I10 + I1 TE = TE + F(K,J)*WE(I1,I2) 160 CONTINUE G(I,J) = G(I,J) + TE 170 CONTINUE 180 CONTINUE RETURN END C DOM 10 C DOM 20 C THESE ARE THE USER SUPPLIED SUBROUTINES IN WHICH INFORMATION DOM 30 C ABOUT THE REGION AND DATA AT THE BOUNDARY ARE GIVEN. DOM 40 C DOM 50 C THESE SUBROUTINES COMPUTE: DOM 60 C DOM 70 C 1. IP - THE NUMBER OF IRREGULAR MESH POINTS, I.E., POINTS DOM 80 C STRICTLY INSIDE THE REGION WHICH HAVE AT LEAST ONE DOM 90 C NEIGHBOR OUTSIDE OR ON THE BOUNDARY. DOM 100 C IP SERVES AS THE LENGTH OF VECTORS IC,JC. DOM 110 C 2. IPP - THE LENGTH OF VECTORS HX,HY,TX,TY,BX,BY. DOM 120 C 3. IC,JC - INTEGER VECTORS OF LENGTH IP CONTAINING DOM 130 C COORDINATES (I,J) OF ALL IRREGULAR POINTS ENUMERATED DOM 140 C IN AN ARBITRARY ORDER. DOM 150 C THE VALUES OF I MUST BE LARGER THAN 1 AND LESS THAN N. DOM 160 C THE VALUES OF J MUST BE LARGER THAN 1 AND LESS THAN M. DOM 170 C 4. HX,HY - REAL VECTORS OF LENGTH IPP WITH SIGNED SHORTEST DOM 180 C DISTANCES IN BOTH COORDINATE DIRECTIONS FROM THE BOUNDARY DOM 190 C TO AN IRREGULAR POINT (I,J) IN UNITS OF THE STEP SIZE H. DOM 200 C IF (I,J) HAS A NEIGHBOR OUTSIDE THE REGION THEN THE VALUES OF DOM 210 C HX OR HY MUST BELONG TO THE INTERVAL -1.0 TO +1.0 (EXCLUDING 0.0DOM 220 C OTHERWISE THE PROPER VALUE OF MAGNITUDE LARGER THAN 1.0 DOM 230 C MUST BE PUT IN. DOM 240 C 5. TX,TY - REAL VECTORS OF LENGTH IPP WITH ABSOLUTE VALUES OF DOM 250 C TANGENTS OF THE ANGLE BETWEEN THE NORMAL AND THE X-AXIS DOM 260 C AT (I+HX,J) ( AND THE Y-AXIS AT (I,J+HY) ). DOM 270 C TX,TY ARE USED ONLY FOR A PROBLEM WITH THE NEUMANN BOUNDARY DOM 280 C CONDITIONS. THE VALUES OF TX,TY MUST BE POSITIVE. DOM 290 C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0, A DUMMY DOM 300 C VALUE OF ARBITRARY MAGNITUDE CAN BE PUT IN. DOM 310 C 6. BX,BY - REAL VECTORS OF LENGTH IPP WITH BOUNDARY VALUES DOM 320 C AT (IC(I)+HX,JC(I)) AND (IC(I),JC(I)+HY). DOM 330 C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0 THEN DOM 340 C A DUMMY VALUE OF BX OR BY SHOULD BE PUT IN. DOM 350 C 7. F - REAL ARRAY OF ORDER NDIM BY M CONTAINING THE VALUES OF DOM 360 C THE RIGHT HAND SIDE OF THE HELMHOLTZ EQUATION FOR THE DOM 370 C GRID POINTS INSIDE THE IRREGULAR REGION AND SUPPLEMENTED DOM 380 C BY ARBITRARY VALUES, FOR INSTANCE BY ZERO VALUES, FOR DOM 390 C THE GRID POINTS OUTSIDE THE IRREGULAR REGION. DOM 400 C 8. H - REAL VARIABLE, THE MESH LENGTH IN BOTH COORDINATE DIRECTIONS.DOM 410 C THE DEFAULT VALUE IS H=(B-A)/N. DOM 420 C DOM 430 C THE SET OF IRREGULAR POINTS CONSISTS OF THE SET (1) OF IP1 DOM 440 C POINTS WITH AT MOST 2 NEIGHBORS OUTSIDE THE REGION (OR ON DOM 450 C THE BOUNDARY) AND THE REMAINING SET (2) OF IP2 POINTS WITH DOM 460 C EXACTLY 3 NEIGHBORS OUTSIDE THE REGION (OR ON THE BOUNDARY). DOM 470 C THUS IP = IP1 + IP2. DOM 480 C IF A POINT BELONGS TO THE SET (2) THEN THE DISTANCES FROM DOM 490 C THE BOUNDARY MUST BE STORED IN A PAIR OF LOCATIONS: DOM 500 C (HX(I),HX(I+1)) AND (HY(I),HY(I+1)). NOTE THAT THE SECOND ENTRY DOM 510 C MUST BE LARGER IN MAGNITUDE THAN THE FIRST. A SIMILAR SITUATION DOM 520 C OCCURS FOR THE BOUNDARY VALUES BX,BY AND THE TANGENTS TX,TY. DOM 530 C THUS THE LENGTH OF THESE VECTORS MUST BE IPP = IP1 + IP2*2. DOM 540 C DOM 550 C DOM 560 SUBROUTINE DOMAIN(N, M, IC, JC, HX, HY, TX, TY, IPMAX, IP, IPP, DOM 570 * H, DIR) REAL TX(IPMAX), TY(IPMAX), HX(IPMAX), HY(IPMAX) INTEGER IC(IPMAX), JC(IPMAX) LOGICAL DIR COMMON /INPUT/ A, B, C, D, XC, YC, R, NN C C THIS SUBROTINE COMPUTES THE MESH LENGTH H, NUMBER OF IRREGULAR C POINTS IP, LENGTH OF DATA VECTORS IPP, COORDINATES OF IRREGULAR C POINTS STORED IN IC,JC, DISTANCES FROM THE BOUNDARY STORED IN HX C AND TANGENTS OF THE NORMALS (FOR THE NEUMANN PROBLEM) STORED IN C C THE ARGUMENTS ARE DEFINED AS: C C C * * * * * * * * * * * * ON INPUT * * * * * * * * * * * * * * * * * C C C N C C THE NUMBER OF MESH LENGTHS IN THE X DIRECTION. C N MUST BE A POWER OF 2 AND NOT LESS THAN 8. C C M C C THE NUMBER OF MESH LENGTHS IN THE Y DIRECTION. C M SHOULD APPROXIMATELY BE M = N * (D-C)/(B-A). C WHERE A,B IS THE RANGE OF X IN THE RECTANGULAR REGION, C I.E., A =< X =< B, AND C,D IS THE RANGE OF Y C IN THE RECTANGULAR REGION, I.E., C =< Y =< D. C M MUST BE GREATER THAN OR EQUAL TO N. C C IPMAX C C AN UPPER ESTIMATE OF IP, THE NUMBER OF IRREGULAR MESH POINTS. C SERVES AS THE DIMENSION OF THE VECTORS IC,JC,HX,HY,TX,TY,BX,BY. C C DIR C C LOGICAL VARIABLE, EQUAL TO .TRUE. FOR A DIRICHLET PROBLEM, C AND EQUAL TO .FALSE. FOR A NEUMANN PROBLEM. C C C * * * * * * * * * * * * ON OUTPUT * * * * * * * * * * * * * * * * C C C IC,JC C C INTEGER VECTORS, EXPLAINED ABOVE. C C HX,HY,TX,TY C C REAL VECTORS, EXPLAINED ABOVE. C C H C C REAL VARIABLE, EXPLAINED ABOVE. C C IPMAX,IP,IPP C C INTEGER VARIABLES, EXPLAINED ABOVE. C C * * * C C BELOW IS A TEST EXAMPLE FOR A CIRCULAR REGION C C 2 2 2 C ( X - XC ) + ( Y - YC ) = R C C THIS REGION IS IMBEDDED IN A RECTANGLE (HERE A SQUARE) C WITH SIDES SLIGHTY LARGER THAN R; SEE BELOW ON R. C C THE EXACT SOLUTION IS: U(X,Y) = X**5 + Y**5. C C IN ADDITION, IN THE TEST EXAMPLE THE FOLLOWING PARAMETERS C ARE USED TO DESCRIBE THE REGION AND THE SOLUTION: C C A,B C C THE RANGE OF X IN THE RECTANGULAR REGION, I.E., C A =< X =< B. A MUST BE LESS THAN B. C C C,D C C THE RANGE OF Y IN THE RECTANGULAR REGION, I.E., C C =< Y =< D. C MUST BE LESS THAN D. C (B-A) MUST BE LESS THAN OR EQUAL TO (D-C). C C XC,YC C C X AND Y COORDINATES OF THE CENTER OF THE CIRCULAR REGION. C THE DEFAULT VALUES ARE XC=(B-A)/2 AND YC=(D-C)/2. C C R C C RADIUS OF THE CIRCULAR REGION. C THE DEFAULT VALUE IS R=0.75*XC. C R MUST BE LESS THAN OR EQUAL TO (N-3)/2N. C C NN C C INTEGER VARIABLE, EXPONENT IN THE TEST EXAMPLE U=X**NN+Y**NN; C NN MUST BE POSITIVE ( AND ODD FOR THE NEUMANN PROBLEM). C C C C C THE TEST EXAMPLE BELOW CAN EASILY BE EXTENDED TO ELLIPSOIDAL REGIONS C C C 2 2 2 C THE REGION IS (X-XC) + (Y-YC) =< R C C C INITIALIZE C A = 0.0 B = 1.0 C = 0.0 D = 1.0 XC = 0.5 YC = 0.5 R = 0.38 NN = 5 H = (B-A)/FLOAT(N) C EPS = 1.0 10 EPS = EPS*0.5 E1 = EPS + 1. IF (E1.GT.1.0) GO TO 10 TEST = 8.0*EPS RSQ = R**2 C IP1 = 0 IP2 = 0 C C TEST EACH POINT IN THE RECTANGLE. FIND THE IRREGULAR POINTS, C I.E., THOSE INTERIOR POINTS WHO HAVE AT LEAST ONE NEIGHBOR C OUTSIDE THE REGION. C DO 40 J=1,M YY = FLOAT(J-1)*H + C T2 = (YY-YC)**2 DO 30 I=1,N XX = FLOAT(I-1)*H + A T1 = (XX-XC)**2 IF (T1+T2-RSQ.GT.(-TEST)) GO TO 30 C C POINT (X,Y) IS INSIDE THE REGION. C C TEST WHETHER (X,Y) IS AN IRREGULAR POINT. C IR = 0 XT = SQRT(ABS(RSQ-T2)) XD1 = XT + XC - XX XD2 = -XT + XC - XX C A MESH POINT IS CONSIDERED TO BE AN EXTERIOR ONE C IF A DISTANCE TO THE BOUNDARY IS LESS THAN TEST. IF (ABS(XD1)-H.LT.TEST) IR = IR + 1 IF (ABS(XD2)-H.LT.TEST) IR = IR + 1 YT = SQRT(ABS(RSQ-T1)) YD1 = YT + YC - YY YD2 = -YT + YC - YY IF (ABS(YD1)-H.LT.TEST) IR = IR + 1 IF (ABS(YD2)-H.LT.TEST) IR = IR + 1 IF (IR.EQ.0) GO TO 30 C C POINT (X,Y) IS AN IRREGULAR POINT. CALCULATE AND STORE C THE COORDINATES AND SIGNED SHORTEST DISTANCES TO THE C BOUNDARY IN COORDINATE DIRECTIONS IN UNITS OF H. C XD = XD1 IF (ABS(XD1).GT.ABS(XD2)) XD = XD2 YD = YD1 IF (ABS(YD1).GT.ABS(YD2)) YD = YD2 XDH = XD/H YDH = YD/H IF (IR.GT.2) GO TO 20 C C POINT (X,Y) HAS AT MOST TWO NEIGHBORS OUTSIDE THE REGION. C IP1 = IP1 + 1 HX(IP1) = XDH HY(IP1) = YDH IC(IP1) = I JC(IP1) = J C C CALCULATE AND STORE TANGENTS FOR THE NEUMANN PROBLEM. C IF (DIR) GO TO 30 C C NEUMANN BOUNDARY CONDITIONS C TX(IP1) = ABS((YY-YC)/(XX+XD-XC)) TY(IP1) = ABS((XX-XC)/(YY+YD-YC)) GO TO 30 20 CONTINUE C C POINT (X,Y) HAS THREE NEIGHBORS OUTSIDE THE REGION. C IP2 = IP2 + 1 XXD = XD2 IF (ABS(XD1).GT.ABS(XD2)) XXD = XD1 YYD = YD2 IF (ABS(YD1).GT.ABS(YD2)) YYD = YD1 HX2 = XXD/H HY2 = YYD/H C C INITIALLY STORE ALL THE INFORMATION IN THE LAST C LOCATIONS OF THE VECTORS IC,JC,HX,HY,TX & TY. C INC = IPMAX - IP2 + 1 IC(INC) = I JC(INC) = J INH = IPMAX - IP2*2 + 1 HX(INH) = XDH HY(INH) = YDH HX(INH+1) = HX2 HY(INH+1) = HY2 C C CALCULATE AND STORE TANGENTS FOR THE NEUMANN PROBLEM. C IF (DIR) GO TO 30 C C NEUMANN BOUNDARY CONDITIONS C TX(INH) = ABS((YY-YC)/(XX+XD-XC)) TY(INH) = ABS((XX-XC)/(YY+YD-YC)) TX(INH+1) = ABS((YY-YC)/(XX+XXD-XC)) TY(INH+1) = ABS((XX-XC)/(YY+YYD-YC)) 30 CONTINUE 40 CONTINUE C IP = IP1 + IP2 IPP = IP + IP2 IF (IP2.EQ.0) GO TO 60 C C SHIFT BACKWARD THE INFORMATION ABOUT IP2 MESH POINTS, C I.E.,THOSE WITH THREE EXTERIOR NEIGHBORS, C SO THAT THERE ARE NO GAPS IN FILLING ALL THE VECTORS. C DO 50 I=1,IP2 II = IP1 + I JJ = IPMAX - IP2 + I IC(II) = IC(JJ) JC(II) = JC(JJ) KK = IP1 + I*2 - 1 LL = IPMAX - (IP2-I)*2 - 1 HX(KK) = HX(LL) HY(KK) = HY(LL) HX(KK+1) = HX(LL+1) HY(KK+1) = HY(LL+1) IF (DIR) GO TO 50 TX(KK) = TX(LL) TY(KK) = TY(LL) TX(KK+1) = TX(LL+1) TY(KK+1) = TY(LL+1) 50 CONTINUE 60 CONTINUE C RETURN END SUBROUTINE CHARGE(F, NDIM, N, M, H, CON) CHA 10 DIMENSION F(NDIM,M) COMMON /INPUT/ A, B, C, D, XC, YC, R, NN C C THIS SUBROUTINE COMPUTES THE RIGHT SIDE C OF THE HELMHOLTZ EQUATION AND STORES IT IN F. C C *** ON INPUT *** C C F C C REAL ARRAY OF ORDER NDIM BY M; UNINITIALIZED. C C NDIM C C THE FIRST DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING CHARGE THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. NDIM MUST BE AT LEAST N. C C CON C C REAL VARIABLE, THE CONSTANT IN THE HELMHOLTZ EQUATION. C C *** ON OUTPUT *** C C F C REAL ARRAY OF ORDER NDIM BY M CONTAINING THE VALUES OF C THE RIGHT HAND SIDE OF THE HELMHOLTZ EQUATION FOR THE C GRID POINTS INSIDE THE IRREGULAR REGION AND SUPPLEMENTED C BY ARBITRARY VALUES, FOR INSTANCE BY ZERO VALUES, FOR C THE GRID POINTS OUTSIDE THE IRREGULAR REGION. C C VALUES IN THE COMMON BLOCK ARE INITIALIZED C IN THE SUBROUTINE DOMAIN. C EPS = 1.0 10 EPS = EPS*0.5 E1 = EPS + 1. IF (E1.GT.1.0) GO TO 10 TEST = 8.0*EPS RSQ = R*R DNN = NN*(NN-1) NN2 = NN - 2 C DO 30 J=1,M YY = FLOAT(J-1)*H + C T2 = (YY-YC)**2 DO 20 I=1,N C SET RIGHT SIDE F TO ZERO F(I,J) = 0.0 XX = FLOAT(I-1)*H + A T1 = (XX-XC)**2 IF (T1+T2-RSQ.GT.(-TEST)) GO TO 20 C C POINT (X,Y) IS INSIDE THE REGION. C CALCULATE AND STORE THE VALUES OF THE RIGHT SIDE C OF THE HELMHOLTZ EQUATION. C IF (NN.GE.2) F(I,J) = -DNN*(XX**NN2+YY**NN2) + * CON*(XX**NN+YY**NN) IF (NN.LT.2) F(I,J) = CON*(XX**NN+YY**NN) 20 CONTINUE 30 CONTINUE C RETURN END SUBROUTINE BNDRY(IC, JC, HX, HY, BX, BY, IP, IPP, H, DIR) BND 10 DIMENSION IC(IP), JC(IP), HX(IPP), HY(IPP), BX(IPP), BY(IPP) LOGICAL DIR COMMON /INPUT/ A, B, C, D, XC, YC, R, NN C C THIS SUBROUTINE CALCULATES THE BOUNDARY VALUES C AND STORES THEM IN BX AND BY. C C BX,BY C REAL VECTORS OF LENGTH IPP WITH BOUNDARY VALUES C AT (IC(I)+HX,JC(I)) AND (IC(I),JC(I)+HY). C IF HX OR HY ARE LARGER IN MAGNITUDE THAN 1.0 THEN C A DUMMY VALUE OF BX OR BY SHOULD BE PUT IN. C C VALUES IN THE COMMON BLOCK ARE INITIALIZED C IN THE SUBROUTINE DOMAIN. C IP1 = IP - (IPP-IP) IP2 = IP - IP1 DN = NN C DO 40 I=1,IP K = IC(I) L = JC(I) XX = FLOAT(K-1)*H + A YY = FLOAT(L-1)*H + C IF (I.GT.IP1) GO TO 20 XD = HX(I)*H YD = HY(I)*H C IF (DIR) GO TO 10 C C NEUMANN BOUNDARY CONDITIONS C C FOR U=X**NN+Y**NN, NN A POSITIVE INTEGER, C NORMAL DERIVATIVE AT THE BOUNDARY OF A CIRCLE C ( X - XC )**2 + ( Y - YC )**2 = R**2 C IS: DU/DN=-NN*((X-XC)*X**(NN-1)+(Y-YC)*Y**(NN-1))/R. C BX(I) = -DN*((XX+XD-XC)*(XX+XD)**(NN-1)+(YY-YC)*YY**(NN-1))/R BY(I) = -DN*((XX-XC)*XX**(NN-1)+(YY+YD-YC)*(YY+YD)**(NN-1))/R GO TO 40 10 CONTINUE C C DIRICHLET BOUNDARY CONDITIONS C BX(I) = (XX+XD)**NN + YY**NN BY(I) = XX**NN + (YY+YD)**NN GO TO 40 C 20 CONTINUE C C POINTS WITH THREE NEIGHBORS OUTSIDE THE REGION C J = IP1 + (I-IP1)*2 - 1 XD = HX(J)*H YD = HY(J)*H XXD = HX(J+1)*H YYD = HY(J+1)*H C IF (DIR) GO TO 30 C C NEUMANN BOUNDARY CONDITIONS C BX(J) = -DN*((XX+XD-XC)*(XX+XD)**(NN-1)+(YY-YC)*YY**(NN-1))/R BY(J) = -DN*((XX-XC)*XX**(NN-1)+(YY+YD-YC)*(YY+YD)**(NN-1))/R BX(J+1) = -DN*((XX+XXD-XC)*(XX+XD)**(NN-1)+(YY-YC)*YY**(NN-1))/R BY(J+1) = -DN*((XX-XC)*XX**(NN-1)+(YY+YYD-YC)*(YY+YD)**(NN-1))/R GO TO 40 30 CONTINUE C C DIRICHLET BOUNDARY CONDITIONS C BX(J) = (XX+XD)**NN + YY**NN BY(J) = XX**NN + (YY+YD)**NN BX(J+1) = (XX+XXD)**NN + YY**NN BY(J+1) = XX**NN + (YY+YYD)**NN C 40 CONTINUE C RETURN END SUBROUTINE OUTPUT(F, NDIM, N, M, H, DIR, IOUT) OUT 10 DIMENSION F(NDIM,M) LOGICAL DIR COMMON /INPUT/ A, B, C, D, XC, YC, R, NN DIMENSION LG(128) C C THIS SUBROUTINE COMPARES THE NUMERICAL SOLUTION WITH THE EXACT C SOLUTION FOR THE GRID POINTS INSIDE THE IRREGULAR REGION FOR C THE TEST EXAMPLE , AND PRINTS OUT THE NUMBER OF CORRECT DIGITS C FOR ALL SUCH POINTS. C FINALLY THE ERROR IN THE SOLUTION IN MAXIMUM AND L2 NORMS C IS COMPUTED AND PRINTED OUT. C C VALUES IN THE COMMON BLOCK ARE INITIALIZED C IN THE SUBROUTINE DOMAIN. C EPS = 1.0 10 EPS = EPS*0.5 E1 = EPS + 1. IF (E1.GT.1.0) GO TO 10 TEST = 8.0*EPS IF (DIR) GO TO 40 C C FOR THE NEUMANN PROBLEM THE SOLUTION IS UNIQUE UP TO A CONSTANT. C KNOWING THE SOLUTION AT ONE POINT WE SHIFT IT ACCORDINGLY. C K = (XC-A)/H + 1.0000001 L = (YC-C)/H + 1.0000001 C IN THIS EXAMPLE THE SOLUTION IS X**NN+Y**NN AT (XC,YC) S = F(K,L) - (XC**NN+YC**NN) WRITE (IOUT,99999) S DO 30 J=1,M DO 20 I=1,N F(I,J) = F(I,J) - S 20 CONTINUE 30 CONTINUE 40 CONTINUE C RSQ = R**2 T = 0.0 V = 0.0 L = 0 DO 60 J=1,M YY = FLOAT(J-1)*H + C T2 = (YY-YC)**2 DO 50 I=1,N XX = FLOAT(I-1)*H + A T1 = (XX-XC)**2 LG(I) = 0 IF (T1+T2-RSQ.GT.-TEST) GO TO 50 C C POINT (X,Y) IS INSIDE THE REGION. C L = L + 1 C C THE EXACT SOLUTION TO THE TEST EXAMPLE. C INSERT A CODE FOR PRINTING OUT THE SOLUTION, IF REQUIRED. C S = (FLOAT(I-1)*H+A)**NN + (FLOAT(J-1)*H+C)**NN S = ABS(F(I,J)-S) IF (S.GT.V) V = S T = T + S*S LG(I) = -ALOG10(S+1E-20) + 0.5 50 CONTINUE WRITE (IOUT,99998) (LG(I),I=1,N) 60 CONTINUE T = SQRT(T/FLOAT(L)) WRITE (IOUT,99997) L, V, T C ON CDC COMPUTERS ONE CAN USE 64I2.0 FORMAT, INSTEAD. 99999 FORMAT (/1X, 50H SYSTEM IS SINGULAR AND THE SOLUTION IS COMPUTED , * /1X, 33H UP TO THE ADDITIVE CONSTANT C = , F8.4/) 99998 FORMAT (1X, 64I2) 99997 FORMAT (/1X, 38H NUMBER OF POINTS INSIDE THE REGION L=, I6/1X, * 38H MAXIMUM ERROR V=, E15.5/1X, 9H MEAN ERR, * 29HOR IN L2 SENSE T=, E15.5) RETURN END C BLK 10 BLOCK DATA BLK 20 COMMON /INPUT/ A, B, C, D, XC, YC, R, NN BLK 30 END BLK 40 SUBROUTINE STRIP(FS, IS, JS, NOSP, XT, IT, JT, NOTP, FIRST, F, STR 10 * NDIM, N, M, CH2, MODE, RX, AIX, SAJN, COSAJN, IB, SI, IIS, IIT, * SUB) DIMENSION F(NDIM,M) DIMENSION FS(NOSP), IS(NOSP), JS(NOSP), IT(NOTP), JT(NOTP), * XT(NOTP) DIMENSION RX(M), AIX(M), IIS(NOSP), IIT(NOTP) DIMENSION SAJN(N), COSAJN(N), IB(N), SI(1) LOGICAL FIRST, SUBSET, SUB C C THIS SUBROUTINE SOLVES THE HELMHOLTZ EQUATION ON A RECTANGLE C WITH PERIODIC BOUNDARY CONDITIONS IN THE X DIRECTION AND F=0 C OUTSIDE THE RECTANGLE IN Y DIRECTION. C TWO DIFFERENT SOLVERS ARE EMPLOYED HERE - ONE USING FFT C AND THE OTHER WHICH EXPLOITS THE SPARSITY OF THE PROBLEM. C ANALYSIS SYNTHESIS C IF MODE=1 SOURCES ARE IN FS AND TARGETS IN FT: SPARSE SPARSE C IF MODE=2 SOURCES ARE IN F AND TARGETS IN FT: FFT SPARSE C IF MODE=3 SOURCES ARE IN FS AND TARGETS IN F : SPARSE FFT C IF MODE=4 SOURCES ARE IN F AND TARGETS IN F : FFT FFT C C NOTE : MODE=3 REQUIRES PROPER INITIALIZATION OF ARRAY F. C C INPUT PARAMETERS C C FIRST - .TRUE. IF SAJN,COSAJN AND SI ARE PRECOMPUTED. C SUB - .TRUE. IF THE COORDINATES (IS,JS) ARE A SUBSET OF (IT,JT) C OR VICE VERSA (USED ONLY WHEN MODE=1), .FALSE. OTHERWISE. C IF MODE = 2 OR 4 C F - REAL ARRAY OF ORDER N BY M CONTAINING CHARGES(SOURCES). C NDIM - THE FIRST DIMENSION OF THE ARRAY F AS IT APPEARS C IN THE PROGRAM CALLING STRIP. C IF MODE = 1 OR 3 C FS - REAL VECTOR OF LENGTH NOSP CONTAINING SOURCES C IS,JS - INTEGER VECTORS OF LENGTH NOSP CONTAINING COORDINATES C OF THE SOURCES. C C OUTPUT PARAMETERS C C IF MODE = 3 OR 4 C F - REAL ARRAY OF ORDER N BY M CONTAINING THE SOLUTION(TARGETS). C IF MODE = 1 OR 2 C XT - REAL VECTOR OF LENGTH NOTP CONTAINING TARGETS. C IT,JT - INTEGER VECTORS OF LENGTH NOSP CONTAINING COORDINATES C OF THE TARGETS. C C WORKSPACE C C FOR THE FFT SOLVER C IB - INTEGER ARRAY OF LENGTH N. C SI - REAL ARRAY OF LENGTH N/4. C FOR THE SPARSE SOLVER C SAJN,COSAJN - REAL ARRAYS OF LENGTH N C RX,AIX - REAL ARRAYS OF LENGTH M C IIS,IIT - INTEGER ARRAYS OF LENGTH NOSP AND NOTP, RESPECTIVELY. C C FOR REFERENCE SEE : W.PROSKUROWSKI - TRANS. ON NUMERICAL SOFTWARE, C VOL.5 (1979), P.36-49 AND W.PROSKUROWSKI AND O.WIDLUND - MATH. C COMP.,VOL.30 (1976), P.433-466. C SUBSET = SUB FNI = 1./FLOAT(N) PIN = 6.2831853071796*FNI LOG2N = 2 10 LOG2N = LOG2N + 1 IF (N.GT.2**LOG2N) GO TO 10 IF (.NOT.FIRST) GO TO 30 C C PRECOMPUTE SINE TABLE C ANG = 0. DO 20 I=1,N SAJN(I) = SIN(ANG) COSAJN(I) = COS(ANG) ANG = ANG + PIN 20 CONTINUE C CALL RFORT(F, NDIM, LOG2N, SI, 0, M, IB) C 30 CONTINUE C IF (MODE.NE.1) SUBSET = .FALSE. IF (MODE.EQ.2*(MODE/2)) GO TO 420 C C ********************************************* C THIS PART IS A SOLVER THAT EXPLOITS SPARSITY. C ********************************************* C C MODE = 1 OR 3 C IF (SUBSET) NOP = MAX0(NOSP,NOTP) IF (.NOT.SUBSET) NOP = NOSP DO 40 K=1,NOTP XT(K) = 0.0 40 CONTINUE DO 50 K=1,NOSP FS(K) = FNI*FS(K) 50 CONTINUE DO 60 K=1,NOP IS(K) = IS(K) - 1 IIS(K) = 1 60 CONTINUE IF (MODE.EQ.3 .OR. SUBSET) GO TO 80 DO 70 K=1,NOTP IT(K) = IT(K) - 1 IIT(K) = 1 70 CONTINUE 80 CONTINUE C C START THE MAIN LOOP C C TAKE THE FOURIER TRANSFORM C ND2P1 = N/2 + 1 DO 350 L=1,ND2P1 LM2 = L - 2 DO 90 J=1,M RX(J) = 0. AIX(J) = 0. 90 CONTINUE IF (L.EQ.2) GO TO 110 IF (L.GT.2) GO TO 130 DO 100 K=1,NOSP JSK = JS(K) RX(JSK) = RX(JSK) + FS(K) 100 CONTINUE D = 2. + CH2 GO TO 160 110 CONTINUE DO 120 K=1,NOSP JSK = JS(K) RX(JSK) = RX(JSK) - FLOAT(2*IS(K)-1-4*(IS(K)/2))*FS(K) 120 CONTINUE D = 6. + CH2 GO TO 160 130 CONTINUE DO 140 K=1,NOP IIS(K) = IIS(K) + IS(K) IF (IIS(K).GT.N) IIS(K) = IIS(K) - N 140 CONTINUE DO 150 K=1,NOSP JSK = JS(K) IISK = IIS(K) RX(JSK) = RX(JSK) + COSAJN(IISK)*FS(K) AIX(JSK) = AIX(JSK) - SAJN(IISK)*FS(K) 150 CONTINUE D = 4. + CH2 - 2.*COS(FLOAT(LM2)*PIN) 160 CONTINUE C C SOLVE THE TRIDIAGONAL SYSTEM WITH C OFF-DIAGONAL ELEMENTS EQUAL TO -1.0 , AND C DIAGONAL ELEMENTS D(I) EQUAL TO D, EXCEPT C D(1)=D(M)=D/2+SQRT((D/2)**2-1.0) C M1 = M - 1 IF (ABS(D).LE.2.0) GO TO 190 C C SIMPLE FACTORIZATION AS ABS(D(I)).GT.2 C ROOT = SQRT(.25*D*D-1.) BEI = .5*D - ROOT C FORWARD SUBSTITUTION DO 170 I=2,M RX(I) = RX(I) + RX(I-1)*BEI AIX(I) = AIX(I) + AIX(I-1)*BEI 170 CONTINUE RX(M) = RX(M)/(2.*ROOT) AIX(M) = AIX(M)/(2.*ROOT) C BACKWARD SUBSTITIUTION DO 180 I=1,M1 MI = M - I RX(MI) = (RX(MI)+RX(MI+1))*BEI AIX(MI) = (AIX(MI)+AIX(MI+1))*BEI 180 CONTINUE GO TO 230 190 CONTINUE C C -2.0 < D(I) < 2.0 C USE A STABLE RECURSION FORMULA, SEE W.PROSKUROWSKI AND C O.WIDLUND - MATH.COMP.,VOL.30 (1976), P.449. C S = D S1 = 1.0 V1 = RX(2) W1 = AIX(2) DO 200 I=3,M1 V1 = V1 + S*RX(I) W1 = W1 + S*AIX(I) S2 = S1 S1 = S S = D*S1 - S2 200 CONTINUE S = D S1 = 1.0 V2 = RX(1) + RX(3) W2 = AIX(1) + AIX(3) DO 210 I=4,M V2 = V2 + S*RX(I) W2 = W2 + S*AIX(I) S2 = S1 S1 = S S = D*S1 - S2 210 CONTINUE TR = RX(2) TI = AIX(2) RX(1) = -(V1+S*RX(M))/2.0 AIX(1) = -(W1+S*AIX(M))/2.0 RX(2) = -V2/2.0 AIX(2) = -W2/2.0 DO 220 I=3,M TR1 = RX(I) TI1 = AIX(I) RX(I) = D*RX(I-1) - RX(I-2) - TR AIX(I) = D*AIX(I-1) - AIX(I-2) - TI TR = TR1 TI = TI1 220 CONTINUE C C END OF THE TRIDAGONAL SOLVER C 230 CONTINUE IF (MODE.NE.3) GO TO 270 C SAVE IN F THE SUM OF THE PRESENT AND PREVIOUS C INTERMEDIATE RESULT. IF (L.GT.2) GO TO 250 DO 240 K=1,M F(L,K) = RX(K) + F(L,K) 240 CONTINUE GO TO 350 250 CONTINUE I = 2*L - 3 DO 260 K=1,M F(I,K) = RX(K) + F(I,K) F(I+1,K) = AIX(K) + F(I+1,K) 260 CONTINUE GO TO 350 C C TAKE THE INVERSE FOURIER TRANSFORM C 270 CONTINUE IF (L.EQ.2) GO TO 290 IF (L.GT.2) GO TO 310 DO 280 K=1,NOTP JTK = JT(K) XT(K) = XT(K) + RX(JTK)*0.5 280 CONTINUE GO TO 350 290 CONTINUE DO 300 K=1,NOTP JTK = JT(K) XT(K) = XT(K) - FLOAT(2*IT(K)-1-4*(IT(K)/2))*RX(JTK)*0.5 300 CONTINUE GO TO 350 310 CONTINUE IF (SUBSET) GO TO 330 DO 320 K=1,NOTP IIT(K) = IIT(K) + IT(K) IF (IIT(K).GT.N) IIT(K) = IIT(K) - N 320 CONTINUE 330 CONTINUE DO 340 K=1,NOTP JTK = JT(K) IITK = IIT(K) XT(K) = XT(K) + COSAJN(IITK)*RX(JTK) - SAJN(IITK)*AIX(JTK) 340 CONTINUE 350 CONTINUE C C THE MAIN LOOP ENDS C DO 360 K=1,NOTP XT(K) = XT(K) + XT(K) 360 CONTINUE C RESTORE IS,IT AND FS TO THEIR ORIGINAL VALUES. DO 370 K=1,NOSP FS(K) = FS(K)*FLOAT(N) 370 CONTINUE DO 380 K=1,NOP IS(K) = IS(K) + 1 380 CONTINUE IF (MODE.NE.3) GO TO 390 C C MODE = 3 C TAKE INVERSE FAST FOURIER TRANSFORM C CALL RFORT(F, NDIM, LOG2N, SI, +2, M, IB) C C EXIT FOR MODE=3 C RETURN C 390 CONTINUE IF (SUBSET) GO TO 410 DO 400 K=1,NOTP IT(K) = IT(K) + 1 400 CONTINUE 410 CONTINUE C C EXIT FOR MODE=1 C RETURN C 420 CONTINUE C C ************************************ C THIS PART IS A SOLVER WHICH USES FFT. C ************************************ C C MODE = 2 OR 4 C C TAKE FAST FOURIER TRANSFORMATION OF F C CALL RFORT(F, NDIM, LOG2N, SI, -2, M, IB) C C SOLVE THE TRIDIAGONAL SYSTEM WITH C OFF-DIAGONAL ELEMENTS EQUAL TO -1.0 , AND C DIAGONAL ELEMENTS D(I) EQUAL TO D, EXCEPT C D(1)=D(M)=D/2+SQRT((D/2)**2-1.0) C N2 = N/2 MM1 = M - 1 DO 520 L=1,N2 CS = COS(PIN*FLOAT(L-1)) DO 510 IJ=1,2 I = 2*L + IJ - 2 DO 430 J=1,M RX(J) = F(I,J) 430 CONTINUE D = 2.0 + CH2 IF (I.EQ.2) D = D + 4.0 IF (I.GT.2) D = D + 2.0 - 2.0*CS IF (ABS(D).LE.2.0) GO TO 460 C C SIMPLE FACTORIZATION AS ABS(D(I)).GT.2 C S = SQRT(.25*D*D-1.0) IF (D.GT.0.0) S = -S B = 0.5*D + S C FORWARD SUBSTITUTION DO 440 J=2,M RX(J) = RX(J) + RX(J-1)*B 440 CONTINUE RX(M) = RX(M)/(-2.0*S) C BACKWARD SUBSTITIUTION DO 450 JJ=1,MM1 J = M - JJ RX(J) = (RX(J)+RX(J+1))*B 450 CONTINUE GO TO 490 460 CONTINUE C C -2.0 < D(I) < 2.0 C USE A STABLE RECURSION FORMULA, SEE W.PROSKUROWSKI AND C O.WIDLUND - MATH.COMP.,VOL.30 (1976), P.449. C S1 = 1.0 S = D V1 = RX(2) V2 = RX(1) + RX(3) DO 470 J=3,MM1 V1 = V1 + S*RX(J) V2 = V2 + S*RX(J+1) S2 = S1 S1 = S S = D*S1 - S2 470 CONTINUE T = RX(2) RX(1) = -(V1+S*RX(M))/2.0 RX(2) = -V2/2.0 DO 480 J=3,M T1 = RX(J) RX(J) = D*RX(J-1) - RX(J-2) - T T = T1 480 CONTINUE 490 CONTINUE C DO 500 J=1,M F(I,J) = RX(J) 500 CONTINUE 510 CONTINUE 520 CONTINUE C C END OF THE TRIDAGONAL SOLVER C C C *** TAKE INVERSE FOURIER TRANSFORM. C IF (MODE.NE.4) GO TO 530 C C MODE=4 C CALL RFORT(F, NDIM, LOG2N, SI, +2, M, IB) C C EXIT FOR MODE=4 C RETURN C 530 CONTINUE C C MODE=2 C DO 540 K=1,NOTP IT(K) = IT(K) - 1 IIT(K) = 1 XT(K) = 0. 540 CONTINUE ND2P1 = N/2 + 1 C C START THE MAIN LOOP C DO 650 L=1,ND2P1 IF (MODE.NE.2) GO TO 580 IF (L.GT.2) GO TO 560 DO 550 K=1,M RX(K) = F(L,K) 550 CONTINUE GO TO 580 560 CONTINUE I = 2*L - 3 DO 570 K=1,M RX(K) = F(I,K) AIX(K) = F(I+1,K) 570 CONTINUE 580 CONTINUE IF (L.GT.2) GO TO 620 IF (L.EQ.2) GO TO 600 DO 590 K=1,NOTP JTK = JT(K) XT(K) = XT(K) + RX(JTK)*0.5 590 CONTINUE GO TO 650 600 CONTINUE DO 610 K=1,NOTP JTK = JT(K) XT(K) = XT(K) - FLOAT(2*IT(K)-1-4*(IT(K)/2))*RX(JTK)*0.5 610 CONTINUE GO TO 650 620 CONTINUE DO 630 K=1,NOTP IIT(K) = IIT(K) + IT(K) IF (IIT(K).GT.N) IIT(K) = IIT(K) - N 630 CONTINUE DO 640 K=1,NOTP JTK = JT(K) IITK = IIT(K) XT(K) = XT(K) + COSAJN(IITK)*RX(JTK) - SAJN(IITK)*AIX(JTK) 640 CONTINUE 650 CONTINUE DO 660 K=1,NOTP XT(K) = XT(K) + XT(K) 660 CONTINUE DO 670 K=1,NOTP IT(K) = IT(K) + 1 670 CONTINUE C C EXIT FOR MODE=2 C RETURN END SUBROUTINE RFORT(A, NDIM, M, S, IFS, MM, IB) RFO 10 DIMENSION A(1), S(1), IB(1) C C THIS SUBROUTINE COMPUTES FFT OF LENGTH N FOR REAL DATA C AND IS A MODIFICATION OF THE REAL FOURIER TRANSFORM PROGRAM C BY DR.J.COOLEY. C VECTOR A IS OF LENGTH NDIM*MM, WHERE NDIM>=N, N=2**M, M=4,5,..., C NDIM IS THE FIRST DIMENSION OF THE TWO-DIMENSIONAL ARRAY A C AS IT APPEARS IN THE PROGRAM CALLING RFORT, AND MM=1,2,... C THE ORIGINAL PROGRAM ALLOWED FOR MM=1, ONLY. C ADDITIONAL CHANGES ALLOW TO STORE ALL FOURIER COEFFICIENTS C IN A VECTOR (1=N) = COS $ SIN AMPLITUDES OF GROWING FREQUENCY C I=1 TO N/2-1 ARE IN POS.2*I+1,2*I+2;COS(0) $ COS(N/2) IN POS 1 $ 2. C IN ADDITION SOME INDEXES USED IN BINARY REORDERING ARE C SAVED IN THE VECTOR IB OF LENGTH N. C VECTOR S OF LENGTH N/4 CONTAINS A TABLE OF SINES. C IPS IS A PARAMETER TO BE USED AS FOLLOWS: C = 0 TO SET UP A SINE TABLE, C = 1 TO SET UP A SINE TABLE, AND DO FOURIER SYNTHESIS, C =-1 TO SET UP A SINE TABLE, AND DO FOURIER ANALYSIS, C = 2 TO DO FOURIER SYNTHESIS, C =-2 TO DO FOURIER ANALYSIS. C C THIS PROGRAM USES THE SUBROUTINE FORT TO COMPUTE COMPLEX C FOURIER TRANSFORM OF REAL DATA. C IF (IFS.NE.0) GO TO 10 C C GENERATE A TABLE OF SINES, S, AND REPRESENTATION OF A C PERMUTATION USED IN BINARY REORDERING OF THE DATA. C CALL FORT(A, NDIM, M, S, 0, MM, IB) C RETURN 10 CONTINUE N = 2**M N2 = 2*N NV2 = N/2 NV2M2 = NV2 - 2 MM1 = M - 1 NP = N MP = M KD = NP/N NPV4 = NP/4 NDIFF = NDIM - N IF (IFS.GT.0) GO TO 40 C C COMPUTE FORWARD FFT SIMULTANEOUSLY ON EACH SUBARRAYS OF LENGTH N. C CALL FORT(A, NDIM, MM1, S, -2, MM, IB) C KMIN = 2 KMAX = NV2M2 LN = N DO 30 L=1,MM KT = KD DO 20 K=KMIN,KMAX,2 J = LN - K A1R = A(K+1) + A(J+1) A1I = A(K+2) - A(J+2) A2R = A(K+2) + A(J+2) A2I = A(J+1) - A(K+1) KKT = NPV4 - KT AWR = A2R*S(KKT) + A2I*S(KT) AWI = A2I*S(KKT) - A2R*S(KT) A(K+1) = (A1R+AWR)*0.25 A(K+2) = (A1I+AWI)*0.25 A(J+1) = (A1R-AWR)*0.25 A(J+2) = (AWI-A1I)*0.25 KT = KT + KD 20 CONTINUE T = A(KMIN-1) A(KMIN-1) = (T+A(KMIN))*0.5 A(KMIN) = (T-A(KMIN))*0.5 NK = NV2 + KMIN NK1 = NV2 + KMIN - 1 A(NK1) = .5*A(NK1) A(NK) = -.5*A(NK) KMIN = KMIN + N + NDIFF KMAX = KMAX + N + NDIFF LN = LN + N2 + NDIFF*2 30 CONTINUE RETURN 40 CONTINUE C C COMPUTE INVERSE FFT SIMULTANEOUSLY ON EACH SUBARRAYS OF LENGTH N. C KMIN = 2 KMAX = NV2M2 LN = N DO 60 L=1,MM KT = KD DO 50 K=KMIN,KMAX,2 J = LN - K A1R = A(K+1) + A(J+1) A1I = A(K+2) - A(J+2) AWR = A(K+1) - A(J+1) AWI = A(K+2) + A(J+2) KKT = NPV4 - KT A2R = AWR*S(KKT) - AWI*S(KT) A2I = AWR*S(KT) + AWI*S(KKT) A(K+1) = A1R - A2I A(K+2) = A1I + A2R A(J+1) = A1R + A2I A(J+2) = A2R - A1I KT = KT + KD 50 CONTINUE T = A(KMIN-1) Z = A(KMIN) A(KMIN-1) = T + Z A(KMIN) = T - Z NK = NV2 + KMIN NK1 = NV2 + KMIN - 1 A(NK1) = 2.0*A(NK1) A(NK) = -2.0*A(NK) KMIN = KMIN + N + NDIFF KMAX = KMAX + N + NDIFF LN = LN + N2 + NDIFF*2 60 CONTINUE C CALL FORT(A, NDIM, MM1, S, 2, MM, IB) C RETURN END SUBROUTINE FORT(A, NDIM, M, S, IFS, MM, IB) FOR 10 DIMENSION A(1), S(1), IB(1) C C THIS SUBROUTINE COMPUTES THE COMPLEX FFT. C SEE THE COMMENTS TO SUBROUTINE RFORT. C N = 2**M IF (IFS.NE.0) GO TO 100 C C GENERATE A TABLE OF SINES, S, AND REPRESENTATION OF A C PERMUTATION USED IN BINARY REORDERING OF THE DATA. C THETA = .78539816339745 NT = N/4 MT = M - 2 IF (MT.LE.0) GO TO 90 JSTEP = NT JDIF = NT/2 S(JDIF) = SIN(THETA) IF (MT.LT.2) GO TO 30 DO 20 L=2,MT THETA = THETA*0.5 JSTEP2 = JSTEP JSTEP = JDIF JDIF = JDIF/2 S(JDIF) = SIN(THETA) JC1 = NT - JDIF S(JC1) = COS(THETA) JLAST = NT - JSTEP2 IF (JLAST.LT.JSTEP) GO TO 20 DO 10 J=JSTEP,JLAST,JSTEP JC = NT - J JD = J + JDIF S(JD) = S(J)*S(JC1) + S(JDIF)*S(JC) 10 CONTINUE 20 CONTINUE 30 CONTINUE DO 40 I=1,N IB(I) = 0 40 CONTINUE N2 = N/2 J = 2 NM2 = N - 2 DO 80 I=2,NM2,2 IF (I.GE.J) GO TO 50 IB(I) = J 50 K = N2 60 IF (K.GE.J) GO TO 70 J = J - K K = K/2 GO TO 60 70 J = J + K 80 CONTINUE 90 CONTINUE RETURN 100 CONTINUE C C TAKE BIT REVERSAL C N2 = 2*N NT = N/2 NDIFF = NDIM - N2 DO 120 I=2,N2,2 IF (IB(I).EQ.0) GO TO 120 IR = 0 DO 110 L=1,MM J = IB(I) + IR K = I + IR T = A(K) A(K) = A(J) A(J) = T T = A(K-1) A(K-1) = A(J-1) A(J-1) = T IR = IR + N2 + NDIFF 110 CONTINUE 120 CONTINUE IF (IFS.GT.0) GO TO 150 C C MAKE FOURIER ANALYSIS C FN = N FN = 1.0/FN IMIN = 2 IMAX = N2 DO 140 J=1,MM DO 130 I=IMIN,IMAX,2 A(I-1) = A(I-1)*FN A(I) = -A(I)*FN 130 CONTINUE IMIN = IMIN + N2 + NDIFF IMAX = IMAX + N2 + NDIFF 140 CONTINUE 150 CONTINUE IMIN = 2 IMAX = N2 DO 170 J=1,MM DO 160 I=IMIN,IMAX,4 T = A(I-1) A(I-1) = T + A(I+1) A(I+1) = T - A(I+1) T = A(I) A(I) = T + A(I+2) A(I+2) = T - A(I+2) 160 CONTINUE IMIN = IMIN + N2 + NDIFF IMAX = IMAX + N2 + NDIFF 170 CONTINUE LEXP1 = 2 LEXP = 8 NPL = 2**(M-1) DO 240 L=2,M IMIN = 2 IMAX = N2 DO 190 J=1,MM DO 180 I=IMIN,IMAX,LEXP I1 = I + LEXP1 I2 = I1 + LEXP1 I3 = I2 + LEXP1 T = A(I-1) A(I-1) = T + A(I2-1) A(I2-1) = T - A(I2-1) T = A(I) A(I) = T + A(I2) A(I2) = T - A(I2) T = -A(I3) TI = A(I3-1) A(I3-1) = A(I1-1) - T A(I3) = A(I1) - TI A(I1-1) = A(I1-1) + T A(I1) = A(I1) + TI 180 CONTINUE IMIN = IMIN + N2 + NDIFF IMAX = IMAX + N2 + NDIFF 190 CONTINUE IF (L.EQ.2) GO TO 230 LEN = N2 + NDIFF MLEN = MM*LEN JMAX = LEXP1 DO 220 JMIN=4,MLEN,LEN KLAST = N2 - LEXP JJ = NPL DO 210 J=JMIN,JMAX,2 NPJJ = NT - JJ UR = S(NPJJ) UI = S(JJ) ILAST = J + KLAST DO 200 I=J,ILAST,LEXP I1 = I + LEXP1 I2 = I1 + LEXP1 I3 = I2 + LEXP1 T = A(I2-1)*UR - A(I2)*UI TI = A(I2-1)*UI + A(I2)*UR A(I2-1) = A(I-1) - T A(I2) = A(I) - TI A(I-1) = A(I-1) + T A(I) = A(I) + TI T = -A(I3-1)*UI - A(I3)*UR TI = A(I3-1)*UR - A(I3)*UI A(I3-1) = A(I1-1) - T A(I3) = A(I1) - TI A(I1-1) = A(I1-1) + T A(I1) = A(I1) + TI 200 CONTINUE JJ = JJ + NPL 210 CONTINUE JMAX = JMAX + N2 + NDIFF 220 CONTINUE 230 LEXP1 = 2*LEXP1 LEXP = 2*LEXP NPL = NPL/2 240 CONTINUE IF (IFS.GT.0) RETURN IMIN = 2 IMAX = N2 DO 260 J=1,MM DO 250 I=IMIN,IMAX,2 A(I) = -A(I) 250 CONTINUE IMIN = IMIN + N2 + NDIFF IMAX = IMAX + N2 + NDIFF 260 CONTINUE RETURN END SUBROUTINE DECOMP(C, NDIM, IP, SCALES, IPS, IFLAG) DEC 10 DIMENSION C(NDIM,IP) DIMENSION SCALES(IP), IPS(IP) C C *** LU DECOMPOSISION OF MATRIX C OF ORDER IP BY IP C *** AFTER FORSYTHE-MOLER, PRENTICE (1967). C *** THE RESULT LU OVERWRITES THE INPUT MATRIX C. C *** NDIM IS THE FIRST DIMENSION OF MATRIX C AS C *** IT APPEARS IN THE PROGRAM CALLING DECOMP. C *** SCALES AND IPS ARE HELP VECTORS WHICH VALUES C *** MUST NOT BE DESTROYED BEFORE CALLING SOLVE. C *** IFLAG IS AN ERROR INDICATOR. IF MATRIX C C *** IS COMPUTATIONALLY SINGULAR THEN IFLAG=1. C IFLAG = 0 DO 50 I=1,IP IPS(I) = I ROWNRM = 0E0 DO 20 J=1,IP IF (ROWNRM-ABS(C(I,J))) 10, 20, 20 10 ROWNRM = ABS(C(I,J)) 20 CONTINUE IF (ROWNRM) 30, 40, 30 30 SCALES(I) = 1E0/ROWNRM GO TO 50 40 IFLAG = 1 SCALES(I) = 0E0 50 CONTINUE IPM1 = IP - 1 DO 140 K=1,IPM1 BIG = 0E0 DO 70 I=K,IP IS = IPS(I) SIZE = ABS(C(IS,K))*SCALES(IS) IF (SIZE-BIG) 70, 70, 60 60 BIG = SIZE IDXPIV = I 70 CONTINUE IF (BIG) 90, 80, 90 80 IFLAG = 1 GO TO 140 90 IF (IDXPIV-K) 100, 110, 100 100 J = IPS(K) IPS(K) = IPS(IDXPIV) IPS(IDXPIV) = J 110 KP = IPS(K) PIVOT = C(KP,K) KP1 = K + 1 DO 130 I=KP1,IP IS = IPS(I) EM = -C(IS,K)/PIVOT C(IS,K) = -EM DO 120 J=KP1,IP C(IS,J) = C(IS,J) + EM*C(KP,J) 120 CONTINUE 130 CONTINUE 140 CONTINUE KP = IPS(IP) IF (C(KP,IP)) 160, 150, 160 150 IFLAG = 1 160 RETURN END SUBROUTINE SOLVE(C, NDIM, IP, X, W, SCALES, IPS) SOL 10 DIMENSION C(NDIM,IP) DIMENSION X(IP), W(IP), SCALES(IP), IPS(IP) C C *** SOLVING SYSTEM OF EQUATIONS WITH THE LU DECOMPOSED C *** MATRIX C OF ORDER IP BY IP. C *** THE SOLUTION IS IN X , THE RHS IN W. C *** NDIM IS THE FIRST DIMENSION OF MATRIX C AS C *** IT APPEARS IN THE PROGRAM CALLING SOLVE. C IP1 = IP + 1 IS = IPS(1) X(1) = W(IS) DO 20 I=2,IP IS = IPS(I) IM1 = I - 1 SUM = 0E0 DO 10 J=1,IM1 SUM = SUM + C(IS,J)*X(J) 10 CONTINUE X(I) = W(IS) - SUM 20 CONTINUE IS = IPS(IP) X(IP) = X(IP)/C(IS,IP) DO 40 IBACK=2,IP I = IP1 - IBACK IS = IPS(I) IPL1 = I + 1 SUM = 0E0 DO 30 J=IPL1,IP SUM = SUM + C(IS,J)*X(J) 30 CONTINUE X(I) = (X(I)-SUM)/C(IS,I) 40 CONTINUE RETURN END C DRIVER PROGRAM TO ILLUSTRATE THE USE OF CMMIMP MAN 10 C MAN 20 REAL F(20,20), W(300), HX(40), HY(40), TX(40), TY(40), BX(40), MAN 30 * BY(40) MAN 40 INTEGER IW(800), IC(40), JC(40) MAN 50 LOGICAL DIR, LAPL, TEST MAN 60 C MAN 70 C NDIM IS THE FIRST DIMENSION OF THE ARRAY F. MAN 80 C THIS PARAMETER IS USED TO SPECIFY THE VARIABLE MAN 90 C DIMENSION OF F. NDIM MUST NOT BE LESS THAN N. MAN 100 C MAN 110 C LEW IS THE TOTAL DIMENSION OF THE WORKSPACE W. MAN 120 C FOR THE DIRICHLET PROBLEM MAN 130 C LEW MUST NOT BE LESS THAN 6*IP; MAN 140 C FOR THE NEUMANN PROBLEM MAN 150 C LEW MUST NOT BE LESS THAN 7*IP. MAN 160 C MAN 170 C LEIW IS THE TOTAL DIMENSION OF THE WORKSPACE IW. MAN 180 C FOR THE DIRICHLET PROBLEM MAN 190 C LEIW MUST NOT BE LESS THAN 19*IP; MAN 200 C FOR THE NEUMANN PROBLEM MAN 210 C LEIW MUST NOT BE LESS THAN 13*IP. MAN 220 C MAN 230 C IPMAX IS THE DIMENSION OF VECTORS IC,JC,HX,HY,TX,TY,BX,BY MAN 240 C AND IS AN UPPER ESTIMATE OF IP, THE NUMBER OF IRREGULAR MAN 250 C MESH POINTS. MAN 260 C MAN 270 C MAN 280 C ******************** MAN 290 C * INITIALIZE * MAN 300 C ******************** MAN 310 C MAN 320 NDIM = 20 MAN 330 LEW = 300 MAN 340 LEIW = 800 MAN 350 IPMAX = 40 MAN 360 CON = 0.0 MAN 370 DIR = .TRUE. MAN 380 LAPL = .FALSE. MAN 390 TEST = .TRUE. MAN 400 ACCY = 1E-3 MAN 410 ITMAX = 20 MAN 420 N = 16 MAN 430 M = N MAN 440 IOUT = 6 MAN 450 C MAN 460 C THE TEST IS PERFORMED IN THE CIRCULAR REGION MAN 470 C MAN 480 C ( X - 0.5 ) ** 2 + ( Y - 0.5 ) ** 2 = R ** 2, R=0.375, MAN 490 C MAN 500 C WITH THE EXACT SOLUTION BEING U(X,Y) = X**5 + Y**5. MAN 510 C MAN 520 C ******************** MAN 530 C * INPUT DATA * MAN 540 C ******************** MAN 550 C MAN 560 C PROCESS THE INFORMATION ABOUT THE REGION AND ITS BOUNDARY. MAN 570 C MAN 580 C INPUT DATA ABOUT THE REGION: DEFINE IP,IPP,H, MAN 590 C IC,JC,HX,HY AND TX,TY (FOR THE NEUMANN PROBLEM). MAN 600 C MAN 610 CALL DOMAIN(N, M, IC, JC, HX, HY, TX, TY, IPMAX, IP, IPP, H, DIR) MAN 620 C MAN 630 C INPUT BOUNDARY DATA: DEFINE BX AND BY. MAN 640 C MAN 650 CALL BNDRY(IC, JC, HX, HY, BX, BY, IP, IPP, H, DIR) MAN 660 C MAN 670 C INPUT RIGHT SIDE: DEFINE F MAN 680 C MAN 690 CALL CHARGE(F, NDIM, N, M, H, CON) MAN 700 C MAN 710 C ******************** MAN 720 C * CALL CMMIMP * MAN 730 C ******************** MAN 740 C MAN 750 CALL CMMIMP(IC, JC, HX, HY, TX, TY, BX, BY, F, NDIM, N, M, H, MAN 760 * CON, DIR, LAPL, TEST, ACCY, ITMAX, IPMAX, IP, IPP, W, LEW, IW, MAN 770 * LEIW, IOUT, IERROR) MAN 780 C MAN 790 IF (IERROR.GT.0) WRITE (IOUT,99999) IERROR MAN 800 C MAN 810 C ******************** MAN 820 C * PRINTOUT * MAN 830 C ******************** MAN 840 C MAN 850 C FOR THE TEST EXAMPLE MAN 860 C PRINT OUT THE NUMBER OF CORRECT DIGITS MAN 870 C FOR ALL POINTS INSIDE THE IRREGULAR REGION. MAN 880 C MAN 890 IF (TEST .AND. IERROR.EQ.0) CALL OUTPUT(F, NDIM, N, M, H, DIR, MAN 900 * IOUT) MAN 910 C MAN 920 C NOTE: SUBROUTINES DOMAIN, BNDRY, CHARGE AND OUTPUT ARE USED MAN 930 C ONLY IN THIS DRIVER AND ARE NOT REQUIRED IN CMMIMP. MAN 940 C MAN 950 99999 FORMAT (//1X, 19H DETECTED ERROR NO., I3//) MAN 960 STOP MAN 970 END MAN 980