SUBROUTINE IESIMP(KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, IES 10 * NUPPER, MUPPER, W, IFLGR2, IER) C THE INTEGRAL EQUATION BEING SOLVED IS C B C X(S) - INT KERNEL(S,T)*X(T)*DT = RHFCN(S) C A C THE METHOD BEING USED IS BASED ON THE NYSTROM METHOD WITH C SIMPSONS RULE, WITH AN ITERATIVE TECHNIQUE OF SOLUTION FOR THE C RESULTING LINEAR SYSTEM. C KERNEL THESE ARE DOUBLE PRECISION FUNCTIONS OF TWO AND ONE C RHFCN VARIABLES, RESPECTIVELY. THEY MUST BE DECLARED IN AN C EXTERNAL STATEMENT IN THE PROGRAM CALLING IESIMP. C NZ THE INITIAL VALUE OF N IN THE PROGRAM. N IS THE ORDER C OF AN INVERSE MATRIX WHICH IS BEING USED TO C ITERATIVELY SOLVE A LARGER SYSTEM OF ORDER M, WHICH C APPROXIMATES THE ABOVE INTEGRAL EQUATION. THE CHOICE C OF NZ MUST ALWAYS BE ODD AND GREATER THAN TWO. FOR A C DEFAULT CHOICE, SET NZ=3. C ON COMPLETION OF THE PROGRAM, NZ IS SET EQUAL TO THE C FINAL NUMBER OF NODE POINTS USED IN OBTAINING THE C SOLUTION. THE ARRAYS W AND X WILL CONTAIN THE NZ C FINAL NODE POINTS AND CORRESPONDING SOLUTION C VALUES. C EPS THE DESIRED ERROR, INTERPRETED ACCORDING TO IFLAG. C THIS VARIABLE IS CHANGED ON COMPLETION OF THE C PROGRAM. SEE THE DISCUSSION OF IFLAG AND IER FOR MORE C INFORMATION. C IFLAG =0 EPS IS TO BE INTERPRETED AS AN ABSOLUTE ERROR. C =1 EPS IS TO BE INTERPRETED AS A RELATIVE ERROR. C X THE ANSWER TO THE INTEGRAL EQUATION, EVALUATED AT THE C FINAL SET OF NODE POINTS, WHICH ARE STORED IN W. THE C DIMENSION OF X SHOULD BE AT LEAST MUPPER. C NUPPER AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM. N C IS THE ORDER OF APPROXIMATE INVERSE WHICH IS BEING C USED TO ITERATIVELY SOLVE A LARGER SYSTEM OF ORDER M. C MUPPER AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM. C M IS THE NUMBER OF NODE POINTS BEING USED IN THE C SIMPSON RULE APPROXIMATION OF THE INTEGRAL OPERATOR. C N AND M ARE ALWAYS ODD INTEGERS, GREATER THAN TWO. C W THIS IS TEMPORARY WORKING STORAGE. IT SHOULD HAVE C DIMENSION 4*NU*NU+10*NU+9*MU, WITH NU=NUPPER, C MU=MUPPER. ON EXIT, W WILL CONTAIN THE FINAL CHOICE C OF NODE POINTS,CORRESPONDING TO THE APPROXIMATE C SOLUTION VALUES STORED IN X. C IFLGR2 =0 THE NORMAL SETTING. C =1 IF THE INTEGRAND K(S,T)*X(T), REGARDED AS A FUNCTION C OF T, IS KNOWN TO BE PERIODIC ALONG WITH A NUMBER OF C ITS DERIVATIVES ON THE INTERVAL (A,B), FOR ALL S IN C THE INTERVAL, THEN THE ERROR ESTIMATES WILL BE MORE C ACCURATE IF IFLGR2=1. THIS WILL GENERALLY RESULT IN C MORE RAPID CONVERGENCE OF THE PROGRAM. C IER =0 THIS ERROR COMPLETION CODE MEANS THAT THE ERROR TEST C WAS SATISFIED. THE PREDICTED ERROR IS STORED IN EPS. C =1 THE ERROR TEST WAS NOT SATISFIED. THE PREDICTED ERROR C IS IN EPS. C =2 THE ERROR TEST WAS NOT SATISFIED. THE VARIABLE EPS C HAS BEEN SET TO ZERO. C REGARDLESS OF THE VALUE OF IER, THE MOST ACCURATE C SOLUTION OBTAINED IS STORED IN X. C IESIMP IS PRESENTLY LIMITED TO NUPPER .LE. 100. THIS IS ENTIRELY C DUE TO A RESTRICTION IN THE PROGRAM LINSYS. TO REMOVE THIS C RESTRICTION, CHANGE THE DIMENSION STATEMENT IN THE C COMMON/LINEAR/ STATEMENT IN LINSYS. DOUBLE PRECISION KERNEL, RHFCN, A, B, EPS, X, W, CUTOFF, * RATIO, ROOTRT DIMENSION X(1), W(1) EXTERNAL KERNEL, RHFCN C ******************************************************************* C * COMMON /INFO/ R1, R2, NFINAL C * C THE VARIABLES IN COMMON/INFO/ GIVE ADDITIONAL INFORMATION ABOUT* C THE FUNCTIONING OF IESIMP. R1 GIVES THE FINAL ITERATIVE RATE * C OF CONVERGENCE FOR SOLVING THE LINEAR SYSTEMS OF ORDER M. R2 * C GIVES THE RATE OF CONVERGENCE OF THE NYSTROM METHOD. FOR A * C SMOOTH KERNEL FUNCTION AND SMOOTH SOLUTION FUNCTION X(S), THE * C VALUE OF R2 SHOULD BE 0.0625 FOR SIMPSONS RULE. NFINAL GIVES * C THE FINAL VALUE OF N USED IN ITERATIVELY SOLVING THE LARGER * C SYSTEMS OF ORDER M. IF THE VALUE OF NFINAL EQUALS THE FINAL * C VALUE OF NZ, THEN ITERATION WAS NOT INVOKED SUCCESSFULLY. * C * C ******************************************************************* DATA CUTOFF /0.5D0/ C ******************************************************************* C * DATA RATIO /0.0625D0/, ROOTRT /0.25D0/ C * C THESE CONSTANTS DEPEND ON THE NUMERICAL INTEGRATION RULE * C BEING USED. SEE THE DISCUSSION OF IESIMP IN ORDER TO SET THEM * C PROPERLY WHEN CHANGING NUMERICAL INTEGRATION RULES. RATIO * C DEPENDS ON THE RATE OF CONVERGENCE OF THE RULE BEING USED. * C ROOTRT USUALLY EQUALS SQRT(RATIO). * C * C ******************************************************************* C SET UP RELATIVE BASE ADDRESSES FOR THE VARIOUS ARRAYS INTO WHICH C W IS TO BE SPLIT. N = NUPPER M = MUPPER NSQ = N*N I1 = 1 I2 = I1 + NSQ I3 = I2 + NSQ I4 = I3 + (NSQ+N)/2 I5 = I4 + (NSQ+N)/2 I6 = I5 + M I7 = I6 + M I8 = I7 + N I9 = I8 + N I10 = I9 + M I11 = I10 + N I12 = I11 + M I13 = I12 + M I14 = I13 + M I15 = I14 + N I16 = I15 + M I17 = I16 + M I18 = I17 + N I19 = I18 + M I20 = I19 + 4*N NHALF = (N+1)/2 CALL IESP(KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, NUPPER, * MUPPER, IFLGR2, IER, RATIO, ROOTRT, CUTOFF, NHALF, W, W(I1), * W(I2), W(I3), W(I4), W(I5), W(I6), W(I7), W(I8), W(I9), * W(I10), W(I11), W(I12), W(I13), W(I14), W(I15), W(I16), * W(I17), W(I18), W(I19), W(I20)) RETURN END SUBROUTINE IESP(KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, NUP, IES 10 * MUP, IFLGR2, IER, RATIO, ROOTRT, CUTOFF, NHALF, T, LUFACT, * KMM, KMN, KNM, RHS, R, RH, DELN, TM, TN, XM, XMZ, WM, WN, * OLDX, SAVE, XN, SAVE2, ASIDE, ASIDE3) C THIS ROUTINE CONTROLS THE SOLUTION OF THE INTEGRAL EQUATION. DOUBLE PRECISION KERNEL, RHFCN, A, B, EPS, X, RATIO, ROOTRT, * CUTOFF, T, LUFACT, KMM, KMN, KNM, RHS, R, RH, DELN, TM, TN, * XM, XMZ, WM, WN, OLDX, SAVE, SAVE2, XN, ASIDE, ASIDE2, * ASIDE3, R2, R1RAT, NUMR2, DENR2, NORM, DMIN1, DSQRT, ERROR, * TEMP, TEST, R1, DENR1, NUMR1, RT, DET INTEGER OLDM, TWICE, FLAG DIMENSION X(MUP), T(MUP), LUFACT(NUP,NUP), KMM(NUP,NUP), * KMN(NUP,NHALF), KNM(NHALF,NUP), RHS(MUP), R(MUP), RH(NUP), * DELN(NUP), TM(MUP), TN(NUP), XM(MUP), XMZ(MUP), WM(MUP), * WN(NUP), OLDX(MUP), SAVE(MUP), XN(NUP), SAVE2(MUP), * ASIDE(NUP,4), ASIDE2(5), ASIDE3(NUP,NUP) COMMON /INFO/ R1, R2, N EXTERNAL KERNEL, RHFCN C ******************************************************************* C * TWICE(KK) = 2*KK - 1 C * C THIS FUNCTION GIVES THE FORMULA FOR INCREASING THE NUMBER OF * C NODE POINTS. WITH ANOTHER INTEGRATION RULE, IT MAY BE * C NECESSARY TO CHANGE THE DEFINITION OF TWICE(KK). * C * C ******************************************************************* C INITIALIZATION. IER = 0 LOOP = 1 N = NZ R2 = 0.5D0 M = TWICE(N) R1RAT = ROOTRT C STAGE A. DIRECT SOLUTION OF LINEAR SYSTEM (I-KN)*XN=RHS, WHILE C TRYING TO FIND A GOOD APPROXIMATE INVERSE TO IMPLEMENT THE C ITERATIVE METHOD OF SOLUTION. C CREATE NODES TN(I) AND WEIGHTS WN(I), I=1,...,N. CALL WANDT(WN, TN, N, A, B) C CREATE SYSTEM (I-KN)*XN=RHFCN. DO 20 I=1,N DO 10 J=1,N LUFACT(I,J) = -WN(J)*KERNEL(TN(I),TN(J)) 10 CONTINUE XN(I) = RHFCN(TN(I)) LUFACT(I,I) = LUFACT(I,I) + 1.0D0 20 CONTINUE GO TO 60 C THIS IS ENTRANCE FOR AN INCREASED VALUE OF N, USING PREVIOUSLY C STORED VALUES IN KMM AND RHS. 30 DO 50 J=1,N DO 40 I=1,N LUFACT(I,J) = -KMM(I,J) 40 CONTINUE WN(J) = WM(J) TN(J) = TM(J) XN(J) = RHS(J) LUFACT(J,J) = LUFACT(J,J) + 1.0D0 50 CONTINUE C THIS IS THE ENTRANCE WHEN ITERATION IN STAGE B FAILS AND WE NEED C TO INCREASE N TO OBTAIN A BETTER ITERATIVE RATE. 60 CONTINUE C SOLVE (I-KN)*XN=RHFCN AT ALL TN(I).ALSO OBTAIN THE LU C DECOMPOSITION FOR LATER USE IN THE STAGE B ITERATIVE METHOD. C ******************************************************************* C * CALL LINSYS(LUFACT, LUFACT, N, XN, XN, 1, DET, NUP) C * C THIS LINEAR EQUATION SOLVER MAY BE REPLACED BY ANOTHER PROGRAM,* C BUT CARE MUST BE TAKEN THAT THE NEW PROGRAM DOES THE SAME JOB * C AS LINSYS. THE PROGRAM LINSYS IS ALSO REFERENCED IN THE * C SUBROUTINE ITERT. * C * C ******************************************************************* IF (LOOP.EQ.1) GO TO 110 IF (LOOP.EQ.2) GO TO 80 C COMPUTE THE RATE OF CONVERGENCE, AND COMPARE WITH THE C THEORTICAL RATIO. NUMR2 = NORM(XN,OLDX,N,1) R2 = DMIN1(0.5D0,NUMR2/DENR2) IF ((R2.LT.RATIO) .AND. (IFLGR2.EQ.0)) R2 = RATIO R1RAT = DMIN1(DSQRT(R2),ROOTRT) IF (IFLGR2.EQ.0) R1RAT = ROOTRT C CHECK FOR ERROR IN XN USING TEST INVOLVING R2 AND OLDX,ACCORDING C TO THEORY FOR ASYMPTOTIC ERROR BOUNDS. 70 ERROR = (R2/(1.0D0-R2))*NUMR2 IF (IFLAG.EQ.1) ERROR = ERROR/NORM(XN,XN,N,0) IF ((IFLGR2.EQ.1) .AND. (R2.LT.RATIO)) ERROR = 2.0*ERROR IF (ERROR.LE.EPS) GO TO 90 DENR2 = NUMR2 GO TO 110 C ENTRANCE FOR LOOP=2. 80 NUMR2 = NORM(XN,OLDX,N,1) DENR2 = 0.0D0 GO TO 70 C SET UP T,X,EPS,NZ FOR SUCCESSFUL RETURN. 90 DO 100 I=1,N X(I) = XN(I) T(I) = TN(I) 100 CONTINUE EPS = ERROR NZ = N RETURN C INITIALIZE FOR SOLVING (I-KM)*XM=RHFCN ITERATIVELY. 110 CALL WANDT(WM, TM, M, A, B) FLAG = 0 CALL INTERP(TM, WM, XMZ, M, TN, WN, XN, N, KERNEL, RHFCN, * RHS, KMN, NHALF, NUP) DO 120 I=1,M OLDX(I) = XMZ(I) 120 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) DENR1 = NORM(XM,XMZ,M,1) FLAG = 1 DO 130 I=1,M XMZ(I) = XM(I) 130 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) NUMR1 = NORM(XM,XMZ,M,1) C CHECK ON THE SPEED OF CONVERGENCE OF ITERATIVE METHOD. IF IT IS C SUFFICIENTLY RAPID, THEN FIX N AND GO TO STAGE B. R1 = NUMR1/DENR1 IF (M.GT.NUP) GO TO 190 IF (R1.LE.R1RAT) GO TO 150 C REINITIALIZE FOR SOLVING (I-KN)*XN=RHFCN AGAIN WITH A LARGER N. 140 N = M LOOP = LOOP + 1 M = TWICE(N) GO TO 30 C SAVE INFORMATION IN CASE STAGE B ABORTS AT A LARGER VALUE OF M C AND STAGE A HAS TO BE RETURNED TO. 150 DO 160 I=1,M ASIDE(I,1) = OLDX(I) ASIDE(I,2) = WM(I) ASIDE(I,3) = TM(I) ASIDE(I,4) = RHS(I) 160 CONTINUE ASIDE2(1) = LOOP ASIDE2(2) = M ASIDE2(3) = R2 ASIDE2(4) = DENR2 ASIDE2(5) = R1RAT DO 180 I=1,M DO 170 J=1,M ASIDE3(I,J) = KMM(I,J) 170 CONTINUE 180 CONTINUE GO TO 240 C STAGE B. ITERATIVE METHOD OF SOLUTION OF (I-KM)*XM=RHS. 190 IF (R1.LE.CUTOFF) GO TO 240 C IF ITERATES ARE CONVERGING VERY SLOWLY OR NOT AT ALL, THEN C RETURN WITHOUT FURTHER ATTEMPTS TO LESSEN ERROR. 200 IF (LOOP.NE.1) GO TO 230 EPS = 0.0D0 IER = 2 210 DO 220 I=1,N X(I) = XN(I) T(I) = TN(I) 220 CONTINUE NZ = N RETURN 230 EPS = ERROR IER = 1 GO TO 210 C TEST TO SEE IF THE CURRENT ITERATE XM IS SUFFICIENTLY ACCURATE C COMPARED TO THE TRUE XM. 240 TEMP = NORM(XM,OLDX,M,1) RT = DMIN1(RATIO,R2) TEST = ((1.0D0-R1)/R1)*(RT/(1.0D0-RT))*TEMP IF (NUMR1.LE.TEST) GO TO 320 C ITERATE NOT ACCURATE. INITIALIZE FOR COMPUTATION OF ANOTHER C ITERATE. DENR1 = NUMR1 DO 250 I=1,M XMZ(I) = XM(I) 250 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) NUMR1 = NORM(XM,XMZ,M,1) R1 = NUMR1/DENR1 IF (R1.LE.CUTOFF) GO TO 240 C THIS IS ENTRANCE FOR CASE WHERE ITERATION FAILS IN STAGE B. C PARAMETERS MUST BE RESET FOR A RETURN TO STAGE A, OR IF N CANNOT C BE INCREASED, FOR AN ERROR EXIT FROM IESIMP. 260 MNEW = ASIDE2(2) IF (MNEW.GT.NUP) GO TO 290 IF (M.EQ.TWICE(N)) GO TO 140 N = MNEW DO 280 J=1,N DO 270 I=1,N LUFACT(I,J) = -ASIDE3(I,J) 270 CONTINUE OLDX(J) = ASIDE(J,1) WN(J) = ASIDE(J,2) TN(J) = ASIDE(J,3) XN(J) = ASIDE(J,4) LUFACT(J,J) = LUFACT(J,J) + 1.0D0 280 CONTINUE M = TWICE(N) LOOP = ASIDE2(1) + 1.0D0 R2 = ASIDE2(3) DENR2 = ASIDE2(4) R1RAT = ASIDE2(5) GO TO 60 C ABORTIVE EXIT FROM STAGE B. N CANNOT BE INCREASED FURTHER, AND C R1 IS NOT SUFFICIENTLY SMALL. 290 IF (M.EQ.TWICE(N)) GO TO 200 DO 300 I=1,OLDM T(I) = SAVE(I) 300 CONTINUE NZ = OLDM IF (LOOP.EQ.1) GO TO 310 EPS = ERROR IER = 1 RETURN 310 EPS = 0.0 IER = 2 RETURN C AN ACCURATE VALUE OF XM HAS BEEN OBTAINED. R2 IS TO BE COMPUTED C AND COMPARED TO RATIO. THEN THE ERROR IN XM IS TO BE ESTIMATED. 320 IF (LOOP.EQ.1) GO TO 340 NUMR2 = TEMP R2 = DMIN1(NUMR2/DENR2,0.5D0) IF ((R2.LT.RATIO) .AND. (IFLGR2.EQ.0)) R2 = RATIO DENR2 = NUMR2 330 ERROR = (R2/(1.0D0-R2))*TEMP IF (IFLAG.EQ.1) ERROR = ERROR/NORM(XM,XM,M,0) IF ((IFLGR2.EQ.1) .AND. (R2.LT.RATIO)) ERROR = 2.0*ERROR IF (ERROR.LE.EPS) GO TO 400 MNEW = TWICE(M) IF (MNEW.LE.MUP) GO TO 350 IER = 1 GO TO 400 340 DENR2 = TEMP LOOP = 2 GO TO 330 C ERROR NOT SUFFICIENTLY SMALL. M IS INCREASED AND TWO MORE C ITERATES ARE COMPUTED WITH A NEW M. 350 DO 360 I=1,M X(I) = XM(I) 360 CONTINUE OLDM = M M = MNEW DO 370 I=1,OLDM SAVE2(I) = WM(I) SAVE(I) = TM(I) 370 CONTINUE CALL WANDT(WM, TM, M, A, B) FLAG = 0 CALL INTERP(TM, WM, XMZ, M, SAVE, SAVE2, XM, OLDM, KERNEL, * RHFCN, RHS, KMN, NHALF, NUP) DO 380 I=1,M OLDX(I) = XMZ(I) 380 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) DENR1 = NORM(XM,XMZ,M,1) FLAG = 1 DO 390 I=1,M XMZ(I) = XM(I) 390 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) NUMR1 = NORM(XM,XMZ,M,1) R1 = NUMR1/DENR1 IF (R1.LE.CUTOFF) GO TO 240 GO TO 260 400 DO 410 I=1,M X(I) = XM(I) T(I) = TM(I) 410 CONTINUE NZ = M EPS = ERROR RETURN END SUBROUTINE WANDT(W, T, N, A, B) WAN 10 C THIS ROUTINE CALCULATES THE WEIGHTS AND NODE POINTS FOR C SIMPSONS RULE ON (A,B), WITH N EVENLY SPACED NODES. DOUBLE PRECISION W, T, A, B, FN, FI, TEMP, H DIMENSION W(N), T(N) C CALCULATE NODE POINTS. FN = N H = (B-A)/(FN-1.0D0) DO 10 I=1,N FI = I T(I) = A + (FI-1.0D0)*H 10 CONTINUE C CALCULATE WEIGHTS. W(1) = H/3.0D0 W(N) = H/3.0D0 NM1 = N - 1 TEMP = 4.0D0*H/3.0D0 DO 20 I=2,NM1,2 W(I) = TEMP 20 CONTINUE IF (N.EQ.3) RETURN TEMP = 2.0D0*H/3.0D0 NM2 = N - 2 DO 30 I=3,NM2,2 W(I) = TEMP 30 CONTINUE RETURN END SUBROUTINE INTERP(TM, WM, XM, M, TN, WN, XN, N, KERNEL, INT 10 * RHFCN, RHS, KMN, NHALF, NUP) C USE THE VALUES OF XN(I), I=1,...,N, TO CALCULATE THE NYSTROM C INTERPOLATES XM(I), I=1,...,M. DOUBLE PRECISION KERNEL, RHFCN, TM, WM, XM, TN, WN, XN, RHS, * KMN DIMENSION TM(M), WM(M), XM(M), TN(N), WN(N), XN(N), RHS(M), * KMN(NUP,NHALF) IF (M.GT.NUP) GO TO 60 C SINCE M .LE. NUPPER, SAVE K(TM(I),TN(J))=KMN(I,J) AND C RHS(I)=RHFCN(TM(I)) FOR LATER USE IN ITERT. DO 20 I=1,M DO 10 J=1,N KMN(I,J) = WN(J)*KERNEL(TM(I),TN(J)) 10 CONTINUE 20 CONTINUE DO 30 I=1,M RHS(I) = RHFCN(TM(I)) XM(I) = RHS(I) 30 CONTINUE C CALCULATE NYSTROM INTERPOLATING FORMULA. DO 50 I=1,M DO 40 J=1,N XM(I) = XM(I) + KMN(I,J)*XN(J) 40 CONTINUE 50 CONTINUE RETURN C M .GT. NUPPER, SO SAVE JUST RHS(I) FOR LATER USE IN ITERT. C CALCULATE NYSTROM INTERPOLATING FORMULA. 60 DO 80 I=1,M RHS(I) = RHFCN(TM(I)) XM(I) = RHS(I) DO 70 J=1,N XM(I) = XM(I) + WN(J)*KERNEL(TM(I),TN(J))*XN(J) 70 CONTINUE 80 CONTINUE RETURN END SUBROUTINE ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, ITE 10 * XMZ, KMM, KMN, KNM, RHS, LUFACT, R, RH, DELN, NUP, NHALF, * IFLG) C THIS ROUTINE CALCULATES ONE ITERATE XM GIVEN THE INITIAL GUESS C XMZ. THE ROUTINE IS DIVIDED ACCORDING TO WHETHER OR NOT C M .GT. NUPPER. DOUBLE PRECISION KERNEL, RHFCN, TN, WN, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, LUFACT, R, RH, DELN, SUM, DET DIMENSION TN(N), WN(N), TM(M), WM(M), XM(M), XMZ(M), * KMM(NUP,NUP), KMN(NUP,NHALF), KNM(NHALF,NUP), RHS(M), * LUFACT(NUP,NUP), R(M), RH(N), DELN(N) C M .GT. NUPPER MEANS THAT THE MATRICES KMM,KMN,KNM CAN NO LONGER C BE STORED DUE TO LACK OF SPACE. IF (M.GT.NUP) GO TO 120 IF (IFLG.EQ.1) GO TO 40 C IF IFLG=0, THEN THE MATRICES KMM AND KNM MUST BE COMPUTED C AND STORED. DO 30 J=1,M DO 10 I=1,M KMM(I,J) = WM(J)*KERNEL(TM(I),TM(J)) 10 CONTINUE DO 20 I=1,N KNM(I,J) = WM(J)*KERNEL(TN(I),TM(J)) 20 CONTINUE 30 CONTINUE C COMPUTE RESIDUALS R(I)=RHFCN(TM(I))-XMZ(I)+KM(TM(I))*XMZ(I) 40 DO 60 I=1,M SUM = 0.0D0 DO 50 J=1,M SUM = SUM + KMM(I,J)*XMZ(J) 50 CONTINUE R(I) = RHS(I) - (XMZ(I)-SUM) 60 CONTINUE C COMPUTE RH=KM*R AT ALL TN(I). DO 80 I=1,N RH(I) = 0.0D0 DO 70 J=1,M RH(I) = RH(I) + KNM(I,J)*R(J) 70 CONTINUE 80 CONTINUE C CALCULATE DELN=((I-KN)**(-1))*KM*R AT ALL TN(I). C ******************************************************************* C * CALL LINSYS(LUFACT, LUFACT, N, RH, DELN, 3, DET, NUP) C * C SEE THE ORIGINAL REFERENCE IN IESP. * C * C ******************************************************************* C CALCULATE NEW XM. DO 110 I=1,M SUM = 0.0D0 DO 90 J=1,M SUM = SUM + KMM(I,J)*R(J) 90 CONTINUE DO 100 J=1,N SUM = SUM + KMN(I,J)*DELN(J) 100 CONTINUE XM(I) = SUM + R(I) + XMZ(I) 110 CONTINUE RETURN C ENTRANCE WHEN M .GT. NUP. C CALCULATE RESIDUALS. 120 DO 140 I=1,M SUM = 0.0D0 DO 130 J=1,M SUM = SUM + WM(J)*KERNEL(TM(I),TM(J))*XMZ(J) 130 CONTINUE R(I) = RHS(I) - (XMZ(I)-SUM) 140 CONTINUE C CALCULATE RH=KM*R. DO 160 I=1,N RH(I) = 0.0D0 DO 150 J=1,M RH(I) = RH(I) + WM(J)*KERNEL(TN(I),TM(J))*R(J) 150 CONTINUE 160 CONTINUE C ******************************************************************* C * CALL LINSYS(LUFACT, LUFACT, N, RH, DELN, 3, DET, NUP) C * C SEE THE ORIGINAL REFERENCE IN IESP. * C * C ******************************************************************* C CALCULATE XM. DO 190 I=1,M SUM = 0.0D0 DO 170 J=1,M SUM = SUM + WM(J)*KERNEL(TM(I),TM(J))*R(J) 170 CONTINUE DO 180 J=1,N SUM = SUM + WN(J)*KERNEL(TM(I),TN(J))*DELN(J) 180 CONTINUE XM(I) = SUM + R(I) + XMZ(I) 190 CONTINUE RETURN END DOUBLE PRECISION FUNCTION NORM(X, Y, N, IFLAG) NOR 10 C IFLAG=0 CALCULATE THE MAXIMUM NORM OF X. C IFLAG=1 CALCULATE THE MAXIMUM NORM OF X-Y. DOUBLE PRECISION X, Y, DMAX1, DABS DIMENSION X(N), Y(N) IF (IFLAG.EQ.1) GO TO 20 C FIND THE NORM OF X. NORM = 0.0D0 DO 10 I=1,N NORM = DMAX1(NORM,DABS(X(I))) 10 CONTINUE RETURN C FIND THE NORM OF X-Y. 20 NORM = 0.0D0 DO 30 I=1,N NORM = DMAX1(NORM,DABS(X(I)-Y(I))) 30 CONTINUE RETURN END SUBROUTINE LINSYS(A, D, N, B, X, OPTION, DET, MACHIN) LIN 10 C SOLVE A*X=B, ORDER(A)=N, DIMENSION OF A=MACHIN. C OPTION=1 SOLVE A*X=B, LEAVE THE LU DECOMPOSITION A=L*U IN D C AND THE PIVOTS IN PIVOT. THE ANSWERS ARE LEFT IN X. C IT IS PERMISSABLE TO LET B=X AND D=A, BUT THEN THE C ORIGINAL CONTENTS OF A AND B ARE LOST. C OPTION=2 CALCULATE DECOMPOSITION A=L*U INCLUDING PIVOTS. SOLVE C A*X=B, AND THEN CALCULATE THE RESIDUAL AND ONE C CORRECTION. THE CORRECTIONS ARE LEFT IN R, THE NEW C VALUE OF X1 IN X, THE RELATIVE ERROR C NORM(X0-X1)/NORM(X1) C IN THE VARIABLE ERROR, AND THE RELATIVE RESIDUAL C NORM(RESIDUAL)/NORM(B) C IN THE VARIABLE RELRSD. THESE VALUES CAN BE OBTAINED C USING THE COMMON/LINEAR/ GIVEN BELOW. C OPTION=3 SAME AS OPTION=1, EXCEPT A=L*U IS KNOWN AND STORED C IN D. C OPTION=4 SAME AS OPTION=2, EXCEPT A=L*U IS KNOWN AND STORED C IN D. C THE DECOMPOSITION OF A INTO L*U USES SCALED PARTIAL PIVOTING IN C THE COLUMNS. FOR OPTIONS 1 AND 2, THE DETERMINANT OF A IS C CALCULATED AND STORED IN DET. IF DET=0, THEN THE ANSWERS ARE C NONSENSE, D AND X DO NOT CONTAIN USEFUL INFORMATION, AND AND A C AND B ARE LEFT UNDISTURBED (UNLESS D=A OR X=B). C LINSYS IS PRESENTLY LIMITED TO N .LE. 100. TO REMOVE THIS C RESTRICTION,CHANGE THE DIMENSION STATEMENT FOR THE ARRAYS R, C SCALE, AND PIVOT, WHICH ARE GIVEN IN COMMON/LINEAR/. INTEGER OPTION, PIVOT DOUBLE PRECISION NORMX, NORME, NORMB, NORMR, A, D, B, X, R, * SCALE, ERROR, RELRSD, DET, C, TEMP, DMAX1, DABS, SUM DIMENSION A(MACHIN,MACHIN), D(MACHIN,MACHIN), B(N), X(N) COMMON /LINEAR/ R(100), SCALE(100), ERROR, RELRSD, PIVOT(100) ISWIT = 1 IF (OPTION.GT.2) GO TO 110 C PRODUCE LU DECOMPOSITION OF A AND DET(A) DET = 0.0D0 DO 20 I=1,N DO 10 J=1,N D(I,J) = A(I,J) 10 CONTINUE 20 CONTINUE C PRODUCE SCALING FACTORS FOR SCALED PARTIAL PIVOTING. DO 40 I=1,N SCALE(I) = 0.0D0 DO 30 J=1,N SCALE(I) = DMAX1(SCALE(I),DABS(D(I,J))) 30 CONTINUE IF (SCALE(I).EQ.0.0D0) RETURN 40 CONTINUE DET = 1.0D0 NM1 = N - 1 C BEGIN GAUSSIAN ELIMINATION. DO 100 K=1,NM1 C SELECT PIVOT ROW. C = DABS(D(K,K))/SCALE(K) INDEX = K KP1 = K + 1 DO 50 I=KP1,N TEMP = DABS(D(I,K))/SCALE(I) IF (TEMP.LE.C) GO TO 50 INDEX = I C = TEMP 50 CONTINUE PIVOT(K) = INDEX IF (INDEX.EQ.K) GO TO 70 C SWITCH ROWS OF D. DO 60 J=K,N TEMP = D(K,J) D(K,J) = D(INDEX,J) D(INDEX,J) = TEMP 60 CONTINUE TEMP = SCALE(K) SCALE(K) = SCALE(INDEX) SCALE(INDEX) = TEMP DET = -DET 70 DET = DET*D(K,K) C ELIMINATE UNKNOWN =K FROM BELOW DIAGONAL IN COLUMN K. IF (DET.EQ.0.0D0) RETURN DO 90 I=KP1,N D(I,K) = D(I,K)/D(K,K) TEMP = D(I,K) DO 80 J=KP1,N D(I,J) = D(I,J) - TEMP*D(K,J) 80 CONTINUE 90 CONTINUE 100 CONTINUE DET = DET*D(N,N) IF (DET.EQ.0.0D0) RETURN C THE DECOMPOSITION A=P*L*U IS COMPLETE. C BEGIN SOLUTION OF LINEAR SYSTEM P*L*U*X=B C SET R AND X FOR SOLUTION OF A*X=B. 110 DO 120 I=1,N R(I) = B(I) X(I) = 0.0D0 120 CONTINUE GO TO 170 C SET R=B-A*X FOR COMPUTATION OF CORRECTION. 130 DO 150 I=1,N SUM = 0.0D0 DO 140 J=1,N SUM = SUM + A(I,J)*X(J) 140 CONTINUE R(I) = B(I) - SUM 150 CONTINUE NORMB = 0.0D0 NORMR = 0.0D0 DO 160 I=1,N NORMB = DMAX1(NORMB,DABS(B(I))) NORMR = DMAX1(NORMR,DABS(R(I))) 160 CONTINUE RELRSD = NORMR/NORMB ISWIT = 2 C SOLVE P*L*Z=R AND STORE IN R. 170 NM1 = N - 1 DO 200 K=1,NM1 IF (PIVOT(K).EQ.K) GO TO 180 KPIV = PIVOT(K) TEMP = R(K) R(K) = R(KPIV) R(KPIV) = TEMP 180 KP1 = K + 1 DO 190 I=KP1,N R(I) = R(I) - D(I,K)*R(K) 190 CONTINUE 200 CONTINUE C SOLVE U*E=R AND STORE IN R. ALSO LET X=X+E. R(N) = R(N)/D(N,N) GO TO (210, 220, 210, 220), OPTION 210 X(N) = R(N) GO TO 230 220 X(N) = X(N) + R(N) 230 DO 270 II=1,NM1 I = N - II SUM = 0.0D0 IP1 = I + 1 DO 240 J=IP1,N SUM = SUM + D(I,J)*R(J) 240 CONTINUE R(I) = (R(I)-SUM)/D(I,I) GO TO (250, 260, 250, 260), OPTION 250 X(I) = R(I) GO TO 270 260 X(I) = X(I) + R(I) 270 CONTINUE C SOLUTION OF LINEAR SYSTEM IS COMPLETE. RETURN FOR OPTION=1,3. GO TO (280, 290, 280, 290), OPTION 280 RETURN 290 GO TO (130, 300), ISWIT C CALCULATE ERRORS BASED ON CORRECTION. 300 NORMX = 0.0D0 NORME = 0.0D0 DO 310 I=1,N NORMX = DMAX1(NORMX,DABS(X(I))) NORME = DMAX1(NORME,DABS(R(I))) 310 CONTINUE ERROR = NORME/NORMX RETURN END SUBROUTINE IEGAUS(KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, IEG 10 * NUPPER, MUPPER, W, IER) C THE INTEGRAL EQUATION BEING SOLVED IS C B C X(S) - INT KERNEL(S,T)*X(T)*DT = RHFCN(S) C A C THE METHOD BEING USED IS BASED ON THE NYSTROM METHOD WITH C GAUSSIAN QUADRATURE, WITH AN ITERATIVE TECHNIQUE OF SOLUTION C FOR THE RESULTING LINEAR SYSTEM. C KERNEL THESE ARE DOUBLE PRECISION FUNCTIONS OF TWO AND ONE C RHFCN VARIABLES, RESPECTIVELY. THEY MUST BE DECLARED IN AN C EXTERNAL STATEMENT IN THE PROGRAM CALLING IEGAUS. C EP THE DESIRED ERROR. THE VARIABLE EP IS CHANGED ON C COMPLETION OF THE PROGRAM. SEE THE DISCUSSION OF IER C AND IFLAG FOR MORE INFORMATION. C IFLAG =0 EP IS INTERPRETED AS AN ABSOLUTE ERROR TOLERANCE. C =1 EP IS INTERPRETED AS A RELATIVE ERROR TOLERANCE. C X THE COMPUTED APPROXIMATE SOLUTION OF THE INTEGRAL C EQUATION, EVALUATED AT THE NODE POINTS IN T, IS C STORED IN X ON COMPLETION OF THE ROUTINE. THIS IS C TRUE IRREGARDLESS OF WHETHER OR NOT THE DESIRED ERROR C TOLERANCE WAS ATTAINED. C T CONTAINS THE NODE POINTS AT WHICH THE SOLUTION OF THE C INTEGRAL EQUATION IS DESIRED. SEE THE VARIABLE NT FOR C MORE INFORMATION. C NT IF NT=0, THEN T AND X WILL BE SET EQUAL TO THE FINAL C GAUSSIAN NODES AND THE CORRESPONDING SOLUTION VALUES, C AND NT WILL BE SET TO THE NUMBER OF THE SOLUTION C VALUES STORED IN X AND T. THE ARRAYS T AND X SHOULD C HAVE DIMENSION AT LEAST MUPPER, ASSIGNED IN THE C CALLING PROGRAM. C IF NT .GT. 0, THEN T CONTAINS NT USER SUPPLIED NODES C AT WHICH THE SOLUTION X IS DESIRED. C NUPPER AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM. C N IS THE ORDER OF A LINEAR SYSTEM WHICH IS BEING C USED TO ITERATIVELY SOLVE A LARGER LINEAR SYSTEM OF C ORDER M WHICH APPROXIMATES THE ABOVE INTEGRAL C EQUATION. C MUPPER AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM. C N AND M ARE ALWAYS POWERS OF TWO. C W TEMPORARY WORKING STORAGE FOR THE PROGRAM. IT MUST C CONTAIN AT LEAST 5*NU*NU+9*(NU+MU) POSITIONS, WITH C NU=NUPPER, MU=MUPPER. C IER =0 THIS ERROR COMPLETION CODE MEANS THE ROUTINE WAS C COMPLETED SATISFACTORILY. EP CONTAINS THE PREDICTED C ERROR. C =1 THE ERROR TEST WAS NOT SATISFIED. EP CONTAINS THE C PREDICTED ERROR. C =2 THE ERROR TEST WAS NOT SATISFIED. EP HAS BEEN SET C TO ZERO. C =3 THE ORIGINAL VALUE OF EP WAS TOO SMALL, DUE TO C POSSIBLE ILL-CONDITIONING PROBLEMS IN THE INTEGRAL C EQUATION. THE VALUE OF EP WAS RESET TO A MORE C REALISTIC VALUE, AND THAT TOLERANCE WAS ATTAINED. C =4 THE ERROR WAS SATISFACTORY AT THE GAUSSIAN NODE C POINTS (IER=0), BUT THE INTERPOLATION PROCESS(DUE TO C NT .GT. 0) MAY NOT PRESERVE THIS ACCURACY. CHECK THE C VALUE OF NORM(K) FOR POSSIBLE INDICATIONS THAT THE C INTEGRAL EQUATION MAY BE ALMOST FIRST KIND. SUCH C EQUATIONS ARE QUITE ILL-CONDITIONED. THE ERROR IN EP C IS THE PREDICTED ERROR FOR THE SOLUTION AT THE C GAUSSIAN NODE POINTS OF ORDER MFINAL. C =5 THE ANALOGUE OF IER=4, BUT WITH IER=1 AT THE C GAUSSIAN NODE POINTS. C =6 THE ANALOGUE OF IER=4, BUT WITH IER=3 AT THE C GAUSSIAN NODE POINTS. C IEGAUS IS PRESENTLY LIMITED TO NUPPER .LE. 100. THIS IS ENTIRELY C DUE TO A RESTRICTION IN THE PROGRAM LINSYS. TO REMOVE THIS C RESTRICTION, CHANGE THE DIMENSION STATEMENTS IN THE C COMMON/LINEAR/ STATEMENTS IN LINSYS AND THE SUBPROGRAM IEGS. DOUBLE PRECISION KERNEL, RHFCN, A, B, EP, X, T, W, CUTOFF, * ROOTRT, UNITRD, R1, R2, FINLEP, NORMK DIMENSION X(1), T(1), W(1) EXTERNAL KERNEL, RHFCN C ******************************************************************* C * COMMON /INFO/ R1, R2, FINLEP, NORMK, NFINAL, MFINAL C * C THE NUMBERS IN INFO GIVE ADDITIONAL INFORMATION ABOUT THE * C FUNCTIONING OF IEGAUS. R1 IS THE ITERATIVE RATE OF CONVERGENCE * C IN THE MOST RECENTLY COMPUTED LINEAR SYSTEM. R2 IS THE RATE OF * C CONVERGENCE OF THE GAUSSIAN QUADRATURE VARIANT OF THE NYSTROM * C METHOD. FINLEP IS THE FINAL VALUE OF EP USED AS THE DESIRED * C ERROR TOLERANCE. USUALLY FINLEP WILL EQUAL THE INPUT VALUE OF * C EP, UNLESS EP WAS MUCH TOO SMALL. NORMK IS AN APPROXIMATE * C VALUE FOR THE NORM OF THE INTEGRAL OPERATOR K, AND IT IS * C CALCULATED ONLY IF NT .GT. 0. * C NFINAL AND MFINAL ARE THE FINAL VALUES OF N AND M USED IN * C IEGAUS. IF NFINAL=MFINAL, THEN ITERATION WAS NOT INVOKED * C SUCCESSFULLY. * C * C ******************************************************************* C * DATA UNITRD /2.22D-16/ C * C UNITRD IS THE SMALLEST NUMBER U FOR WHICH * C 1+U .GT. 1 . * C UNITRD VARIES WITH THE COMPUTER AND WITH THE ARITHMETIC BEING * C USED. TO CHANGE TO ANOTHER ARITHMETIC, UNITRD MUST BE CHANGED. * C BUT UNITRD ALSO REFLECTS THE ACCURACY OF THE CONSTANTS USED * C IN SUBROUTINE WANDT. * C * C ******************************************************************* DATA CUTOFF /0.5D0/, ROOTRT /0.1D0/ C SET UP THE RELATIVE BASE ADDRESSES FOR THE VARIOUS ARRAYS INTO C WHICH W IS TO BE SPLIT. N = NUPPER M = MUPPER NSQ = N*N I1 = 1 I2 = I1 + NSQ I3 = I2 + NSQ I4 = I3 + NSQ/2 I5 = I4 + NSQ/2 I6 = I5 + M I7 = I6 + M I8 = I7 + N I9 = I8 + N I10 = I9 + M I11 = I10 + N I12 = I11 + M I13 = I12 + M I14 = I13 + M I15 = I14 + N I16 = I15 + M I17 = I16 + M I18 = I17 + N I19 = I18 + M I20 = I19 + 4*N I21 = I20 + NSQ NHALF = N/2 CALL IEGS(KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, NUPPER, * MUPPER, IER, CUTOFF, ROOTRT, UNITRD, NHALF, W(I1), W(I2), * W(I3), W(I4), W(I5), W(I6), W(I7), W(I8), W(I9), W(I10), * W(I11), W(I12), W(I13), W(I14), W(I15), W(I16), W(I17), * W(I18), W(I19), W(I20), W(I21)) RETURN END SUBROUTINE IEGS(KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, IEG 10 * NUP, MUP, IER, CUTOFF, ROOTRT, UNITRD, NHALF, LUFACT, KMM, * KMN, KNM, RHS, R, RH, DELN, TM, TN, XM, XMZ, WM, WN, OLDX, * SAVE, XN, SAVE2, ASIDE, ASIDE3, IMKNN) C THIS ROUTINE CONTROLS THE SOLUTION OF THE INTEGRAL EQUATION. DOUBLE PRECISION KERNEL, RHFCN, A, B, EP, X, T, CUTOFF, * ROOTRT, UNITRD, LUFACT, KMM, KMN, KNM, RHS, R, RH, DELN, TM, * TN, XM, XMZ, WM, WN, OLDX, SAVE, XN, SAVE2, ASIDE, ASIDE2, * ASIDE3, IMKNN, RESID, SCALE, ELINSY, RELRSD, XNORM, PASTC, * PASTRE, R1, R2, FINLEP, NORMK, DFLOAT, NORM, CONEW, RMIN, * R1RAT, COND, AVERR, EPS, DET, RELMIN, DMIN1, DMAX1, DSQRT, * ERROR, NUMR2, DENR2, TEMP2, RELERR, NUMR1, DENR1, RATE, * TEMP, RT, ESTERR, TEST INTEGER FLAG, OLDM DIMENSION X(1), T(1), LUFACT(NUP,NUP), KMM(NUP,NUP), * RHS(MUP), KNM(NHALF,NUP), KMN(NUP,NHALF), R(MUP), RH(NUP), * DELN(NUP), TM(MUP), TN(NUP), XM(MUP), XMZ(MUP), WM(MUP), * WN(NUP), OLDX(MUP), SAVE(MUP), XN(NUP), SAVE2(MUP), * ASIDE(NUP,4), ASIDE2(5), ASIDE3(NUP,NUP), IMKNN(NUP,NUP) COMMON /LINEAR/ RESID(100), SCALE(100), ELINSY, RELRSD, * IPIVOT(100) COMMON /INFO/ R1, R2, FINLEP, NORMK, NFINAL, MFINAL EXTERNAL KERNEL, RHFCN C INITIALIZATION LOOP = 1 N = 2 R2 = 0.5D0 M = 2*N R1RAT = ROOTRT COND = 1.0D0 PASTC = 1.0D0 PASTRE = 0.0D0 EPS = EP C STAGE A. DIRECT SOLUTION OF LINEAR SYSTEM (I-KN)*XN=RHS, WHILE C TRYING TO FIND A GOOD APPROXIMATE INVERSE TO IMPLEMENT C ITERATIVE METHOD OF SOLUTION. C CREATE THE NODES AND WEIGHTS TN(I) AND WN(I), I=1,...,N CALL WANDT(WN, TN, N, A, B) C SET UP MATRIX FOR (I-KN)*XN=RHFCN DO 20 J=1,N DO 10 I=1,N IMKNN(I,J) = -WN(J)*KERNEL(TN(I),TN(J)) 10 CONTINUE XMZ(J) = RHFCN(TN(J)) IMKNN(J,J) = IMKNN(J,J) + 1.0D0 20 CONTINUE GO TO 60 C THIS IS ENTRANCE FOR AN INCREASED VALUE OF N, USING PREVIOUSLY C STORED VALUES IN KMM TO DEFINE MATRIX FOR (I-KN)*XN=RHFCN WITH C NEW VALUE OF N. 30 DO 50 J=1,N DO 40 I=1,N IMKNN(I,J) = -KMM(I,J) 40 CONTINUE WN(J) = WM(J) TN(J) = TM(J) XMZ(J) = RHS(J) IMKNN(J,J) = IMKNN(J,J) + 1.0D0 50 CONTINUE C THIS IS THE ENTRANCE WHEN ITERATION IN STAGE B FAILS AND WE NEED C TO INCREASE N TO OBTAIN A BETTER ITERATIVE RATE. 60 CONTINUE C SOLVE (I-KN)*XN=RHFCN AT ALL TN(I).ALSO OBTAIN THE LU C DECOMPOSITION FOR LATER USE IN THE STAGE B ITERATIVE METHOD. C ******************************************************************* C * CALL LINSYS(IMKNN, LUFACT, N, XMZ, XN, 2, DET, NUP) C * C LINSYS IS A GENERAL LINEAR EQUATION SOLVER. IT HAS SPECIAL * C OPTIONS WHICH ARE USED IN THE FOLLOWING PROGRAM, AND THUS * C LINSYS SHOULD NOT BE REPLACED BY ANOTHER LINEAR EQUATION * C PROGRAM. LINSYS IS ALSO USED IN THE SUBROUTINE ITERT. * C * C ******************************************************************* COND = CONEW(COND,ELINSY,RELRSD,AVERR,PASTC,PASTRE) RELMIN = RMIN(N,N,COND,UNITRD,AVERR) IF (LOOP.EQ.1) GO TO 100 IF (LOOP.EQ.2) GO TO 80 C SET UP APPROXIMATE RATE OF CONVERGENCE OF SOLUTIONS XN TO TRUE C SOLUTION X. ALSO SET UP DESIRED RATIO FOR ITERATIVE METHOD. NUMR2 = NORM(XN,OLDX,N,1) R2 = DMIN1(0.5D0,DMAX1(NUMR2/DENR2,1.0D-4)) R1RAT = DMIN1(ROOTRT,DSQRT(R2)) C CHECK FOR ERROR IN XN USING TEST INVOLVING R2 AND OLDX,ACCORDING C TO THEORY FOR ASYMPTOTIC ERROR BOUNDS. MODIFY ERROR IF IT IS C OUTSIDE PRECISION RANGE OF COMPUTER, POSSIBLY DUE TO C ILL-CONDITIONING. 70 ERROR = (R2/(1.0D0-R2))*NUMR2 XNORM = NORM(XN,XN,N,0) RELERR = ERROR/XNORM IF (IFLAG.EQ.0) EPS = DMAX1(EP,XNORM*RELMIN) IF (IFLAG.EQ.1) EPS = DMAX1(EP,RELMIN) IF (IFLAG.EQ.1) ERROR = DMAX1(RELERR,RELMIN) IF ((IFLAG.EQ.0) .AND. (RELERR.LT.RELMIN)) ERROR = * RELMIN*XNORM IF (ERROR.LE.EPS) GO TO 90 DENR2 = NUMR2 GO TO 100 C ENTRANCE FOR LOOP=2. 80 NUMR2 = NORM(XN,OLDX,N,1) DENR2 = 0.0D0 GO TO 70 C EXIT FOR SUCCESSFUL RETURN. ITERATION WAS NOT NECESSARY. 90 CALL LEAVE(0, N, N, XN, TN, WN, ERROR, KERNEL, RHFCN, EP, * IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, * XNORM) RETURN C ATTEMPT TO SOLVE (I-KM)*XM=RHFCN ITERATIVELY, CHECKING TO SEE IF C THE RATE OF CONVERGENCE IS SUFFICIENTLY FAST SO AS TO ENTER C STAGE B. C CALCULATE TM(I) AND WM(I), I=1,...,M. 100 CALL WANDT(WM, TM, M, A, B) FLAG = 0 C CALCULATE INITIAL GUESS XMZ FOR ITERATION METHOD. CALL INTERP(TM, WM, XMZ, M, TN, WN, XN, N, KERNEL, RHFCN, * RHS, KMN, NHALF, NUP) DO 110 I=1,M OLDX(I) = XMZ(I) 110 CONTINUE C CALCULATE FIRST ITERATE. CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) COND = CONEW(COND,ELINSY,RELRSD,AVERR,PASTC,PASTRE) DENR1 = NORM(XM,XMZ,M,1) FLAG = 1 DO 120 I=1,M XMZ(I) = XM(I) 120 CONTINUE C CALCULATE SECOND ITERATE. CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) COND = CONEW(COND,ELINSY,RELRSD,AVERR,PASTC,PASTRE) NUMR1 = NORM(XM,XMZ,M,1) C CHECK ON THE SPEED OF CONVERGENCE OF ITERATIVE METHOD. IF IT IS C SUFFICIENTLY RAPID, THEN FIX N AND GO TO STAGE B. R1 = NUMR1/DENR1 RATE = R1 IF (M.GT.NUP) GO TO 170 IF (R1.LE.R1RAT) GO TO 130 C THE ITERATION DID NOT WORK WELL ENOUGH, AND STAGE A IS TO BE C REPEATED. RE-INITIALIZE FOR SOLVING (I-KN)*XN=RHFCN AGAIN C WITH A LARGER N. N = M LOOP = LOOP + 1 M = 2*N GO TO 30 C THE ITERATIVE RATE IS SUFFICIENTLY RAPID, AND CONTROL WILL GO TO C STAGE B. SAVE INFORMATION IN CASE STAGE B ABORTS AT A LARGER C VALUE OF M AND STAGE A HAS TO BE RETURNED TO. 130 DO 140 I=1,M ASIDE(I,1) = OLDX(I) ASIDE(I,2) = WM(I) ASIDE(I,3) = TM(I) ASIDE(I,4) = RHS(I) 140 CONTINUE ASIDE2(1) = LOOP ASIDE2(3) = R2 ASIDE2(4) = DENR2 ASIDE2(5) = R1RAT DO 160 J=1,M DO 150 I=1,M ASIDE3(I,J) = KMM(I,J) 150 CONTINUE 160 CONTINUE C STAGE B. ITERATIVE METHOD OF SOLUTION OF (I-KM)*XM=RHS. 170 OLDM = N ASIDE2(2) = M IF (R1.LE.CUTOFF) GO TO 200 C THE ITERATES ARE CONVERGING VERY SLOWLY OR NOT AT ALL. THUS C RETURN WITHOUT FURTHER ATTEMPTS TO LESSEN THE ERROR. IF (LOOP.NE.1) GO TO 190 180 CALL LEAVE(2, N, N, XN, TN, WN, 0.0D0, KERNEL, RHFCN, EP, * IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, * XNORM) RETURN 190 CALL LEAVE(1, N, N, XN, TN, WN, ERROR, KERNEL, RHFCN, EP, * IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, * XNORM) RETURN C TEST TO SEE IF THE CURRENT ITERATE XM IS SUFFICIENTLY ACCURATE C COMPARED TO THE TRUE XM. 200 RATE = R1*RATE TEMP = NORM(XM,OLDX,M,1) IF (LOOP.EQ.1) TEMP2 = 0.5D0 IF (LOOP.GT.1) TEMP2 = TEMP/DENR2 RT = DMIN1(0.01D0,DMAX1(TEMP2,0.0001D0))/2.0D0 XNORM = NORM(XM,XM,M,0) ESTERR = (RT/(1.0D0-RT))*TEMP/XNORM IF (ESTERR.LT.RELMIN) ESTERR = RELMIN ESTERR = ESTERR*XNORM TEST = ((1.0D0-R1)/R1)*ESTERR IF (NUMR1.LE.TEST) GO TO 260 C ITERATE NOT SUFFICIENTLY ACCURATE. INITIALIZE FOR COMPUTATION C OF ANOTHER ITERATE. DENR1 = NUMR1 DO 210 I=1,M XMZ(I) = XM(I) 210 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) COND = CONEW(COND,ELINSY,RELRSD,AVERR,PASTC,PASTRE) NUMR1 = NORM(XM,XMZ,M,1) R1 = NUMR1/DENR1 IF (R1.LE.CUTOFF) GO TO 200 C THIS IS ENTRANCE FOR CASE WHERE ITERATION FAILS IN STAGE B. C PARAMETERS MUST BE RESET FOR A RETURN TO STAGE A OR FOR AN C ABORTIVE EXIT IF N CANNOT BE INCREASED ANY FURTHER. 220 MNEW = ASIDE2(2) IF (MNEW.GT.NUP) GO TO 250 N = MNEW DO 240 J=1,N DO 230 I=1,N IMKNN(I,J) = -ASIDE3(I,J) 230 CONTINUE OLDX(J) = ASIDE(J,1) WN(J) = ASIDE(J,2) TN(J) = ASIDE(J,3) XMZ(J) = ASIDE(J,4) IMKNN(J,J) = IMKNN(J,J) + 1.0D0 240 CONTINUE M = 2*N LOOP = ASIDE2(1) + 1.0D0 R2 = ASIDE2(3) DENR2 = ASIDE2(4) R1RAT = ASIDE2(5) GO TO 60 C ABORTIVE EXIT FROM STAGE B. N CANNOT BE INCREASED FURTHER, AND C R1 IS NOT SUFFICIENTLY SMALL. 250 IF (LOOP.EQ.1) GO TO 180 CALL WANDT(WM, TM, OLDM, A, B) CALL LEAVE(1, N, OLDM, SAVE, TM, WM, ERROR, KERNEL, RHFCN, * EP, IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, * XMZ, KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, * NHALF, XNORM) RETURN C AN ACCURATE VALUE OF XM HAS BEEN OBTAINED. R2 IS TO BE TESTED AS C TO WHETHER IT SHOULD BE RESET. THEN CHECK ERROR IN XM COMPARED C WITH THE TRUE SOLUTION X. 260 IF (LOOP.EQ.1) GO TO 290 NUMR2 = TEMP R2 = DMAX1(1.0D-4,RATE,DMIN1(NUMR2/DENR2,0.5D0)) DENR2 = NUMR2 270 ERROR = (R2/(1.0D0-R2))*TEMP XNORM = NORM(XM,XM,M,0) RELERR = ERROR/XNORM RELMIN = RMIN(N,M,COND,UNITRD,AVERR) IF (IFLAG.EQ.0) EPS = DMAX1(EP,XNORM*RELMIN) IF (IFLAG.EQ.1) EPS = DMAX1(EP,RELMIN) IF (IFLAG.EQ.1) ERROR = DMAX1(RELERR,RELMIN) IF ((IFLAG.EQ.0) .AND. (RELERR.LT.RELMIN)) ERROR = * RELMIN*XNORM IF (ERROR.GT.EPS) GO TO 280 CALL LEAVE(0, N, M, XM, TM, WM, ERROR, KERNEL, RHFCN, EP, * IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, * XNORM) RETURN 280 MNEW = 2*M IF (MNEW.LE.MUP) GO TO 300 C M CANNOT BE INCREASED ANY FURTHER. CALL LEAVE(1, N, M, XM, TM, WM, ERROR, KERNEL, RHFCN, EP, * IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, * XNORM) RETURN 290 DENR2 = TEMP LOOP = 2 GO TO 270 C ERROR NOT SUFFICIENTLY SMALL. M IS INCREASED AND TWO MORE C MORE ITERATES ARE COMPUTED WITH THE NEW M. 300 OLDM = M M = MNEW DO 310 I=1,OLDM SAVE2(I) = WM(I) SAVE(I) = TM(I) 310 CONTINUE CALL WANDT(WM, TM, M, A, B) FLAG = 0 CALL INTERP(TM, WM, XMZ, M, SAVE, SAVE2, XM, OLDM, KERNEL, * RHFCN, RHS, KMN, NHALF, NUP) DO 320 I=1,OLDM SAVE(I) = XM(I) 320 CONTINUE DO 330 I=1,M OLDX(I) = XMZ(I) 330 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) COND = CONEW(COND,ELINSY,RELRSD,AVERR,PASTC,PASTRE) DENR1 = NORM(XM,XMZ,M,1) FLAG = 1 DO 340 I=1,M XMZ(I) = XM(I) 340 CONTINUE CALL ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, FLAG) COND = CONEW(COND,ELINSY,RELRSD,AVERR,PASTC,PASTRE) NUMR1 = NORM(XM,XMZ,M,1) R1 = NUMR1/DENR1 RATE = R1 IF (R1.LE.CUTOFF) GO TO 200 GO TO 220 END SUBROUTINE WANDT(WV, TV, N, A, B) WAN 10 C INTEGRATION WEIGHTS AND NODES ARE TO BE CALCULATED AND STORED IN C WV AND TV, RESPECTIVELY. N IS ASSUMED TO BE A POWER OF TWO. IF C 2 .LE. N .LE. 256, THEN GAUSSIAN QUADRATURE IS USED. IF N .GT. C 256, THEN THE INTERVAL (A,B) IS DIVIDED N/256 TIMES AND THE 256 C POINT FORMULA IS APPLIED TO EACH SUBINTERVAL. DOUBLE PRECISION WV, TV, A, B, W, T, FLOOP, H, SCALE, AL, BL, * S, R, FL DIMENSION WV(N), TV(N) DIMENSION W(255), T(255) DATA T(1), T(2), T(3), T(4), T(5), T(6), T(7), T(8), T(9), * T(10), T(11), T(12), T(13), T(14), T(15) * /.577350269189626D0,.861136311594053D0,.339981043584856D0, * .960289856497536D0,.796666477413627D0,.525532409916329D0, * .183434642495650D0,.989400934991650D0,.944575023073233D0, * .865631202387832D0,.755404408355003D0,.617876244402644D0, * .458016777657227D0,.281603550779259D0,.950125098376374D-1/ DATA T(16), T(17), T(18), T(19), T(20), T(21), T(22), T(23), * T(24), T(25), T(26), T(27), T(28), T(29), T(30), T(31) * /.997263861849482D0,.985611511545268D0,.964762255587506D0, * .934906075937740D0,.896321155766052D0,.849367613732570D0, * .794483795967942D0,.732182118740290D0,.663044266930215D0, * .587715757240762D0,.506899908932229D0,.421351276130635D0, * .331868602282128D0,.239287362252137D0,.144471961582796D0, * .483076656877383D-1/ DATA T(32), T(33), T(34), T(35), T(36), T(37), T(38), T(39), * T(40), T(41), T(42), T(43), T(44), T(45), T(46), T(47) * /.999305041735772D0,.996340116771955D0,.991013371476744D0, * .983336253884626D0,.973326827789911D0,.961008799652054D0, * .946411374858403D0,.929569172131940D0,.910522137078503D0, * .889315445995114D0,.865999398154093D0,.840629296252580D0, * .813265315122798D0,.783972358943341D0,.752819907260532D0, * .719881850171611D0/ DATA T(48), T(49), T(50), T(51), T(52), T(53), T(54), T(55), * T(56), T(57), T(58), T(59), T(60), T(61), T(62), T(63) * /.685236313054233D0,.648965471254657D0,.611155355172393D0, * .571895646202634D0,.531279464019895D0,.489403145707053D0, * .446366017253464D0,.402270157963992D0,.357220158337668D0, * .311322871990211D0,.264687162208767D0,.217423643740007D0, * .169644420423993D0,.121462819296121D0,.729931217877990D-1, * .243502926634244D-1/ DATA T(64), T(65), T(66), T(67), T(68), T(69), T(70), T(71), * T(72), T(73), T(74), T(75), T(76), T(77), T(78), T(79) * /.999824887947132D0,.999077459977376D0,.997733248625514D0, * .995792758534981D0,.993257112900213D0,.990127818491734D0, * .986406742724586D0,.982096108435719D0,.977198491463907D0, * .971716818747137D0,.965654366431965D0,.959014757853700D0, * .951801961341264D0,.944020287830220D0,.935674388277916D0, * .926769250878948D0/ DATA T(80), T(81), T(82), T(83), T(84), T(85), T(86), T(87), * T(88), T(89), T(90), T(91), T(92), T(93), T(94), T(95), * T(96), T(97) /.917310198080961D0,.907302883401757D0, * .896753288049158D0,.885667717345397D0,.874052796958032D0, * .861915468939548D0,.849262987577969D0,.836102915060907D0, * .822443116955644D0,.808291757507914D0,.793657294762193D0, * .778548475506412D0,.762974330044095D0,.746944166797062D0, * .730467566741909D0,.713554377683587D0,.696214708369514D0, * .678458922447719D0/ DATA T(98), T(99), T(100), T(101), T(102), T(103), T(104), * T(105), T(106), T(107), T(108), T(109), T(110), T(111), * T(112), T(113), T(114) /.660297632272646D0, * .641741692562308D0,.622802193910585D0,.603490456158549D0, * .583818021628763D0,.563796648226618D0,.543438302412810D0, * .522755152051175D0,.501759559136144D0,.480464072404172D0, * .458881419833552D0,.437024501037104D0,.414906379552275D0, * .392540275033267D0,.369939555349859D0,.347117728597636D0, * .324088435024413D0/ DATA T(115), T(116), T(117), T(118), T(119), T(120), T(121), * T(122), T(123), T(124), T(125), T(126), T(127) * /.300865438877677D0,.277462620177904D0,.253893966422694D0, * .230173564226660D0,.206315590902079D0,.182334305985337D0, * .158244042714225D0,.134059199461188D0,.109794231127644D0, * .854636405045155D-1,.610819696041396D-1,.366637909687335D-1, * .122236989606158D-1/ DATA T(128), T(129), T(130), T(131), T(132), T(133), T(134), * T(135), T(136), T(137), T(138), T(139), T(140), T(141), * T(142) /.999956050018992D0,.999768437409263D0, * .999430937466261D0,.998943525843409D0,.998306266473006D0, * .997519252756721D0,.996582602023382D0,.995496454481096D0, * .994260972922410D0,.992876342608822D0,.991342771207583D0, * .989660488745065D0,.987829747564861D0,.985850822286126D0, * .983724009760315D0/ DATA T(143), T(144), T(145), T(146), T(147), T(148), T(149), * T(150), T(151), T(152), T(153), T(154), T(155), T(156), * T(157) /.981449629025464D0,.979028021257622D0, * .976459549719234D0,.973744599704370D0,.970883578480743D0, * .967876915228489D0,.964725060975706D0,.961428488530732D0, * .957987692411178D0,.954403188769716D0,.950675515316628D0, * .946805231239127D0,.942792917117462D0,.938639174837814D0, * .934344627502003D0/ DATA T(158), T(159), T(160), T(161), T(162), T(163), T(164), * T(165), T(166), T(167), T(168), T(169), T(170), T(171), * T(172) /.929909919334006D0,.925335715583316D0, * .920622702425146D0,.915771586857490D0,.910783096595065D0, * .905657979960145D0,.900397005770304D0,.895000963223085D0, * .889470661777611D0,.883806931033158D0,.878010620604707D0, * .872082599995488D0,.866023758466555D0,.859835004903376D0, * .853517267679503D0/ DATA T(173), T(174), T(175), T(176), T(177), T(178), T(179), * T(180), T(181), T(182), T(183), T(184), T(185), T(186), * T(187) /.847071494517296D0,.840498652345763D0, * .833799727155505D0,.826975723850813D0,.820027666098917D0, * .812956596176432D0,.805763574812999D0,.798449681032171D0, * .791016011989546D0,.783463682808184D0,.775793826411326D0, * .768007593352446D0,.760106151642655D0,.752090686575492D0, * .743962400549112D0/ DATA T(188), T(189), T(190), T(191), T(192), T(193), T(194), * T(195), T(196), T(197), T(198), T(199), T(200), T(201), * T(202) /.735722512885918D0,.727372259649652D0, * .718912893459971D0,.710345683304543D0,.701671914348685D0, * .692892887742577D0,.684009920426076D0,.675024344931163D0, * .665937509182049D0,.656750776292973D0,.647465524363725D0, * .638083146272911D0,.628605049469015D0,.619032655759261D0, * .609367401096334D0/ DATA T(203), T(204), T(205), T(206), T(207), T(208), T(209), * T(210), T(211), T(212), T(213), T(214), T(215), T(216), * T(217) /.599610735362968D0,.589764122154454D0, * .579829038559083D0,.569806974936569D0,.559699434694481D0, * .549507934062719D0,.539234001866059D0,.528879179294822D0, * .518445019673674D0,.507933088228616D0,.497344961852181D0, * .486682228866890D0,.475946488786983D0,.465139352078479D0, * .454262439917590D0/ DATA T(218), T(219), T(220), T(221), T(222), T(223), T(224), * T(225), T(226), T(227), T(228), T(229), T(230), T(231), * T(232) /.443317383947527D0,.432305826033741D0, * .421229418017624D0,.410089821468717D0,.398888707435459D0, * .387627756194516D0,.376308656998716D0,.364933107823654D0, * .353502815112970D0,.342019493522372D0,.330484865662417D0, * .318900661840106D0,.307268619799319D0,.295590484460136D0, * .283868007657082D0/ DATA T(233), T(234), T(235), T(236), T(237), T(238), T(239), * T(240), T(241), T(242), T(243), T(244), T(245), T(246), * T(247) /.272102947876337D0,.260297069991943D0, * .248452145001057D0,.236569949758284D0,.224652266709132D0, * .212700883622626D0,.200717593323127D0,.188704193421389D0, * .176662486044902D0,.164594277567554D0,.152501378338656D0, * .140385602411376D0,.128248767270607D0,.116092693560333D0, * .103919204810509D0/ DATA T(248), T(249), T(250), T(251), T(252), T(253), T(254), * T(255) /.917301271635196D-1,.795272891002330D-1, * .673125211657164D-1,.550876556946340D-1,.428545265363791D-1, * .306149687799790D-1,.183708184788137D-1,.612391237518953D-2/ DATA W(1), W(2), W(3), W(4), W(5), W(6), W(7), W(8), W(9), * W(10), W(11), W(12), W(13), W(14), W(15) * /1.0D0,.347854845137454D0,.652145154862546D0, * .101228536290376D0,.222381034453374D0,.313706645877887D0, * .362683783378362D0,.271524594117541D-1,.622535239386479D-1, * .951585116824928D-1,.124628971255534D0,.149595988816577D0, * .169156519395003D0,.182603415044924D0,.189450610455068D0/ DATA W(16), W(17), W(18), W(19), W(20), W(21), W(22), W(23), * W(24), W(25), W(26), W(27), W(28), W(29), W(30), W(31) * /.701861000947010D-2,.162743947309057D-1,.253920653092621D-1, * .342738629130214D-1,.428358980222267D-1,.509980592623762D-1, * .586840934785355D-1,.658222227763618D-1,.723457941088485D-1, * .781938957870703D-1,.833119242269468D-1,.876520930044038D-1, * .911738786957639D-1,.938443990808046D-1,.956387200792749D-1, * .965400885147278D-1/ DATA W(32), W(33), W(34), W(35), W(36), W(37), W(38), W(39), * W(40), W(41), W(42), W(43), W(44), W(45), W(46), W(47) * /.178328072169643D-2,.414703326056247D-2,.650445796897836D-2, * .884675982636395D-2,.111681394601311D-1,.134630478967186D-1, * .157260304760247D-1,.179517157756973D-1,.201348231535302D-1, * .222701738083833D-1,.243527025687109D-1,.263774697150547D-1, * .283396726142595D-1,.302346570724025D-1,.320579283548516D-1, * .338051618371416D-1/ DATA W(48), W(49), W(50), W(51), W(52), W(53), W(54), W(55), * W(56), W(57), W(58), W(59), W(60), W(61), W(62), W(63) * /.354722132568824D-1,.370551285402400D-1,.385501531786156D-1, * .399537411327203D-1,.412625632426235D-1,.424735151236536D-1, * .435837245293235D-1,.445905581637566D-1,.454916279274181D-1, * .462847965813144D-1,.469681828162100D-1,.475401657148303D-1, * .479993885964583D-1,.483447622348030D-1,.485754674415034D-1, * .486909570091397D-1/ DATA W(64), W(65), W(66), W(67), W(68), W(69), W(70), W(71), * W(72), W(73), W(74), W(75), W(76), W(77), W(78), W(79) * /.449380960292090D-3,.104581267934035D-2,.164250301866903D-2, * .223828843096262D-2,.283275147145799D-2,.342552604091022D-2, * .401625498373864D-2,.460458425670296D-2,.519016183267633D-2, * .577263754286570D-2,.635166316170719D-2,.692689256689881D-2, * .749798192563473D-2,.806458989048606D-2,.862637779861675D-2, * .918300987166087D-2/ DATA W(80), W(81), W(82), W(83), W(84), W(85), W(86), W(87), * W(88), W(89), W(90), W(91), W(92), W(93), W(94), W(95) * /.973415341500681D-2,.102794790158322D-1,.108186607395031D-1, * .113513763240804D-1,.118773073727403D-1,.123961395439509D-1, * .129075627392673D-1,.134112712886163D-1,.139069641329520D-1, * .143943450041668D-1,.148731226021473D-1,.153430107688651D-1, * .158037286593993D-1,.162550009097852D-1,.166965578015892D-1, * .171281354231114D-1/ DATA W(96), W(97), W(98), W(99), W(100), W(101), W(102), * W(103), W(104), W(105), W(106), W(107), W(108), W(109), * W(110), W(111), W(112) /.175494758271177D-1, * .179603271850087D-1,.183604439373313D-1,.187495869405447D-1, * .191275236099509D-1,.194940280587066D-1,.198488812328309D-1, * .201918710421300D-1,.205227924869601D-1,.208414477807511D-1, * .211476464682213D-1,.214412055392085D-1,.217219495380521D-1, * .219897106684605D-1,.222443288937998D-1,.224856520327450D-1, * .227135358502365D-1/ DATA W(113), W(114), W(115), W(116), W(117), W(118), W(119), * W(120), W(121), W(122), W(123), W(124), W(125), W(126), * W(127) /.229278441436868D-1,.231284488243870D-1, * .233152299940628D-1,.234880760165359D-1,.236468835844476D-1, * .237915577810034D-1,.239220121367035D-1,.240381686810241D-1, * .241399579890193D-1,.242273192228152D-1,.243002001679719D-1, * .243585572646906D-1,.244023556338496D-1,.244315690978500D-1, * .244461801962625D-1/ DATA W(128), W(129), W(130), W(131), W(132), W(133), W(134), * W(135), W(136), W(137), W(138), W(139), W(140), W(141), * W(142) /.112789017822272D-3,.262534944296446D-3, * .412463254426176D-3,.562348954031410D-3,.712154163473321D-3, * .861853701420089D-3,.101142439320844D-2,.116084355756772D-2, * .131008868190250D-2,.145913733331073D-2,.160796713074933D-2, * .175655573633073D-2,.190488085349972D-2,.205292022796614D-2, * .220065164983991D-2/ DATA W(143), W(144), W(145), W(146), W(147), W(148), W(149), * W(150), W(151), W(152), W(153), W(154), W(155), W(156), * W(157) /.234805295632731D-2,.249510203470371D-2, * .264177682542749D-2,.278805532532771D-2,.293391559082972D-2, * .307933574119934D-2,.322429396179420D-2,.336876850731555D-2, * .351273770505631D-2,.365617995814250D-2,.379907374876626D-2, * .394139764140883D-2,.408313028605267D-2,.422425042138154D-2, * .436473687796806D-2/ DATA W(158), W(159), W(160), W(161), W(162), W(163), W(164), * W(165), W(166), W(167), W(168), W(169), W(170), W(171), * W(172) /.450456858144790D-2,.464372455568006D-2, * .478218392589269D-2,.491992592181387D-2,.505692988078684D-2, * .519317525086928D-2,.532864159391593D-2,.546330858864431D-2, * .559715603368291D-2,.573016385060144D-2,.586231208692265D-2, * .599358091911534D-2,.612395065556793D-2,.625340173954240D-2, * .638191475210788D-2/ DATA W(173), W(174), W(175), W(176), W(177), W(178), W(179), * W(180), W(181), W(182), W(183), W(184), W(185), W(186), * W(187) /.650947041505366D-2,.663604959378107D-2, * .676163330017380D-2,.688620269544632D-2,.700973909296982D-2, * .713222396107539D-2,.725363892583391D-2,.737396577381235D-2, * .749318645480588D-2,.761128308454566D-2,.772823794738156D-2, * .784403349893971D-2,.795865236875435D-2,.807207736287350D-2, * .818429146643827D-2/ DATA W(188), W(189), W(190), W(191), W(192), W(193), W(194), * W(195), W(196), W(197), W(198), W(199), W(200), W(201), * W(202) /.829527784623523D-2,.840501985322154D-2, * .851350102502249D-2,.862070508840101D-2,.872661596169881D-2, * .883121775724875D-2,.893449478375821D-2,.903643154866287D-2, * .913701276045081D-2,.923622333095630D-2,.933404837762327D-2, * .943047322573775D-2,.952548341062928D-2,.961906467984073D-2, * .971120299526628D-2/ DATA W(203), W(204), W(205), W(206), W(207), W(208), W(209), * W(210), W(211), W(212), W(213), W(214), W(215), W(216), * W(217) /.980188453525733D-2,.989109569669583D-2, * .997882309703491D-2,.100650535763064D-1,.101497741990949D-1, * .102329722564782D-1,.103146352679340D-1,.103947509832117D-1, * .104733073841704D-1,.105502926865815D-1,.106256953418966D-1, * .106995040389798D-1,.107717077058046D-1,.108422955111148D-1, * .109112568660490D-1/ DATA W(218), W(219), W(220), W(221), W(222), W(223), W(224), * W(225), W(226), W(227), W(228), W(229), W(230), W(231), * W(232) /.109785814257296D-1,.110442590908139D-1, * .111082800090098D-1,.111706345765534D-1,.112313134396497D-1, * .112903074958755D-1,.113476078955455D-1,.114032060430392D-1, * .114570935980906D-1,.115092624770395D-1,.115597048540436D-1, * .116084131622531D-1,.116553800949452D-1,.117005986066207D-1, * .117440619140606D-1/ DATA W(233), W(234), W(235), W(236), W(237), W(238), W(239), * W(240), W(241), W(242), W(243), W(244), W(245), W(246), * W(247) /.117857634973434D-1,.118256971008240D-1, * .118638567340711D-1,.119002366727665D-1,.119348314595636D-1, * .119676359049059D-1,.119986450878058D-1,.120278543565826D-1, * .120552593295601D-1,.120808558957245D-1,.121046402153405D-1, * .121266087205273D-1,.121467581157945D-1,.121650853785355D-1, * .121815877594818D-1/ DATA W(248), W(249), W(250), W(251), W(252), W(253), W(254), * W(255) /.121962627831147D-1,.122091082480372D-1, * .122201222273040D-1,.122293030687103D-1,.122366493950402D-1, * .122421601042728D-1,.122458343697479D-1,.122476716402898D-1/ LOOP = MAX0(1,N/256) FLOOP = LOOP H = (B-A)/FLOOP SCALE = H/2.0D0 M = MIN0(128,N/2) MT = 2*M NPLACE = M - 1 DO 20 L=1,LOOP FL = L AL = A + (FL-1.0D0)*H BL = A + FL*H K = 256*(L-1) DO 10 I=1,M NPI = NPLACE + I S = T(NPI) R = W(NPI)*SCALE I1 = K + I I2 = K + MT + 1 - I TV(I1) = (AL*(1.D0+S)+(1.0D0-S)*BL)/2.D0 TV(I2) = (AL*(1.D0-S)+(1.0D0+S)*BL)/2.D0 WV(I1) = R WV(I2) = R 10 CONTINUE 20 CONTINUE RETURN END DOUBLE PRECISION FUNCTION RMIN(N, M, COND, UNITRD, AVERR) RMI 10 C FOR A LINEAR SYSTEM (I-KMM)*XM=RHFCN OF ORDER M, THIS IS THE C VALUE OF RELMIN USED IN IEGS. THE VARIABLE UNITRD IS DEFINED IN C IEGAUS, AND THE VARIABLES COND AND AVERR ARE DEFINED IN IEGS C USING CONEW. C IT IS UNLIKELY THAT A SOLUTION X CAN BE FOUND FOR THE ORIGINAL C INTEGRAL EQUATION WITH A SMALLER RELATIVE ERROR THAN RMIN. DOUBLE PRECISION FLOAT1, FLOAT2, COND, AVERR, UNITRD, DMAX1 FLOAT1 = M FLOAT2 = M/N RMIN = DMAX1((FLOAT1**1.5D0)*COND*UNITRD,(FLOAT2**1.5D0)* * AVERR) RETURN END DOUBLE PRECISION FUNCTION CONEW(COND, ELINSY, RELRSD, AVERR, CON 10 * PASTC, PASTRE) C THIS IS USED IN UPDATING THE VALUE OF THE CONDITION NUMBER C IN IEGS. DOUBLE PRECISION COND, ELINSY, RELRSD, AVERR, PASTC, PASTRE, * C, DSQRT, DMAX1 AVERR = DSQRT(ELINSY*PASTRE) PASTRE = ELINSY IF (RELRSD.EQ.0.0D0) GO TO 10 C = DMAX1(1.0D0,ELINSY/RELRSD) CONEW = DSQRT(C*PASTC) PASTC = C RETURN 10 CONEW = COND RETURN END SUBROUTINE LEAVE(IERSET, NF, MF, XV, TV, WV, ERROR, KERNEL, LEA 10 * RHFCN, EP, IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, * WM, XM, XMZ, KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, * NUP, NHALF, XNORM) C THIS ROUTINE SETS ALL NECESSARY PARAMETERS FOR LEAVING IEGAUS. C IF NT .GT. 0, IT ALSO PERFORMS THE NECESSARY NYSTROM C INTERPOLATION AT THE NODES GIVEN IN T. DOUBLE PRECISION KERNEL, RHFCN, EP, X, T, XV, TV, WV, ERROR, * EPS, ELINSY, TN, WN, TM, WM, XM, XMZ, KMM, KMN, KNM, RHS, * IMKNN, LUFACT, R, RH, DELN, NORMK, R1, R2, FINLEP, SAVEP, * SUM, DERROR, NUMR1, NORM, TEMP, DABS, DMAX1, XNORM DIMENSION X(1), T(1), XV(MF), TV(MF), WV(MF), TN(NF), WN(NF), * WM(MF), XM(MF), XMZ(MF), KMM(NUP,NUP), KMN(NUP,NHALF), * KNM(NHALF,NUP), RHS(MF), IMKNN(NUP,NUP), LUFACT(NUP,NUP), * R(MF), RH(NF), TM(MF), DELN(NF) COMMON /INFO/ R1, R2, FINLEP, NORMK, NFINAL, MFINAL EXTERNAL KERNEL, RHFCN C SET ERROR PARAMETERS FOR RETURN. NORMK = 0.0D0 NFINAL = NF MFINAL = MF FINLEP = EPS IF ((EPS.GT.EP) .AND. (ERROR.LE.EPS)) GO TO 10 IER = IERSET EP = ERROR IF (NT.EQ.0) GO TO 20 GO TO 40 10 IER = 3 C SINCE EPS IS THE SMALLEST ERROR POSSIBLE, SET EP=EPS FOR THE C RETURN ERROR ESTIMATE. EP = EPS IF (NT.GT.0) GO TO 40 C NO NYSTROM INTERPOLATION IS DESIRED. RETURN THE VALUES AT THE C GAUSSIAN NODE POINTS. 20 DO 30 I=1,MF X(I) = XV(I) T(I) = TV(I) 30 CONTINUE NT = MF RETURN C CALCULATE NORM(K). 40 SAVEP = EP DO 50 I=1,NF IMKNN(I,I) = IMKNN(I,I) - 1.0D0 50 CONTINUE NORMK = 0.0D0 DO 70 I=1,NF SUM = 0.0D0 DO 60 J=1,NF SUM = SUM + DABS(IMKNN(I,J)) 60 CONTINUE NORMK = DMAX1(NORMK,SUM) 70 CONTINUE DO 80 I=1,NF IMKNN(I,I) = IMKNN(I,I) + 1.0D0 80 CONTINUE IF (NF.EQ.MF) GO TO 130 C ITERATE TO DECREASE THE NOISE LEVEL IN X. THIS SHOULD REDUCE C POSSIBL ERRORS IN NYSTROM INTERPOLATION. DERROR = ((1.0D0-R1)/R1)*EPS/NORMK IF (IFLAG.EQ.1) DERROR = DERROR*XNORM ITLOOP = 0 DO 90 I=1,MF XM(I) = XV(I) 90 CONTINUE 100 DO 110 I=1,MF XMZ(I) = XM(I) 110 CONTINUE CALL ITERT(KERNEL, RHFCN, NF, TN, WN, MF, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, NHALF, * 1) NUMR1 = NORM(XM,XMZ,MF,1) ITLOOP = ITLOOP + 1 IF ((NUMR1.GT.DERROR) .AND. (ITLOOP.LT.5)) GO TO 100 DO 120 I=1,MF XV(I) = XM(I) 120 CONTINUE C ESTIMATE NEW ERROR BOUND FOR NYSTROM INTERPOLATES. TEMP = NORMK*(R1/(1.0D0-R1))*NUMR1 IF (IFLAG.EQ.1) TEMP = TEMP/XNORM EP = DMAX1(EP,TEMP) GO TO 140 C NO ITERATION USED IN COMPUTING X. JUST COMPUTE ERROR ESTIMATE IN C NYSTROM INTERPOLATE. 130 TEMP = NORMK*ELINSY IF (IFLAG.EQ.0) TEMP = TEMP*XNORM IF (IER.NE.2) EP = DMAX1(EP,TEMP) C COMPUTE NYSTROM INTERPOLATES AT THE NODES IN T. 140 DO 160 I=1,NT SUM = 0.0D0 DO 150 J=1,MF SUM = SUM + WV(J)*KERNEL(T(I),TV(J))*XV(J) 150 CONTINUE X(I) = RHFCN(T(I)) + SUM 160 CONTINUE IF ((IER.EQ.0) .AND. (EP.GT.EPS)) IER = 4 IF ((IER.EQ.1) .AND. (EP.GT.ERROR)) IER = 5 IF ((IER.EQ.3) .AND. (EP.GT.EPS)) IER = 6 EP = SAVEP RETURN END SUBROUTINE INTERP(TM, WM, XM, M, TN, WN, XN, N, KERNEL, INT 10 * RHFCN, RHS, KMN, NHALF, NUP) C USE THE VALUES OF XN(I), I=1,...,N, TO CALCULATE THE NYSTROM C INTERPOLATES XM(I), I=1,...,M. DOUBLE PRECISION KERNEL, RHFCN, TM, WM, XM, TN, WN, XN, RHS, * KMN DIMENSION TM(M), WM(M), XM(M), TN(N), WN(N), XN(N), RHS(M), * KMN(NUP,NHALF) IF (M.GT.NUP) GO TO 60 C SINCE M .LE. NUPPER, SAVE K(TM(I),TN(J))=KMN(I,J) AND C RHS(I)=RHFCN(TM(I)) FOR LATER USE IN ITERT. DO 20 I=1,M DO 10 J=1,N KMN(I,J) = WN(J)*KERNEL(TM(I),TN(J)) 10 CONTINUE 20 CONTINUE DO 30 I=1,M RHS(I) = RHFCN(TM(I)) XM(I) = RHS(I) 30 CONTINUE C CALCULATE NYSTROM INTERPOLATING FORMULA. DO 50 I=1,M DO 40 J=1,N XM(I) = XM(I) + KMN(I,J)*XN(J) 40 CONTINUE 50 CONTINUE RETURN C M .GT. NUPPER, SO SAVE JUST RHS(I) FOR LATER USE IN ITERT. C CALCULATE NYSTROM INTERPOLATING FORMULA. 60 DO 80 I=1,M RHS(I) = RHFCN(TM(I)) XM(I) = RHS(I) DO 70 J=1,N XM(I) = XM(I) + WN(J)*KERNEL(TM(I),TN(J))*XN(J) 70 CONTINUE 80 CONTINUE RETURN END SUBROUTINE ITERT(KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, ITE 10 * XMZ, KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, NUP, * NHALF, IFLG) C THIS ROUTINE CALCULATES ONE ITERATE XM GIVEN THE INITIAL GUESS C XMZ. THE ROUTINE IS DIVIDED ACCORDING TO WHETHER OR NOT C M .GT. NUPPER. DOUBLE PRECISION KERNEL, RHFCN, TN, WN, TM, WM, XM, XMZ, KMM, * KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, SUM, DET DIMENSION TN(N), WN(N), TM(M), WM(M), XM(M), XMZ(M), * KMM(NUP,NUP), KMN(NUP,NHALF), KNM(NHALF,NUP), RHS(M), * IMKNN(NUP,NUP), LUFACT(NUP,NUP), R(M), RH(N), DELN(N) C M .GT. NUPPER MEANS THAT THE MATRICES KMM,KMN,KNM CAN NO LONGER C BE STORED DUE TO LACK OF SPACE. IF (M.GT.NUP) GO TO 120 IF (IFLG.EQ.1) GO TO 40 C IF IFLG=0, THEN THE MATRICES KMM AND KNM MUST BE COMPUTED C AND STORED. DO 30 J=1,M DO 10 I=1,M KMM(I,J) = WM(J)*KERNEL(TM(I),TM(J)) 10 CONTINUE DO 20 I=1,N KNM(I,J) = WM(J)*KERNEL(TN(I),TM(J)) 20 CONTINUE 30 CONTINUE C COMPUTE RESIDUALS R(I)=RHFCN(TM(I))-XMZ(I)+KM(TM(I))*XMZ(I) 40 DO 60 I=1,M SUM = 0.0D0 DO 50 J=1,M SUM = SUM + KMM(I,J)*XMZ(J) 50 CONTINUE R(I) = RHS(I) - (XMZ(I)-SUM) 60 CONTINUE C COMPUTE RH=KM*R AT ALL TN(I). DO 80 I=1,N RH(I) = 0.0D0 DO 70 J=1,M RH(I) = RH(I) + KNM(I,J)*R(J) 70 CONTINUE 80 CONTINUE C CALCULATE DELN=((I-KN)**(-1))*KM*R AT ALL TN(I). C ******************************************************************* C * CALL LINSYS(IMKNN, LUFACT, N, RH, DELN, 4, DET, NUP) C * C SEE THE ORIGINAL REFERENCE IN IEGS. * C ******************************************************************* C CALCULATE NEW XM. DO 110 I=1,M SUM = 0.0D0 DO 90 J=1,M SUM = SUM + KMM(I,J)*R(J) 90 CONTINUE DO 100 J=1,N SUM = SUM + KMN(I,J)*DELN(J) 100 CONTINUE XM(I) = SUM + R(I) + XMZ(I) 110 CONTINUE RETURN C ENTRANCE WHEN M .GT. NUP. C CALCULATE RESIDUALS. 120 DO 140 I=1,M SUM = 0.0D0 DO 130 J=1,M SUM = SUM + WM(J)*KERNEL(TM(I),TM(J))*XMZ(J) 130 CONTINUE R(I) = RHS(I) - (XMZ(I)-SUM) 140 CONTINUE C CALCULATE RH=KM*R. DO 160 I=1,N RH(I) = 0.0D0 DO 150 J=1,M RH(I) = RH(I) + WM(J)*KERNEL(TN(I),TM(J))*R(J) 150 CONTINUE 160 CONTINUE C ******************************************************************* C * CALL LINSYS(IMKNN, LUFACT, N, RH, DELN, 4, DET, NUP) C * C SEE THE ORIGINAL REFERENCE IN IEGS. * C ******************************************************************* C CALCULATE XM. DO 190 I=1,M SUM = 0.0D0 DO 170 J=1,M SUM = SUM + WM(J)*KERNEL(TM(I),TM(J))*R(J) 170 CONTINUE DO 180 J=1,N SUM = SUM + WN(J)*KERNEL(TM(I),TN(J))*DELN(J) 180 CONTINUE XM(I) = SUM + R(I) + XMZ(I) 190 CONTINUE RETURN END DOUBLE PRECISION FUNCTION NORM(X, Y, N, IFLAG) NOR 10 C IFLAG=0 CALCULATE THE MAXIMUM NORM OF X. C IFLAG=1 CALCULATE THE MAXIMUM NORM OF X-Y. DOUBLE PRECISION X, Y, DMAX1, DABS DIMENSION X(N), Y(N) IF (IFLAG.EQ.1) GO TO 20 C FIND THE NORM OF X. NORM = 0.0D0 DO 10 I=1,N NORM = DMAX1(NORM,DABS(X(I))) 10 CONTINUE RETURN C FIND THE NORM OF X-Y. 20 NORM = 0.0D0 DO 30 I=1,N NORM = DMAX1(NORM,DABS(X(I)-Y(I))) 30 CONTINUE RETURN END SUBROUTINE LINSYS(A, D, N, B, X, OPTION, DET, MACHIN) LIN 10 C SOLVE A*X=B, ORDER(A)=N, DIMENSION OF A=MACHIN. C OPTION=1 SOLVE A*X=B, LEAVE THE LU DECOMPOSITION A=L*U IN D C AND THE PIVOTS IN PIVOT. THE ANSWERS ARE LEFT IN X. C IT IS PERMISSABLE TO LET B=X AND D=A, BUT THEN THE C ORIGINAL CONTENTS OF A AND B ARE LOST. C OPTION=2 CALCULATE DECOMPOSITION A=L*U INCLUDING PIVOTS. SOLVE C A*X=B, AND THEN CALCULATE THE RESIDUAL AND ONE C CORRECTION. THE CORRECTIONS ARE LEFT IN R, THE NEW C VALUE OF X1 IN X, THE RELATIVE ERROR C NORM(X0-X1)/NORM(X1) C IN THE VARIABLE ERROR, AND THE RELATIVE RESIDUAL C NORM(RESIDUAL)/NORM(B) C IN THE VARIABLE RELRSD. THESE VALUES CAN BE OBTAINED C USING THE COMMON/LINEAR/ GIVEN BELOW. C OPTION=3 SAME AS OPTION=1, EXCEPT A=L*U IS KNOWN AND STORED C IN D. C OPTION=4 SAME AS OPTION=2, EXCEPT A=L*U IS KNOWN AND STORED C IN D. C THE DECOMPOSITION OF A INTO L*U USES SCALED PARTIAL PIVOTING IN C THE COLUMNS. FOR OPTIONS 1 AND 2, THE DETERMINANT OF A IS C CALCULATED AND STORED IN DET. IF DET=0, THEN THE ANSWERS ARE C NONSENSE, D AND X DO NOT CONTAIN USEFUL INFORMATION, AND AND A C AND B ARE LEFT UNDISTURBED (UNLESS D=A OR X=B). C LINSYS IS PRESENTLY LIMITED TO N .LE. 100. TO REMOVE THIS C RESTRICTION,CHANGE THE DIMENSION STATEMENT FOR THE ARRAYS R, C SCALE, AND PIVOT, WHICH ARE GIVEN IN COMMON/LINEAR/. INTEGER OPTION, PIVOT DOUBLE PRECISION NORMX, NORME, NORMB, NORMR, A, D, B, X, R, * SCALE, ERROR, RELRSD, DET, C, TEMP, DMAX1, DABS, SUM DIMENSION A(MACHIN,MACHIN), D(MACHIN,MACHIN), B(N), X(N) COMMON /LINEAR/ R(100), SCALE(100), ERROR, RELRSD, PIVOT(100) ISWIT = 1 IF (OPTION.GT.2) GO TO 110 C PRODUCE LU DECOMPOSITION OF A AND DET(A) DET = 0.0D0 DO 20 I=1,N DO 10 J=1,N D(I,J) = A(I,J) 10 CONTINUE 20 CONTINUE C PRODUCE SCALING FACTORS FOR SCALED PARTIAL PIVOTING. DO 40 I=1,N SCALE(I) = 0.0D0 DO 30 J=1,N SCALE(I) = DMAX1(SCALE(I),DABS(D(I,J))) 30 CONTINUE IF (SCALE(I).EQ.0.0D0) RETURN 40 CONTINUE DET = 1.0D0 NM1 = N - 1 C BEGIN GAUSSIAN ELIMINATION. DO 100 K=1,NM1 C SELECT PIVOT ROW. C = DABS(D(K,K))/SCALE(K) INDEX = K KP1 = K + 1 DO 50 I=KP1,N TEMP = DABS(D(I,K))/SCALE(I) IF (TEMP.LE.C) GO TO 50 INDEX = I C = TEMP 50 CONTINUE PIVOT(K) = INDEX IF (INDEX.EQ.K) GO TO 70 C SWITCH ROWS OF D. DO 60 J=K,N TEMP = D(K,J) D(K,J) = D(INDEX,J) D(INDEX,J) = TEMP 60 CONTINUE TEMP = SCALE(K) SCALE(K) = SCALE(INDEX) SCALE(INDEX) = TEMP DET = -DET 70 DET = DET*D(K,K) C ELIMINATE UNKNOWN =K FROM BELOW DIAGONAL IN COLUMN K. IF (DET.EQ.0.0D0) RETURN DO 90 I=KP1,N D(I,K) = D(I,K)/D(K,K) TEMP = D(I,K) DO 80 J=KP1,N D(I,J) = D(I,J) - TEMP*D(K,J) 80 CONTINUE 90 CONTINUE 100 CONTINUE DET = DET*D(N,N) IF (DET.EQ.0.0D0) RETURN C THE DECOMPOSITION A=P*L*U IS COMPLETE. C BEGIN SOLUTION OF LINEAR SYSTEM P*L*U*X=B C SET R AND X FOR SOLUTION OF A*X=B. 110 DO 120 I=1,N R(I) = B(I) X(I) = 0.0D0 120 CONTINUE GO TO 170 C SET R=B-A*X FOR COMPUTATION OF CORRECTION. 130 DO 150 I=1,N SUM = 0.0D0 DO 140 J=1,N SUM = SUM + A(I,J)*X(J) 140 CONTINUE R(I) = B(I) - SUM 150 CONTINUE NORMB = 0.0D0 NORMR = 0.0D0 DO 160 I=1,N NORMB = DMAX1(NORMB,DABS(B(I))) NORMR = DMAX1(NORMR,DABS(R(I))) 160 CONTINUE RELRSD = NORMR/NORMB ISWIT = 2 C SOLVE P*L*Z=R AND STORE IN R. 170 NM1 = N - 1 DO 200 K=1,NM1 IF (PIVOT(K).EQ.K) GO TO 180 KPIV = PIVOT(K) TEMP = R(K) R(K) = R(KPIV) R(KPIV) = TEMP 180 KP1 = K + 1 DO 190 I=KP1,N R(I) = R(I) - D(I,K)*R(K) 190 CONTINUE 200 CONTINUE C SOLVE U*E=R AND STORE IN R. ALSO LET X=X+E. R(N) = R(N)/D(N,N) GO TO (210, 220, 210, 220), OPTION 210 X(N) = R(N) GO TO 230 220 X(N) = X(N) + R(N) 230 DO 270 II=1,NM1 I = N - II SUM = 0.0D0 IP1 = I + 1 DO 240 J=IP1,N SUM = SUM + D(I,J)*R(J) 240 CONTINUE R(I) = (R(I)-SUM)/D(I,I) GO TO (250, 260, 250, 260), OPTION 250 X(I) = R(I) GO TO 270 260 X(I) = X(I) + R(I) 270 CONTINUE C SOLUTION OF LINEAR SYSTEM IS COMPLETE. RETURN FOR OPTION=1,3. GO TO (280, 290, 280, 290), OPTION 280 RETURN 290 GO TO (130, 300), ISWIT C CALCULATE ERRORS BASED ON CORRECTION. 300 NORMX = 0.0D0 NORME = 0.0D0 DO 310 I=1,N NORMX = DMAX1(NORMX,DABS(X(I))) NORME = DMAX1(NORME,DABS(R(I))) 310 CONTINUE ERROR = NORME/NORMX RETURN END