C TITLE: TEST OF IESIMP: C --------------------- C C THIS IS A TEST PROGRAM FOR THE INTEGRAL EQUATION PROGRAM C 'DIESMP'. IT SOLVES THE INTEGRAL EQUATION C C B C X(S) - LAMBDA*INTEGRAL FCNK(S,T)*X(T)*DT = FCNY(S) C A C C THERE ARE SEVERAL KERNEL FUNCTIONS 'FCNK' AND RIGHT-HAND FUNCTIONS C 'FCNY' WHICH CAN BE CALLED UPON BY SETTING 'NUMKER' AND 'NUMFNY', C RESPECTIVELY. THE FUNCTIONS FCNK AND FCNY USUALLY INVOLVE A C PARAMETER, CALLED 'CONST' IN THE BELOW INPUT STATEMENT. THERE IS C ALSO A PARAMETER 'LAMBDA', WHICH IS MULTIPLIED TIMES THE INTEGRAL C OPERATOR IN THE EQUATION. C 'EPS' IS THE DESIRED ERROR, AND 'IFLAG' SIGNIFIES WHETHER C EPS IS INTERPRETED AS AN ABSOLUTE ERROR(IFLAG=0) OR A RELATIVE C ERROR(IFLAG=1). C IN DIESMP, SIMPSON'S RULE IS USED TO APPROXIMATE THE C INTEGRAL OPERATOR, AND THEN THE RESULTING EQUATION IS REDUCED C TO AN EQUIVALENT LINEAR SYSTEM. THE PARAMETER 'N' IS THE INITIAL C NUMBER OF NODE POINTS TO USE IN THE SIMPSON RULE. AS A DEFAULT, C SET N=3. N MUST BE AN ODD INTEGER. C 'R2FLAG' INDICATES THE POSSIBLE PRESENCE OF PERIODICITY C IN THE EQUATION. IF R2FLAG=0, THEN NO PERIODICITY IS ASSUMED. C IF HOWEVER THE FUNCTION FCNK(S,T)*X(T) IS PERIODIC IN T ON C [A,B], FOR ALL S, THEN SET R2FLAG=1. THIS SHOULD RESULT IN C IMPROVED ERROR ESTIMATES. C SEE THE INTRODUCTORY COMMENTS OF IESIMP FOR MORE C INFORMATION. C THE PROGRAM 'DIESMP' USES THE LINPACK PROGRAMS 'DGECO' C 'DGESL'. IF THESE ARE REPLACED BY OTHER PROGRAMS, THEN CARE C SHOULD BE TAKEN THAT THE REPLACEMENT PROGRAMS FUNCTION IN C THE SAME WAY AS THE LINPACK PROGRAMS. C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER Q DOUBLE PRECISION LAMBDA,INFO INTEGER R2FLAG COMMON/PARAM/LAMBDA,CONST,NUMKER,NUMFNY PARAMETER (MAXN=33, MAXM=257, MAXWRK=MAXN*(4*MAXN+11)+9*MAXM) PARAMETER (ZERO=0.0D0, ONE=1.0D0) DIMENSION X(MAXM),WORK(MAXWRK),IWORK(MAXN),INFO(3) EXTERNAL FCNK,FCNY C C******************************************************************* C * C FOLLOWING IS THE OUTPUT UNIT NUMBER TO BE USED IN THE * C OUTPUT OF THE RESULTS TO THE FILE 'IESIMP.OUT'. * C IWRITE=8 OPEN (IWRITE,FILE='IESIMP.OUT') C * C FOLLOWING IS THE INPUT UNIT NUMBER FOR THE TERMINAL. * C * ITERML=5 C * C******************************************************************* C C READ INPUT PARAMETERS. 10 PRINT *, 'GIVE THE INDICES NUMKER AND NUMFNY.' READ *, NUMKER,NUMFNY PRINT *, 'GIVE THE PROBLEM PARAMETERS LAMBDA, CONST, AND EPS.' READ *, LAMBDA,CONST,EPS PRINT *, 'GIVE THE INTERVAL [A,B].' READ *, A,B PRINT *, 'GIVE THE PARAMETERS N, IFLAG, AND R2FLAG.' READ *, N,IFLAG,R2FLAG C C PRINT THE INPUT PARAMETERS OF THE PROBLEM. WRITE(IWRITE,101) NUMKER,NUMFNY,LAMBDA,CONST,A,B, * EPS,N,IFLAG,R2FLAG C C CALL DIESMP TO SOLVE INTEGRAL EQUATION, AND PRINT ERROR RESULTS. CALL DIESMP(FCNK,FCNY,A,B,N,EPS,IFLAG,X,MAXN,MAXM,WORK, * IWORK,R2FLAG,INFO,IER) NFINAL = INFO(3) WRITE(IWRITE,104)EPS,IER,N,NFINAL,INFO(1),INFO(2) C C PRINT COMPUTED RESULTS AND PRINT THE TRUE ABSOLUTE ERRORS. WRITE(IWRITE,102) XNORM=ZERO DO 20 I=1,N 20 XNORM=MAX(XNORM,ABS(X(I))) C ERRMAX=ZERO DO 40 I=1,N ANS=TRUE(WORK(I)) ERROR=ANS-X(I) ERRMAX=MAX(ERRMAX,ABS(ERROR)) IF(IFLAG .EQ. 0) THEN C ERRORS ARE TO BE ABSOLUTE ERRORS. WRITE(IWRITE,103) WORK(I),X(I),ANS,ERROR ELSE C ERRORS ARE TO BE RELATIVE ERRORS. WRITE(IWRITE,103) WORK(I),X(I),ANS,ERROR/XNORM END IF 40 CONTINUE C C PRINT MAXIMUM ABSOLUTE ERROR. WRITE(IWRITE,105) ERRMAX WRITE(IWRITE,106) ERRMAX/XNORM C C CHECK ON WHETHER TO CONTINUE WITH ANOTHER EQUATION. PRINT *, 'DO YOU WANT TO INPUT ANOTHER SET OF PARAMETERS? (Y/N)' READ(*,'(A1)') Q IF((Q .EQ. 'Y') .OR. (Q .EQ. 'y')) THEN GO TO 10 ELSE CLOSE(IWRITE) STOP END IF C C PRINT FORMATS FOR ABOVE OUTPUT STATEMENTS. 100 FORMAT(I2,I3,3D10.0,I3,3I2) 101 FORMAT('1KERNEL NO.',I2,5X,'RHFCN NO.',I2,/,5X,'LAMBDA=', * 1PD17.10,5X,'CONSTANT=',D17.10,/,5X,'A,B = ',D17.10, * D20.10,/,5X,'DESIRED ERROR=',D17.10,5X,'INITIAL N=',I2, * /,5X,'IFLAG=',I1,5X,'R2FLAG=',I1,//) 102 FORMAT(//,4X,1HT,19X,9HAPPROX(T),11X,7HTRUE(T),13X,5HERROR,/) 103 FORMAT(1PD20.10,2D20.10,D12.2) 104 FORMAT(17H PREDICTED ERROR=,1PD8.2,5X,4HIER=,I1,5X,8HFINAL M=, * I3,/,10X,8HFINAL N=,I2,5X,9HFINAL R1=,D8.2,5X,9HFINAL R2=, * D8.2) 105 FORMAT(//,' ACTUAL ERROR=',1PD8.2) 106 FORMAT(' ACTUAL RELATIVE ERROR=',1PD8.2,//) END FUNCTION FCNK(S,T) C ------------- C C DEFINE THE KERNEL FUNCTION OF THE INTEGRAL EQUATION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LAMBDA COMMON/PARAM/LAMBDA,C,NUMKER,NUMFNY C THE VARIABLE 'C' IS 'CONST' IN THE MAIN PROGRAM. PARAMETER (ONE=1.0D0, TWO=2.0D0) DATA PI/3.14159265358979D0/ C GO TO (10,20,30,40,50),NUMKER C C USE WITH NUMFNY=1. USE [A,B]=[0,C]. 10 FCNK=LAMBDA*COS(PI*S*T) RETURN C C USE WITH NUMFNY=5. USE [A,B]=[0,1]. 20 FCNK=C/(C*C+(S-T)*(S-T))*LAMBDA RETURN C C USE WITH NUMFNY=2 OR 3. USE [A,B]=[0,1]. 30 C2=C*C FCNK=(ONE-C2)/(ONE+C2-TWO*C*COS(TWO*PI*(S+T)))*LAMBDA RETURN C C USE WITH NUMFNY=6. USE [A,B]=[0,1]. 40 IF(S .LE. T) THEN FCNK=-S*(ONE-T)*LAMBDA ELSE FCNK=-T*(ONE-S)*LAMBDA END IF RETURN C C USE WITH NUMFNY=4. USE [A,B]=[0,1]. 50 FCNK=LAMBDA/(S+S+T+C) RETURN END FUNCTION FCNY(S) C ------------- C C DEFINE THE FORCING TERM IN THE INTEGRAL EQUATION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LAMBDA COMMON/PARAM/LAMBDA,C,NUMKER,NUMFNY C THE VARIABLE 'C' IS 'CONST' IN THE MAIN PROGRAM. PARAMETER (ONE=1.0D0, TWO=2.0D0, SIX=6.0D0, SEVEN=7.0D0, * THREE=3.0D0, P8=.8D0, P06=.06D0, P5=.5D0) DATA PI/3.14159265358979D0/ C GO TO (10,20,30,40,50,60),NUMFNY C C USE WITH NUMKER=1. 10 T1=SEVEN+PI*S T2=SEVEN-PI*S EB=DEXP(C) TEM1=(EB*(COS(T1*C)+T1*SIN(T1*C))-ONE)/(ONE+T1*T1) TEM2=(EB*(COS(T2*C)+T2*SIN(T2*C))-ONE)/(ONE+T2*T2) FCNY=EXP(S)*COS(SEVEN*S)-P5*LAMBDA*(TEM1+TEM2) RETURN C C USE WITH NUMKER=3. 20 FCNY =(ONE-LAMBDA*C)*COS(TWO*PI*S) RETURN C C USE WITH NUMKER=3. 30 FCNY =(ONE-LAMBDA*C*C*C)*COS(SIX*PI*S) RETURN C C USE WITH NUMKER=5. 40 TEMP=SQRT(C+S+S) TEMP=TWO*(ONE-TEMP*ATAN(ONE/TEMP)) FCNY=SQRT(S)-LAMBDA*TEMP RETURN C C USE WITH NUMKER=2. 50 A=-P8 B=P06 C2=C*C T1=LOG((C2+(ONE-S)**2)/(C2+S*S)) T2=ATAN((ONE-S)/C)+ATAN(S/C) FCNY =C+C*(S+A/TWO)*T1+(B-C2+S*A+S*S)*T2 FCNY =-LAMBDA*FCNY +B+S*(A+S) RETURN C C USE WITH NUMKER=4. 60 I=C T1=S**I T2=(ONE-S*T1)/(C+ONE)-(ONE-S*S*T1)/(C+THREE) T1=T1*(ONE-S)+LAMBDA*S*T2/(C+TWO) FCNY =C*C*T1 RETURN END FUNCTION TRUE(S) C ------------- C C THIS GIVES THE TRUE SOLUTION OF THE INTEGRAL EQUATION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LAMBDA COMMON/PARAM/LAMBDA,C,NUMKER,NUMFNY C THE VARIABLE 'C' IS 'CONST' IN THE MAIN PROGRAM. PARAMETER (ONE=1.0D0, TWO=2.0D0, SIX=6.0D0, SEVEN=7.0D0, * P8=.8D0, P06=.06D0) DATA PI/3.14159265358979D0/ C GO TO (10,20,30,40,50,60),NUMFNY C C USE WITH NUMKER=1. 10 TRUE=EXP(S)*COS(SEVEN*S) RETURN C C USE WITH NUMKER=3. 20 TRUE=COS(TWO*PI*S) RETURN C C USE WITH NUMKER=3. 30 TRUE=COS(SIX*PI*S) RETURN C C USE WITH NUMKER=5. 40 TRUE=SQRT(S) RETURN C C USE WITH NUMKER=2. 50 TRUE=P06+S*(-P8+S) RETURN C C USE WITH NUMKER=4. 60 IEXP=C TRUE=(S**IEXP)*(ONE-S)*C*C RETURN END SUBROUTINE DIESMP (KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, NUPPER, * MUPPER, W, IW, IFLGR2, INFO, IER) C***BEGIN PROLOGUE DIESIMP C***PURPOSE THIS ROUTINE IS FOR SOLVING LINEAR FREDHOLM INTEGRAL C EQUATIONS OF THE SECOND KIND. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DIESMP-D) C***KEYWORDS NUMERICAL ANALYSIS, LINEAR INTEGRAL EQUATIONS, C AUTOMATIC ALGORITHM, NYSTROM METHOD C***AUTHOR ATKINSON, KENDALL, (UNIVERSITY OF IOWA) C ADDRESS: DEPARTMENT OF MATHEMATICS, C UNIVERSITY OF IOWA, IOWA CITY, IOWA 52242. C E_MAIL: BLAKEAPD@UIAMVS.BITNET C (ASSISTED BY DAVID CHIEN, GRADUATE STUDENT, U. OF IOWA) C***DESCRIPTION C C *USAGE: C C INTEGER NZ, IFLAG, NUPPER, MUPPER, IW(NUPPER), IFLGR2, IER C DOUBLE PRECISION KERNEL, RHFCN, A, B, EPS, X(MUPPER), W(*), C INFO(3) C EXTERNAL KERNEL, RHFCN C C CALL DIESIMP (KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, NUPPER, C MUPPER, W, IW, IFLGR2, INFO, IER) C C *ARGUMENTS: C C KERNEL :EXT THIS IS A DOUBLE PRECISION FUNCTION OF TWO VARIABLES. C IT MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C PROGRAM CALLING DIESMP. C C RHFCN :EXT THIS IS A DOUBLE PRECISION FUNCTION OF ONE VARIABLE. C IT MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C PROGRAM CALLING DIESMP. C C A :IN IS THE LOWER INTEGRATION LIMIT. C C B :IN IS THE UPPER INTEGRATION LIMIT. C C NZ :INOUT THE INITIAL VALUE OF N IN THE PROGRAM. N IS THE ORDER OF C AN INVERSE MATRIX WHICH IS BEING USED TO ITERATIVELY C SOLVE A LARGER SYSTEM OF ORDER M, WHICH APPROXIMATES THE C INTEGRAL EQUATION. THE CHOICE OF NZ MUST ALWAYS BE ODD C AND GREATER THAN TWO. FOR A DEFAULT CHOICE, SET NZ=3. ON C COMPLETION OF THE PROGRAM, NZ IS SET EQUAL TO THE FINAL C NUMBER OF NODE POINTS USED IN OBTAINING THE SOLUTION. C THE ARRAYS W AND X WILL CONTAIN THE NZ FINAL NODE POINTS C AND CORRESPONDING SOLUTION VALUES. C C EPS :INOUT THE DESIRED ERROR, INTERPRETED ACCORDING TO IFLAG. THIS C VARIABLE IS CHANGED ON COMPLETION OF THE PROGRAM. SEE C THE DISCUSSION OF IFLAG AND IER FOR MORE INFORMATION. C C IFLAG :IN IFLAG = 0: C EPS IS TO BE INTERPRETED AS AN ABSOLUTE ERROR. C IFLAG = 1: C EPS IS TO BE INTERPRETED AS A RELATIVE ERROR. C C X :OUT 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 C NUPPER :IN AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM. N IS C THE ORDER OF APPROXIMATE INVERSE WHICH IS BEING USED TO C ITERATIVELY SOLVE A LARGER SYSTEM OF ORDER M. C C MUPPER :IN AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM. M IS C THE NUMBER OF NODE POINTS BEING USED IN THE SIMPSON RULE C APPROXIMATION OF THE INTEGRAL OPERATOR. N AND M ARE C ALWAYS ODD INTEGERS, GREATER THAN TWO. C C W :OUT THIS IS TEMPORARY WORKING STORAGE. IT SHOULD HAVE C DIMENSION 4*NU*NU+11*NU+9*MU, WITH NU=NUPPER, C MU=MUPPER. ON EXIT, W WILL CONTAIN THE FINAL CHOICE OF C NODE POINTS,CORRESPONDING TO THE APPROXIMATE SOLUTION C VALUES STORED IN X. C C IW :WORK THIS IS TEMPORARY WORKING STORAGE FOR INTEGER VARIABLES. C IT SHOULD HAVE DIMENSION = NUPPER. C C IFLGR2 :IN IFLGR2 = 0: C THE NORMAL SETTING. C IFLGR2 = 1: C IF THE INTEGRAND K(S,T)*X(T), REGARDED AS A FUNCTION OF C T, IS KNOWN TO BE PERIODIC ALONG WITH A NUMBER OF ITS C DERIVATIVES ON THE INTERVAL (A,B), FOR ALL S IN THE C INTERVAL, THEN THE ERROR ESTIMATES WILL BE MORE ACCURATE C IF IFLGR2=1. THIS WILL GENERALLY RESULT IN MORE RAPID C CONVERGENCE OF THE PROGRAM. C C INFO :OUT GIVES ADDITIONAL INFORMATION ABOUT THE FUNCTIONING OF C DIESMP. INFO(1) GIVES THE FINAL ITERATIVE RATE OF C CONVERGENCE FOR SOLVING THE LINEAR SYSTEMS OF ORDER M. C INFO(2) GIVES THE RATE OF CONVERGENCE OF THE NYSTROM C METHOD. FOR A SMOOTH KERNEL FUNCTION AND SMOOTH SOLUTION C FUNCTION X(S), THE VALUE OF INFO(2) SHOULD BE 0.0625 FOR C SIMPSONS RULE. INFO(3) GIVES THE FINAL VALUE OF N USED C IN ITERATIVELY SOLVING THE LARGER SYSTEMS OF ORDER M. IF C THE VALUE OF INFO(3) EQUALS THE FINAL VALUE OF NZ, THEN C ITERATION WAS NOT INVOKED SUCCESSFULLY. C C IER :OUT IS AN ERROR COMPLETION CODE WITH THE FOLLOWING POSSIBLE C VALUES: C IER = 0: C THE ERROR TEST WAS SATISFIED. THE PREDICTED ERROR IS C STORED IN EPS. C IER = 1: C THE ERROR TEST WAS NOT SATISFIED. THE PREDICTED ERROR IS C IN EPS. C IER = 2: C THE ERROR TEST WAS NOT SATISFIED. THE VARIABLE EPS HAS C BEEN SET TO ZERO. C REGARDLESS OF THE VALUE OF IER, THE MOST ACCURATE C SOLUTION OBTAINED IS STORED IN X. C IER = 3: C IF THE LINEAR SYSTEM BEING SOLVED IS EXACTLY SINGULAR, C THEN THE ROUTINE 'DIESMP' IS ABORTED. NO MEANINGFUL C RESULTS CAN BE EXPECTED FROM ANY OF THE OUTPUT VARIABLES C . THIS CASE SHOULD NOT OCCUR. C C *DESCRIPTION: C C THE INTEGRAL EQUATION BEING SOLVED IS C C B C X(S) - INT KERNEL(S,T)*X(T)*DT = RHFCN(S) C A C C THE METHOD BEING USED IS BASED ON THE NYSTROM METHOD WITH SIMPSON'S C RULE, WITH AN ITERATIVE TECHNIQUE OF SOLUTION FOR THE RESULTING C LINEAR SYSTEM. C C***REFERENCES 1."AN AUTOMATIC PROGRAM FOR FREDHOLM INTEGRAL EQUATIONS C OF THE SECOND KIND", ACM TRANSACTIONS ON MATH. C SOFTWARE 2, (1976), PP.154-171. C 2."A SURVEY OF NUMERICAL METHODS FOR THE SOLUTION OF C FREDHOLM INTEGRAL EQUATIONS OF THE SECOND KIND", SIAM C PUB., 1976, PART II, CHAP. 5. C***ROUTINES CALLED DGECO, DGESL, DIESP, DINTERP, DITERT, DNORM, DWANDT C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880901 UPGRADE TO FORTRAN 77. C 901001 MODIFIED TO MEET SLATEC PROLOGUE STANDARS AND TO CORRECT C SMALL ERROR. C***END PROLOGUE DIESMP C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL,INFO DIMENSION X(*),W(*),IW(*),INFO(3) EXTERNAL KERNEL,RHFCN C C DATA CUTOFF/0.5D0/ C C********************************************************************* C * DATA RATIO/0.0625D0/, ROOTRT/0.25D0/ C * C THESE CONSTANTS DEPEND ON THE NUMERICAL INTEGRATION RULE * C BEING USED. SEE THE REFERENCES OF DESMP 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 C SET UP RELATIVE BASE ADDRESSES FOR THE VARIOUS ARRAYS INTO WHICH C W IS TO BE SPLIT. C C***FIRST EXECUTABLE STATEMENT DIESMP C 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 I21=I20+NSQ NHALF=(N+1)/2 C CALL DIESP(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), * W(I21),INFO,IW) RETURN END SUBROUTINE DIESP (KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, NUP, * MUP, IFLGR2, IER, RATIO, ROOTRT, CUTOFF, NHALF, T, LUFACT, * KMM, KMN, KNM, RHS, R, RH, DELN, TM, TN, XM, XMZ, WM, WN, * OLDX, SAVETM, XN, SAVEWM, ASIDE, ASIDE3, Z, INFO, IPIVOT) C***BEGIN PROLOGUE DIESP C***SUBSIDIARY C***PURPOSE THIS ROUTINE CONTROLS THE SOLUTION OF THE INTEGRAL EQUATION C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DIESP-D) C***AUTHOR ATKINSON, KENDALL, (UNIVERSITY OF IOWA) C ADDRESS: DEPARTMENT OF MATHEMATICS, C UNIVERSITY OF IOWA, IOWA CITY, IOWA 52242. C E_MAIL: BLAKEAPD@UIAMVS.BITNET C (ASSISTED BY DAVID CHIEN, GRADUATE STUDENT, U. OF IOWA) C***DESCRIPTION C C *USAGE: C C INTEGER NZ, IFLAG, NUP, MUP, IFLGR2, IER, NHALF, IPIVOT(NUP) C DOUBLE PRECISION KERNEL, RHFCN, A, B, EPS, X(MUP), RATIO, ROOTRT, C CUTOFF, T(MUP), LUFACT(NUP,NUP), KMM(NUP,NUP), C KMN(NUP,NHALF), KNM(NHALF,NUP), RHS(MUP), R(MUP), C RH(NUP), DELN(NUP), TM(MUP), TN(NUP), XM(MUP), C XMZ(MUP), WM(MUP), WN(NUP), OLDX(MUP), C SAVETM(MUP), XN(NUP), SAVEWM(MUP), ASIDE(NUP,4), C ASIDE3(NUP,NUP), Z(NUP), INFO(3) C EXTERNAL KERNEL, RHFCN C PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, P5 = 0.5D0) C C CALL DIESP (KERNEL, RHFCN, A, B, NZ, EPS, IFLAG, X, NUP, MUP, C IFLGR2, IER, RATIO, ROOTRT, CUTOFF, NHALF, T, LUFACT, C KMM, KMN, KNM, RHS, R, RH, DELN, TM, TN, XM, XMZ, WM, C WN, OLDX,SAVETM, XN, SAVEWM, ASIDE, ASIDE3, Z, INFO, C IPIVOT) C C *ARGUMENT: C C KERNEL :EXT THIS IS A DOUBLE PRECISION FUNCTIONS OF TWO VARIABLES. C IT MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C PROGRAM CALLING DIESMP. C C RHFCN :EXT THIS IS A DOUBLE PRECISION FUNCTIONS OF ONE VARIABLE. C IT MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C PROGRAM CALLING DIESMP. C C A :IN IS THE LOWER INTEGRATION LIMIT. C C B :IN IS THE UPPER INTEGRATION LIMIT. C C NZ :INOUT THE INITIAL VALUE OF N IN THE PROGRAM. N IS THE ORDER C OF AN INVERSE MATRIX WHICH IS BEING USED TO ITERATIVELY C SOLVE A LARGER SYSTEM OF ORDER M, WHICH APPROXIMATES C THE INTEGRAL EQUATION. THE CHOICE OF NZ MUST ALWAYS BE C ODD AND GREATER THAN TWO. FOR A DEFAULT CHOICE, SET C NZ=3. ONCOMPLETION OF THE PROGRAM, NZ IS SET EQUAL TO C THE FINAL NUMBER OF NODE POINTS USED IN OBTAINING THE C SOLUTION. THE ARRAYS W AND X WILL CONTAIN THE NZ FINAL C NODE POINTS AND CORRESPONDING SOLUTION VALUES. C C EPS :INOUT THE DESIRED ERROR, INTERPRETED ACCORDING TO IFLAG. THIS C VARIABLE IS CHANGED ON COMPLETION OF THE PROGRAM. SEE C THE DISCUSSION OF IFLAG AND IER FOR MORE INFORMATION. C C IFLAG :IN IFLAG = 0: C EPS IS TO BE INTERPRETED AS AN ABSOLUTE ERROR. C IFLAG = 1: C EPS IS TO BE INTERPRETED AS A RELATIVE ERROR. C C X :OUT THE ANSWER TO THE INTEGRAL EQUATION, EVALUATED AT THE C FINAL SET OF NODE POINTS, WHICH ARE STORED IN T. THE C DIMENSION OF X SHOULD BE AT LEAST MUPPER. C C NUP :IN AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM. N IS C THE ORDER OF APPROXIMATE INVERSE WHICH IS BEING USED TO C ITERATIVELY SOLVE A LARGER SYSTEM OF ORDER M.SOLUTION. C C MUP :IN AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM. M IS C THE NUMBER OF NODE POINTS BEING USED IN THE SIMPSON C RULE APPROXIMATION OF THE INTEGRAL OPERATOR. N AND M C ARE ALWAYS ODD INTEGERS, GREATER THAN TWO. C C IFLGR2 :IN IFLGR2 = 0: C THE NORMAL SETTING. C IFLGR2 = 1: C IF THE INTEGRAND K(S,T)*X(T), REGARDED AS A FUNCTION OF C T, IS KNOWN TO BE PERIODIC ALONG WITH A NUMBER OF ITS C DERIVATIVES ON THE INTERVAL (A,B), FOR ALL S IN THE C 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 C C IER :OUT IS AN ERROR COMPLETION CODE WITHTHE FOLLOWING POSSIBLE C VALUES: C IER = 0: C THE ERROR TEST WAS SATISFIED. THE PREDICTED ERROR IS C STORED IN EPS. C IER = 1: C THE ERROR TEST WAS NOT SATISFIED. THE PREDICTED ERROR C IS IN EPS. C IER = 2: C THE ERROR TEST WAS NOT SATISFIED. THE VARIABLE EPS HAS C BEEN SET TO ZERO. C REGARDLESS OF THE VALUE OF IER, THE MOST ACCURATE C SOLUTION OBTAINED IS STORED IN X. C IER = 3: C IF THE LINEAR SYSTEM BEING SOLVED IS EXACTLY SINGULAR, C THEN THE ROUTINE 'DIESMP' IS ABORTED. NO MEANINGFUL C RESULTS CAN BE EXPECTED FROM ANY OF THE OUTPUT C VARIABLES. THIS CASE SHOULD NOT OCCUR. C C RATIO :IN DEPENDS ON THE NUMERICAL INTEGRATION RULE BEING USED. C SEE THE DISCUSSION OF DIESMP IN ORDER TO SET THEM C PROPERLY WHEN CHANGING NUMERICAL INTEGRATION RULES. C RATIO DEPENDS ON THE RATE OF CONVERGENCE OF THE RULE C BEING USED. C C ROOTRT :IN DEPENDS ON THE NUMERICAL INTEGRATION RULE BEING USED. C SEE THE DISCUSSION OF DIESMP IN ORDER TO SET THEM C PROPERLY WHEN CHANGING NUMERICAL INTEGRATION RULES. C ROOTRT USUALLY EQUALS SQRT(RATIO). C C CUTOFF :IN THE DESIRED ITERATIVE RATE OF CONVERGENCE FOR SOLVING C THE LINEAR SYSTEMS OF ORDER M. WHEN INFO(1) IS SMALLER C THEN CUTOFF, THE ITERATES ARE CONVERGING VERY SLOWLY OR C NOT AT ALL, AND THE VALUE OF N CANNOT BE INCREASED ANY C FURTHER. THUS, RETURN WITHOUT ANY FURTHER ATTEMPTS TO C LESSEN THE ERROR. C C NHALF :IN NHALF = (NUP + 1)/ 2. C C T :OUT THIS IS TEMPORARY WORKING STORAGE. IT SHOULD HAVE C DIMENSION 4*NUP*NUP+11*NUP+9*MUP. ON EXIT, T WILL C CONTAIN THE FINAL CHOICE OF NODE POINTS,CORRESPONDING C TO THE APPROXIMATE SOLUTION VALUES STORED IN X. T IS C THE SAME AS W IN DIESMP. C C LUFACT :OUT THE LINEAR SYSTEM OF THE INTEGRAL EQUATION. IT IS AN C NUP BY NUP ARRAY, IS THE SAME AS W(I1) IN DIESMP. C C KMM :OUT IS THE MULTIPLICATION OF WM AND KERNEL AT (TM(I),TM(J)) C . IT IS AN NUP BY NUP ARRAY, IS THE SAME AS W(I2) IN C DIESMP. C C KMN :OUT IS THE MULTIPLICATION OF WN AND KERNEL WHEN M IS C SMALLER THAN NUP.IT IS AN NUP BY NHALF ARRAY, IS THE C SAME AS W(I3) IN DIESMP. C C KNM :OUT IS THE MULTIPLICATION OF WM AND KERNEL AT (TN(I),TM(J)) C . IT IS AN NHALF BY NUP ARRAY, IS THE SAME AS W(I4) IN C DIESMP. C C RHS :OUT ARE VALUES OF RIGHT-HAND SIDE FUNCTION AT TM(I). IT IS C AN MUP ARRAY, IS THE SAME AS W(I5) IN DIESMP. C C R :WORK IS AN MUP ARRAY, IS THE SAME AS W(I6) IN DIESMP. C C RH :WORK IS AN NUP ARRAY, IS THE SAME AS W(I7) IN DIESMP. C C DELN :WORK IS AN NUP ARRAY, IS THE SAME AS W(I8) IN DIESMP. C C TM :OUT THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH M C EVENLY SPACED NODES.IT IS AN MUP ARRAY, IS THE SAME AS C W(I9) IN DIESMP. C C TN :OUT THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH N C EVENLY SPACED NODES. IT IS AN NUP ARRAY, IS THE SAME C AS W(I10) IN DIESMP. C C XM :OUT THE NYSTROM INTERPOLATES. IT IS AN MUP ARRAY, IS THE C SAME AS W(I11) IN DIESMP. C C XMZ :OUT IS THE INITIAL GUESS FOR SOLVING XM ITERATIVELY. IT IS C AN MUP ARRAY, IS THE SAME AS W(I12) IN DIESMP. C C WM :OUT THE WEIGHTS FOR SIMPSON'S RULE ON (A, B), WITH M EVENLY C SPACED NODES. IT IS AN MUP ARRAY, IS THE SAME AS W(I13) C IN DIESMP. C C WN :OUT THE WEIGHTS FOR SIMPSON'S RULE ON (A, B), WITH N EVENLY C SPACED NODES. IT IS AN NUP ARRAY, IS THE SAME AS W(I14) C IN DIESMP. C C OLDX :INOUT IS AN MUP ARRAY, IS THE SAME AS W(I15) IN DIESMP. C C SAVETM :WORK IS AN MUP ARRAY, IS THE SAME AS W(I16) IN DIESMP. C C XN :OUT THE NYSTROM INTERPOLATES. IT IS AN NUP ARRAY, IS THE C SAME AS W(I17) IN DIESMP. C C SAVEWM :WORK IS AN MUP ARRAY, IS THE SAME AS W(I18) IN DIESMP. C C ASIDE :WORK IS AN NUP BY 4 ARRAY, IS THE SAME AS W(I19) IN DIESMP. C C ASIDE3 :WORK IS AN NUP BY NUP ARRAY, IS THE SAME AS W(I20) IN C DIESMP. C C Z :WORK IS AN NUP ARRAY, IS THE SAME AS W(I21) IN DIESMP. C C INFO :OUT GIVES ADDITIONAL INFORMATION ABOUT THE FUNCTIONING OF C DIESMP. INFO(1) GIVES THE FINAL ITERATIVE RATE OF C CONVERGENCE FOR SOLVING THE LINEAR SYSTEMS OF ORDER M. C INFO(2) GIVES THE RATE OF CONVERGENCE OF THE NYSTROM C METHOD. FOR A SMOOTH KERNEL FUNCTION AND SMOOTH C SOLUTION FUNCTION X(S), THE VALUE OF INFO(2) SHOULD BE C 0.0625 FOR SIMPSONS RULE. INFO(3) GIVES THE FINAL C VALUE OF N USED IN ITERATIVELY SOLVING THE LARGER C SYSTEMS OF ORDER M. IF THE VALUE OF INFO(3) EQUALS THE C FINAL VALUE OF NZ, THEN ITERATION WAS NOT INVOKED C SUCCESSFULLY. C C IPIVOT :WORK THIS IS TEMPORARY WORKING STORAGE FOR INTEGER VARIABLES C . IT SHOULD HAVE DIMENSION = NUP, IS THE SAME AS IW IN C DIESMP. C C***SEE ALSO DIESMP C***ROUTINES CALLED DGECO, DGESL, DINTRP, DITERT, DNORM, DWANDT C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880901 UPGRADE TO FORTRAN 77. C 901001 MODIFIED TO MEET SLATEC PROLOGUE STANDARS AND TO CORRECT C SMALL ERROR. C***END PROLOGUE DIESP C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL,LUFACT,KMM,KMN,KNM,NUMR2,DNORM,NUMR1, * INFO 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),SAVETM(MUP),XN(NUP),SAVEWM(MUP),ASIDE(NUP,4), * ASIDE2(5),ASIDE3(NUP,NUP),Z(NUP),IPIVOT(NUP),INFO(3) PARAMETER(ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0, P5=0.5D0) EXTERNAL KERNEL,RHFCN C C***FIRST EXECUTABLE STATEMENT DIESP C 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 C INITIALIZATION. C IER=0 LOOP=1 N=NZ INFO(3)=N M=TWICE(N) C C *** STAGE A *** C C DIRECT SOLUTION OF LINEAR SYSTEM (I-KN)*XN=RHS, WHILE C TRYING TO FIND A GOOD APPROXIMATE INVERSE TO IMPLEMENT C THE ITERATIVE METHOD OF SOLUTION. C C CREATE NODES TN(I) AND WEIGHTS WN(I), I=1,...,N. CALL DWANDT(WN,TN,N,A,B) C CREATE SYSTEM (I-KN)*XN=RHFCN. DO 20 I=1,N DO 10 J=1,N 10 LUFACT(I,J)=-WN(J)*KERNEL(TN(I),TN(J)) XN(I)=RHFCN(TN(I)) 20 LUFACT(I,I)=LUFACT(I,I)+ONE C C THIS IS THE ENTRANCE TO STAGE A WHEN ITERATION FAILS. WE NEED C TO INCREASE N IN ORDER TO OBTAIN A BETTER ITERATION RATE. 30 CONTINUE C C SOLVE (I-KN)*XN=RHFCN AT ALL TN(I). ALSO OBTAIN THE LU C FACTORIZATION FOR LATER USE IN THE STAGE B ITERATIVE METHOD. C C **************************************************************** C * CALL DGECO(LUFACT,NUP,N,IPIVOT,RCOND,Z) IF (RCOND .EQ. ZERO) THEN IER=3 RETURN END IF CALL DGESL(LUFACT,NUP,N,IPIVOT,XN,0) C * C THIS LINEAR EQUATION SOLVER IS FROM LINPACK. IT MAY BE * C REPLACED BY ANOTHER PROGRAM, BUT CARE MUST BE TAKEN THAT * C THE NEW PROGRAM DOES THE SAME JOB. THIS LINEAR EQUATION * C SOLVER IS ALSO CALLED IN THE SUBROUTINE DITERT. * C * C **************************************************************** C IF (LOOP .EQ. 1) THEN R2=P5 INFO(2)=R2 R1RAT=ROOTRT GO TO 50 ELSEIF (LOOP .EQ. 2) THEN NUMR2=DNORM(XN,OLDX,N,1) DENR2=ZERO R2=P5 INFO(2)=R2 ELSE C LOOP .GT. 2. C ESTIMATE THE RATE OF CONVERGENCE OF XN TO X. COMPARE WITH C THE THEORETICAL RATE OF CONVERGENCE GIVEN IN RATIO. NUMR2=DNORM(XN,OLDX,N,1) R2=MIN(P5,NUMR2/DENR2) INFO(2)=R2 IF ((R2 .LT. RATIO) .AND. (IFLGR2 .EQ. 0)) THEN R2=RATIO INFO(2)=R2 END IF IF (IFLGR2 .EQ. 0) THEN R1RAT=ROOTRT ELSE C ALLOW FOR A FASTER RATE OF CONVERGENCE IN THE C ITERATION, DUE TO THE PERIODICITY OF THE INTEGRAND. R1RAT=MIN(SQRT(R2),ROOTRT) END IF END IF C C CHECK FOR ERROR IN XN USING TEST INVOLVING R2 AND OLDX, C BASED ON THE THEORY FOR ASYMPTOTIC ERROR BOUNDS. ERROR=(R2/(ONE-R2))*NUMR2 IF (IFLAG .EQ. 1) ERROR=ERROR/DNORM(XN,XN,N,0) C C THE FOLLOWING IS AN EMPIRICALLY BASED ADJUSTMENT. IF ((IFLGR2 .EQ. 1) .AND. (R2 .LT. RATIO)) ERROR=TWO*ERROR IF (ERROR .LE. EPS) THEN C SET UP T,X,EPS,NZ FOR SUCCESSFUL RETURN. DO 40 I=1,N X(I)=XN(I) 40 T(I)=TN(I) EPS=ERROR NZ=N RETURN END IF C C XN WAS NOT SUFFICIENTLY ACCURATE. PREPARE TO FIND XM. DENR2=NUMR2 C C INITIALIZE FOR SOLVING (I-KM)*XM=RHFCN ITERATIVELY. 50 CALL DWANDT(WM,TM,M,A,B) FLAG=0 CALL DINTRP(TM,XMZ,M,TN,WN,XN,N,KERNEL,RHFCN,RHS) DO 60 I=1,M 60 OLDX(I)=XMZ(I) C C COMPUTE TWO ITERATES. CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM, * RHS,LUFACT,R,RH,DELN,IPIVOT,NUP,NHALF,FLAG) DENR1=DNORM(XM,XMZ,M,1) FLAG=1 DO 70 I=1,M 70 XMZ(I)=XM(I) CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM, * RHS,LUFACT,R,RH,DELN,IPIVOT,NUP,NHALF,FLAG) NUMR1=DNORM(XM,XMZ,M,1) R1=NUMR1/DENR1 INFO(1)=R1 C C CHECK ON THE SIZE OF M. IF (M .LE. NUP) GO TO 100 C N CANNOT BE INCREASED ANY FURTHER, AND THEREFORE, A WEAKER C TOLERANCE IS ALLOWED FOR THE ITERATION RATE. C IF (R1 .LE. CUTOFF) GO TO 200 C C THE ITERATES ARE CONVERGING VERY SLOWLY OR NOT AT ALL, AND C THE VALUE OF N CANNOT BE INCREASED ANY FURTHER. THUS, C RETURN WITHOUT ANY FURTHER ATTEMPTS TO LESSEN THE ERROR. 80 IF (LOOP .EQ. 1) THEN EPS=ZERO IER=2 ELSE EPS=ERROR IER=1 END IF DO 90 I=1,N X(I)=XN(I) 90 T(I)=TN(I) NZ=N RETURN C C CHECK ON THE SPEED OF CONVERGENCE OF ITERATIVE METHOD. IF C IT IS SUFFICIENTLY RAPID, THEN FIX N AND GO TO STAGE B. 100 IF (R1 .GT. R1RAT) THEN C THE ITERATION IS NOT CONVERGING RAPIDLY ENOUGH. C REINITIALIZE FOR SOLVING (I-KN)*XN=RHFCN AGAIN C WITH A LARGER N. N=M INFO(3)=N LOOP=LOOP+1 M=TWICE(N) DO 120 J=1,N DO 110 I=1,N 110 LUFACT(I,J)=-KMM(I,J) WN(J)=WM(J) TN(J)=TM(J) XN(J)=RHS(J) 120 LUFACT(J,J)=LUFACT(J,J)+ONE GO TO 30 END IF C C THE ITERATION SEEMS TO BE CONVERGING SUFFICIENTLY RAPIDLY. C SAVE INFORMATION IN CASE STAGE B ABORTS AT A LARGER VALUE OF M C AND CONTROL MUST RETURN TO STAGE A. DO 130 I=1,M ASIDE(I,1)=OLDX(I) ASIDE(I,2)=WM(I) ASIDE(I,3)=TM(I) 130 ASIDE(I,4)=RHS(I) ASIDE2(1)=LOOP ASIDE2(2)=M ASIDE2(3)=R2 ASIDE2(4)=DENR2 ASIDE2(5)=R1RAT DO 140 I=1,M DO 140 J=1,M 140 ASIDE3(I,J)=KMM(I,J) C C *** STAGE B *** C C *** ITERATIVE METHOD OF SOLUTION OF (I-KM)*XM=RHS. *** C C TEST TO SEE IF THE CURRENT ITERATE XM IS SUFFICIENTLY ACCURATE C COMPARED TO THE TRUE XM. 200 TEMP=DNORM(XM,OLDX,M,1) RT=MIN(RATIO,R2) TEST=((ONE-R1)/R1)*(RT/(ONE-RT))*TEMP IF (NUMR1 .LE. TEST) GO TO 280 C C THE ITERATE WAS NOT SUFFICIENTLY ACCURATE. INITIALIZE FOR C COMPUTATION OF ANOTHER ITERATE. DENR1=NUMR1 DO 210 I=1,M 210 XMZ(I)=XM(I) CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM, * RHS,LUFACT,R,RH,DELN,IPIVOT,NUP,NHALF,FLAG) NUMR1=DNORM(XM,XMZ,M,1) R1=NUMR1/DENR1 INFO(1)=R1 IF (R1 .LE. CUTOFF) GO TO 200 C C THE FOLLOWING IS FOR THE 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 DIESMP. 220 MNEW=ASIDE2(2) IF (MNEW .GT. NUP) THEN C ABORTIVE EXIT FROM STAGE B. N CANNOT BE INCREASED C FURTHER, AND R1 IS NOT SUFFICIENTLY SMALL. IF (M .EQ. TWICE(N)) GO TO 80 DO 230 I=1,OLDM 230 T(I)=SAVETM(I) NZ=OLDM EPS=ERROR IER=1 RETURN END IF C IF (M .EQ. TWICE(N)) THEN C THERE IS NO NEED TO USE THE SAVED VARIABLES. N = M INFO(3)=N LOOP = LOOP+1 M = TWICE(N) DO 250 J=1,N DO 240 I=1,N 240 LUFACT(I,J) = -KMM(I,J) WN(J) = WM(J) TN(J) = TM(J) XN(J) = XM(J) 250 LUFACT(J,J) = LUFACT(J,J)+ONE GO TO 30 END IF C C RECONSTRUCT THE NEEDED VARIABLES FROM THOSE SAVED PREVIOUSLY, C FOR AN ABORTIVE RETURN TO STAGE A. N=MNEW INFO(3)=N DO 270 J=1,N DO 260 I=1,N 260 LUFACT(I,J)=-ASIDE3(I,J) OLDX(J)=ASIDE(J,1) WN(J)=ASIDE(J,2) TN(J)=ASIDE(J,3) XN(J)=ASIDE(J,4) 270 LUFACT(J,J)=LUFACT(J,J)+ONE M=TWICE(N) LOOP=ASIDE2(1)+ONE R2=ASIDE2(3) INFO(2)=R2 DENR2=ASIDE2(4) R1RAT=ASIDE2(5) GO TO 30 C 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. 280 IF (LOOP .EQ. 1) THEN DENR2=TEMP LOOP=2 ELSE LOOP=LOOP+1 NUMR2=TEMP R2=MIN(NUMR2/DENR2,P5) INFO(2)=R2 IF ((R2 .LT. RATIO) .AND. (IFLGR2 .EQ. 0)) THEN R2=RATIO INFO(2)=R2 END IF DENR2=NUMR2 END IF C ERROR=(R2/(ONE-R2))*TEMP IF (IFLAG .EQ. 1) ERROR=ERROR/DNORM(XM,XM,M,0) IF ((IFLGR2 .EQ. 1) .AND. (R2 .LT. RATIO)) ERROR=TWO*ERROR IF (ERROR .LE. EPS) THEN C THE ANSWER XM IS SUFFICIENTLY ACCURATE. IER=0 NZ=M EPS=ERROR DO 290 I=1,M X(I)=XM(I) 290 T(I)=TM(I) RETURN END IF C C XM IS NOT SUFFICIENTLY ACCURATE RELATIVE TO THE TRUE X. C CHECK TO SEE WHETHER M CAN BE INCREASED FURTHER. MNEW=TWICE(M) IF (MNEW .GT. MUP) THEN C M CANNOT BE INCREASED ANY FURTHER, AND THE SOLUTION C CANNOT BE FOUND WITH THE NEEDED ACCURACY. IER=1 DO 300 I=1,M X(I)=XM(I) 300 T(I)=TM(I) NZ=M EPS=ERROR RETURN END IF C C THE ESTIMATED ERROR IN XM WAS NOT SUFFICIENTLY SMALL. C M IS INCREASED AND TWO MORE ITERATES ARE COMPUTED WITH A NEW M. DO 310 I=1,M 310 X(I)=XM(I) OLDM=M M=MNEW DO 320 I=1,OLDM SAVEWM(I)=WM(I) 320 SAVETM(I)=TM(I) C CALL DWANDT(WM,TM,M,A,B) FLAG=0 CALL DINTRP(TM,XMZ,M,SAVETM,SAVEWM,XM,OLDM,KERNEL,RHFCN,RHS) DO 330 I=1,M 330 OLDX(I)=XMZ(I) CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM, * RHS,LUFACT,R,RH,DELN,IPIVOT,NUP,NHALF,FLAG) DENR1=DNORM(XM,XMZ,M,1) FLAG=1 DO 340 I=1,M 340 XMZ(I)=XM(I) CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM, * RHS,LUFACT,R,RH,DELN,IPIVOT,NUP,NHALF,FLAG) NUMR1=DNORM(XM,XMZ,M,1) C R1=NUMR1/DENR1 INFO(1)=R1 IF (R1 .LE. CUTOFF) THEN GO TO 200 ELSE GO TO 220 END IF END SUBROUTINE DITERT (KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, LUFACT, R, RH, DELN, IPIVOT, NUP, NHALF, * IFLG) C***BEGIN PROLOGUE DITERT C***SUBSIDIARY C***PURPOSE THIS ROUTINE CALCULATES ONE ITERATE XM GIVEN THE INITIAL C GUESS XMZ. THE ROUTINE IS DIVIDED ACCORDING TO WHETHER OR C NOT M .GT. NUPPER. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DITERT-D) C***AUTHOR ATKINSON, KENDALL, (UNIVERSITY OF IOWA) C ADDRESS: DEPARTMENT OF MATHEMATICS, C UNIVERSITY OF IOWA, IOWA CITY, IOWA 52242. C E_MAIL: BLAKEAPD@UIAMVS.BITNET C (ASSISTED BY DAVID CHIEN, GRADUATE STUDENT, U. OF IOWA) C***DESCRIPTION C C *USAGE: C C INTEGER N, M, NUP, NHALF, IPIVOT(N), IFLG C DOUBLE PRECISION KERNEL, RHFCN, TN(N), WN(N), TM(M), WM(M), XM(M), C XMZ(M), KMM(NUP,NUP), KMN(NUP.NHALF), C KNM(NHALF,NUP), RHS(M), LUFACT(NUP,NUP), R(M), C RH(N), DELN(N) C PARAMETER (ZERO = 0.0D0) C C CALL DITERT (KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, KMM, C KMN, KNM, RHS, LUFACT, R, RH, DELN, IPIVOT, NUP, C NHALF, IFLG) C C *ARGUMENTS: C C KERNEL :IN IS THE KERNEL FUNTION OF THE INTEGRAL EQUTION. C C RHFCN :IN IS THE RIGHT-HAND SIDE FUNCTION OF THE INTEGRAL C EQUTION. C C N :IN IS THE NUMBER OF SUBINTERVALS ON (A, B). C C TN :IN THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH N C EVENLY SPACED NODES. C C WN :IN THE WEIGHTS FOR SIMPSON'S RULE ON (A, B), WITH N C EVENLY SPACED NODES. C C M :IN THE NUMBER OF SUBINTERVALS ON (A, B), USUALLY, M IS C GREATER THEN N. C C TM :IN THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH M C EVENLY SPACED NODES. C C WM :IN THE WEIGHTS FOR SIMPSON'S RULE ON (A, B), WITH M C EVENLY SPACED NODES. C C XM :OUT THE NYSTROM INTERPOLATES. C C XMZ :IN IS THE INITIAL GUESS FOR SOLVING XM ITERATIVELY. C C KMM :OUT IS THE MULTIPLICATION OF WM AND KERNEL AT C (TM(I),TM(J)). C C KMN :OUT IS THE MULTIPLICATION OF WN AND KERNEL WHEN M IS C SMALLER THAN NUP. C C KNM :OUT IS THE MULTIPLICATION OF WM AND KERNEL AT C (TN(I),TM(J)). C C RHS :IN ARE VALUES OF RIGHT-HAND SIDE FUNCTION AT TM(I). C C LUFACT :INOUT THE LINEAR SYSTEM OF THE INTEGRAL EQUATION. C C R :WORK RESIDUALS R(I) = RHS(I) - (XMZ(I) - SUM). C C RH :WORK RH = KM*R. C C DELN :WORK THIS IS TEMPORARY WORKING STORAGE. IT SHOULD HAVE C DIMENSION = N. C C IPIVOT :WORK THIS IS TEMPORARY WORKING STORAGE FOR INTEGER C VARIABLES. IT SHOULD HAVE DIMENSION = N. C C NUP :IN AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM( C DIESMP. N IS THE ORDER OF APPROXIMATE INVERSE WHICH IS C BEING USED TO ITERATIVELY SOLVE A LARGER SYSTEM OF C ORDER M. C C NHALF :IN NHALF = (NUP + 1)/ 2 C C IFLG :IN IFLG = 0: C THE MATRICES KMM AND KNM MUST BE COMPUTED AND STORED. C C***ROUTINES CALL DGESL C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880901 UPGRADE TO FORTRAN 77. C 901001 MODIFIED TO MEET SLATEC PROLOGUE STANDARS AND TO CORRECT C SMALL ERROR. C***END PROLOGUE DITERT C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL,KMM,KMN,KNM,LUFACT PARAMETER (ZERO=0.0D0) 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),IPIVOT(N) C C***FIRST EXECUTABLE STATEMENT DITERT C 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 130 C IF(IFLG .EQ. 0) THEN C THE MATRICES KMM, KMN AND KNM MUST BE COMPUTED AND STORED. DO 20 J=1,M DO 10 I=1,M 10 KMM(I,J)=WM(J)*KERNEL(TM(I),TM(J)) DO 20 I=1,N KMN(J,I)=WN(I)*KERNEL(TM(J),TN(I)) 20 KNM(I,J)=WM(J)*KERNEL(TN(I),TM(J)) END IF C C COMPUTE RESIDUALS R(I)=RHFCN(TM(I))-XMZ(I)+KM(TM(I))*XMZ(I) DO 50 I=1,M SUM=ZERO DO 40 J=1,M 40 SUM=SUM+KMM(I,J)*XMZ(J) 50 R(I)=RHS(I)-(XMZ(I)-SUM) C C COMPUTE RH=KM*R AT ALL TN(I). DO 80 I=1,N RH(I)=ZERO DO 70 J=1,M 70 RH(I)=RH(I)+ KNM(I,J)*R(J) 80 DELN(I)=RH(I) C C CALCULATE DELN=((I-KN)**(-1))*KM*R AT ALL TN(I). C C ************************************************************** C * CALL DGESL(LUFACT,NUP,N,IPIVOT,DELN,0) C * C SEE THE ORIGINAL REFERENCE IN SUBROUTINE DIESP. * C * C ************************************************************** C C CALCULATE NEW XM. DO 120 I=1,M SUM=ZERO DO 100 J=1,M 100 SUM=SUM+KMM(I,J)*R(J) DO 110 J=1,N 110 SUM=SUM+KMN(I,J)*DELN(J) 120 XM(I)=SUM+R(I)+XMZ(I) RETURN C C ENTRANCE WHEN M .GT. NUP. C CALCULATE RESIDUALS. 130 DO 150 I=1,M SUM=ZERO DO 140 J=1,M 140 SUM=SUM+WM(J)*KERNEL(TM(I),TM(J))*XMZ(J) 150 R(I)=RHS(I)-(XMZ(I)-SUM) C CALCULATE RH=KM*R. DO 180 I=1,N RH(I)=ZERO DO 170 J=1,M 170 RH(I)=RH(I)+WM(J)*KERNEL(TN(I),TM(J))*R(J) 180 DELN(I)=RH(I) C C ************************************************************** C * CALL DGESL(LUFACT,NUP,N,IPIVOT,DELN,0) C * C SEE THE ORIGINAL REFERENCE IN SUBROUTINE DIESP. * C * C ************************************************************** C C CALCULATE XM. DO 220 I=1,M SUM=ZERO DO 200 J=1,M 200 SUM=SUM+WM(J)*KERNEL(TM(I),TM(J))*R(J) DO 210 J=1,N 210 SUM=SUM+WN(J)*KERNEL(TM(I),TN(J))*DELN(J) 220 XM(I)=SUM+R(I)+XMZ(I) RETURN END SUBROUTINE DINTRP (TM, XM, M, TN, WN, XN, N, KERNEL, RHFCN, RHS) C***BEGIN PROLOGUE DINTRP C***SUBSIDIARY C***PURPOSE USE THE VALUES OF XN(I), I=1,...,N, TO CALCULATE THE C NYSTROM INTERPOLATES XM(I), I=1,...,M. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DINTRP-D) C***AUTHOR ATKINSON, KENDALL, (UNIVERSITY OF IOWA) C ADDRESS: DEPARTMENT OF MATHEMATICS, C UNIVERSITY OF IOWA, IOWA CITY, IOWA 52242. C E_MAIL: BLAKEAPD@UIAMVS.BITNET C (ASSISTED BY DAVID CHIEN, GRADUATE STUDENT, U. OF IOWA) C***DESCRIPTION C C *USAGE: C C INTEGER M, N, C DOUBLE PRECISION WM(M), XM(M), TN(N), WN(N), XN(N), KERNEL, C RHFCN, RHS(M) C C CALL DINTRP (TM, XM, M, TN, WN, XN, N, KERNEL, RHFCN, RHS) C C *ARGUMENT: C C TM :IN THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH M C EVENLY SPACED NODES. C C XM :OUT THE NYSTROM INTERPOLATES. C C M :IN THE NUMBER OF SUBINTERVALS ON (A, B), USUALLY, M IS C GREATER THEN N. C C TN :IN THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH N C EVENLY SPACED NODES. C C WN :IN THE WEIGHTS FOR SIMPSON'S RULE ON (A, B), WITH N EVENLY C SPACED NODES. C C XN :IN THE NYSTROM INTERPOLATES. C C N :IN THE NUMBER OF SUBINTERVALS ON (A, B). C C KERNEL :IN THIS IS A DOUBLE PRECISION FUNCTION OF TWO VARIABLES. IT C MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM C CALLING DIESMP. C C RHFCN :IN THIS IS A DOUBLE PRECISION FUNCTION OF ONE VARIABLE. IT C MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM C CALLING DIESMP. C C RHS :OUT VALUES OF RHFCN AT TM(I). C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880901 UPGRADE TO FORTRAN 77. C 901001 MODIFIED TO MEET SLATEC PROLOGUE STANDARS AND TO CORRECT C SMALL ERROR. C***END PROLOGUE DINTRP C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL DIMENSION TM(M),XM(M),TN(N),WN(N),XN(N),RHS(M) C C***FIRST EXECUTABLE STATEMENT DINTRP C C CALCULATE NYSTROM INTERPOLATING FORMULA. DO 50 I=1,M RHS(I)=RHFCN(TM(I)) XM(I)=RHS(I) DO 50 J=1,N 50 XM(I)=XM(I)+WN(J)*KERNEL(TM(I),TN(J))*XN(J) RETURN END DOUBLE PRECISION FUNCTION DNORM (X, Y, N, IFLAG) C***BEGIN PROLOGUE DNORM C***SUBSIDIARY C***PURPOSE FIND THE MAXIMUM NORM OF A VECTOR. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DNORM-D) C***AUTHOR ATKINSON, KENDALL, (UNIVERSITY OF IOWA) C ADDRESS: DEPARTMENT OF MATHEMATICS, C UNIVERSITY OF IOWA, IOWA CITY, IOWA 52242. C E_MAIL: BLAKEAPD@UIAMVS.BITNET C (ASSISTED BY DAVID CHIEN, GRADUATE STUDENT, U. OF IOWA) C***DESCRIPTION C C *USAGE: C C INTEGER N, IFLAG C DOUBLE PRECISION X(N), Y(N) C C VALUE = DNORM(X, Y, N, IFLAG) C C *ARGUMENT: C C X :IN IS A VECTOR WITH LENGTH N. C C Y :IN IS A VECTOR WITH LENGTH N. C C N :IN IS THE LENGTH OF VECTORS A AND B. C C IFLAG :IN IFLAG = 0: C CALCULATE THE MAXIMUM NORM OF X. C IFLAG = 1: C CALCULATE THE MAXIMUM NORM OF X-Y. C C *FUNCTION RETURN VALUES: C C VALUE : IFLAG = 0: C VALUE IS THE MAXIMUM NORM OF VECTOR X. C IFLAG = 1: C VALUE IS THE MAXIMUM NORM OF VECTOR X-Y. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880901 UPGRADE TO FORTRAN 77. C 901001 MODIFIED TO MEET SLATEC PROLOGUE STANDARS AND TO CORRECT C SMALL ERROR. C***END PROLOGUE DNORM C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (ZERO=0.0D0) DIMENSION X(N),Y(N) C C***FIRST EXECUTABLE STATEMENT DNORM IF (IFLAG .EQ. 0) THEN C FIND THE NORM OF X. DNORM=ZERO DO 10 I=1,N 10 DNORM=MAX(DNORM,ABS(X(I))) ELSE C FIND THE NORM OF X-Y. DNORM=ZERO DO 30 I=1,N 30 DNORM=MAX(DNORM,ABS(X(I)-Y(I))) END IF RETURN END SUBROUTINE DWANDT (W, T, N, A, B) C***BEGIN PROLOGUE DWANDT C***SUBSIDIARY C***PURPOSE THIS ROUTINE CALCULATES THE WEIGHTS AND NODE POINTS FOR C SIMPSON'S RULE ON (A,B), WITH N EVENLY SPACED NODES. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DWANDT-D) C***AUTHOR ATKINSON, KENDALL, (UNIVERSITY OF IOWA) C ADDRESS: DEPARTMENT OF MATHEMATICS, C UNIVERSITY OF IOWA, IOWA CITY, IOWA 52242. C E_MAIL: BLAKEAPD@UIAMVS.BITNET C (ASSISTED BY DAVID CHIEN, GRADUATE STUDENT, U. OF IOWA) C***DESCRIPTION C C *USAGE: C C INTEGER N C DOUBLE PRECISION W(N), T(N), A, B C PARAMETER(R3 = 1.0D0/3.0D0, TWOR3 = 2.0D0*R3, FOURR3 = 4.0D0*R3) C C CALL DWANDT (W, T, N, A, B) C C *ARGUMENT: C C W :OUT THE WEIGHTS FOR SIMPSON'S RULE ON (A, B), WITH N EVENLY C SPACED NODES. C C T :OUT THE NODE POINTS FOR SIMPSON'S RULE ON (A, B), WITH N EVENLY C SPACED NODES. C C N :IN IS THE NUMBER OF SUBINTERVALS ON (A, B). C C A :IN IS THE LOWER INTEGRATION LIMIT. C C B :IN IS THE UPPER INTEGRATION LIMIT. C C***SEE ALSO DIESMP C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 880901 UPGRADE TO FORTRAN 77. C 901001 MODIFIED TO MEET SLATEC PROLOGUE STANDARS AND TO CORRECT C SMALL ERROR. C***END PROLOGUE DWANDT IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(R3=1.0D0/3.0D0,TWOR3=2.0D0*R3,FOURR3=4.0D0*R3) DIMENSION W(N),T(N) C C***FIRST EXECUTABLE STATEMENT DWANDT C CALCULATE NODE POINTS. H=(B-A)/(N-1) DO 10 I=1,N 10 T(I)=A+(I-1)*H C C CALCULATE WEIGHTS. W(1)=R3*H W(N)=R3*H TEMP=FOURR3*H DO 20 I=2,N-1,2 20 W(I)=TEMP IF(N .EQ. 3) RETURN TEMP=TWOR3*H DO 30 I=3,N-2,2 30 W(I)=TEMP RETURN END SUBROUTINE LINSYS(A,LUFACT,N,B,SOLN,OPTION,RCOND,IPIVOT, * MACHIN,Z,ERRMAX) C ----------------- C C SOLVE A*X=B. C C THE LINPACK SUBROUTINES DGECO AND DGESL ARE USED. C C FORMAL PARAMETERS: C C A THIS CONTAINS THE COEFFICIENT MATRIX. C LUFACT THIS CONTAINS (OR WILL CONTAIN) THE LU FACTORIZATION C OF THE MATRIX IN A. SEE OPTION BELOW. C N THE ORDER OF THE LINEAR SYSTEM. C B THE RIGHT SIDE OF THE SYSTEM TO BE SOLVED. C SOLN ON EXIT, THIS CONTAINS THE SOLUTION OF THE LINEAR C SYSTEM. C OPTION C OPTION=1: IT IS ASSUMED THAT THE ARRAYS LUFACT AND A ARE THE C SAME. THE MATRIX IS FACTORED, AND THEN THE LINEAR C SYSTEM IS SOLVED. THE ORIGINAL CONTENTS OF A ARE C DESTROYED. IF B AND SOLN ARE THE SAME ARRAY, THEN C THE ORIGINAL CONTENTS OF B ARE LOST. C OPTION=2: THE LU FACTORIZATION OF A IS STORED IN LUFACT, C AND THE ORIGINAL CONTENTS OF A ARE LEFT C UNDISTURBED. THE LINEAR SYSTEM IS SOLVED. C THE RESIDUAL OF THE SOLUTION IS CALCULATED, AND C THEN THE ERROR IS ESTIMATED BY SOLVING THE C RESIDUAL ERROR EQUATION. DO NOT LET B=SOLN C OR A=LUFACT WHEN CALLING LINSYS, SINCE THE C ORIGINAL VALUES IN A AND B ARE NEEDED IN THE C RESIDUAL CORRECTION PROCESS. C OPTION=3: PROCEED AS WITH OPTION 1, BUT THE LU FACTORI- C ZATION IS ASSUMED TO HAVE ALREADY BEEN STORED C IN LUFACT. C OPTION=4: PROCEED AS WITH OPTION 2, BUT THE LU FACTORI- C ZATION IS ASSUMED TO HAVE ALREADY BEEN STORED C IN LUFACT. C C RCOND ON EXIT, THIS WILL CONTAIN AN ESTIMATE OF THE C RECIPROCAL OF A CONDITION NUMBER OF THE SYSTEM. C IF RCOND=0, THEN THE MATRIX IS EXACTLY SINGULAR. C IPIVOT THIS IS A TEMPORARY INTEGER WORKING VECTOR, OF LENGTH N, C TO RECORD THE ROW INTERCHANGES IN THE FACTORIZATION C PROCESS. C MACHIN THIS IS THE ACTUAL NUMBER OF ROWS THAT A IS DIMENSIONED C AS HAVING IN THE CALLING PROGRAM. C Z THIS IS A TEMPORARY WORKING VECTOR OF LENGTH N. C ERRMAX THIS WILL BE SET TO AN ESTIMATE OF THE RELATIVE ERROR C IN THE SOLUTION SOLN. IT IS OBTAINED BY THE RESIDUAL C CORRECTION METHOD. SINCE EXTENDED PRECISION IS NOT C BEING USED TO CALCULATE THE RESIDUALS, IT IS UNLIKELY C TO BE ACCURATE IN MORE THAN ITS MAGNITUDE. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LUFACT INTEGER OPTION DIMENSION A(MACHIN,*),LUFACT(MACHIN,*),B(*),SOLN(*), * IPIVOT(*),Z(*) PARAMETER(ZERO=0.0D0) SAVE TRCOND C IF(OPTION .EQ. 2) THEN C COPY THE MATRIX A INTO THE MATRIX LUFACT. DO 10 I=1,N DO 10 J=1,N 10 LUFACT(I,J)=A(I,J) END IF C C CALCULATE THE LU FACTORIZATION OF LUFACT. IF(OPTION .LT. 3) THEN CALL DGECO(LUFACT,MACHIN,N,IPIVOT,TRCOND,Z) RCOND = TRCOND IF(TRCOND .EQ. ZERO) RETURN END IF C C SOLVE THE GIVEN LINEAR SYSTEM. DO 20 I=1,N 20 SOLN(I)=B(I) CALL DGESL(LUFACT,MACHIN,N,IPIVOT,SOLN,0) IF((OPTION .EQ. 1) .OR. (OPTION .EQ. 3)) RETURN C C CALCULATE THE RESIDUAL OF THE LINEAR SYSTEM. DO 40 I=1,N SUM=ZERO DO 30 J=1,N 30 SUM=SUM+A(I,J)*SOLN(J) 40 Z(I)=B(I)-SUM R_MAX=ZERO DO 50 I=1,N 50 R_MAX=MAX(R_MAX,ABS(Z(I))) C C ESTIMATE THE ERROR IN THE SOLUTION OBTAINED ABOVE. CALL DGESL(LUFACT,MACHIN,N,IPIVOT,Z,0) ERRMAX=ZERO X_NORM=ZERO DO 60 I=1,N ERRMAX=MAX(ERRMAX,ABS(Z(I))) 60 X_NORM=MAX(X_NORM,ABS(SOLN(I))) ERRMAX=ERRMAX/X_NORM RETURN END SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) 00000010 Cc ---------------- INTEGER LDA,N,IPVT(*) 00000020 DOUBLE PRECISION A(LDA,*),Z(*) 00000030 DOUBLE PRECISION RCOND 00000040 C 00000050 C DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION 00000060 C AND ESTIMATES THE CONDITION OF THE MATRIX. 00000070 C 00000080 C IF RCOND IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER. 00000090 C TO SOLVE A*X = B , FOLLOW DGECO BY DGESL. 00000100 C TO COMPUTE INVERSE(A)*C , FOLLOW DGECO BY DGESL. 00000110 C TO COMPUTE DETERMINANT(A) , FOLLOW DGECO BY DGEDI. 00000120 C TO COMPUTE INVERSE(A) , FOLLOW DGECO BY DGEDI. 00000130 C 00000140 C ON ENTRY 00000150 C 00000160 C A DOUBLE PRECISION(LDA, N) 00000170 C THE MATRIX TO BE FACTORED. 00000180 C 00000190 C LDA INTEGER 00000200 C THE LEADING DIMENSION OF THE ARRAY A . 00000210 C 00000220 C N INTEGER 00000230 C THE ORDER OF THE MATRIX A . 00000240 C 00000250 C ON RETURN 00000260 C 00000270 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS 00000280 C WHICH WERE USED TO OBTAIN IT. 00000290 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE 00000300 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER 00000310 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. 00000320 C 00000330 C IPVT INTEGER(N) 00000340 C AN INTEGER VECTOR OF PIVOT INDICES. 00000350 C 00000360 C RCOND DOUBLE PRECISION 00000370 C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . 00000380 C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS 00000390 C IN A AND B OF SIZE EPSILON MAY CAUSE 00000400 C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . 00000410 C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION 00000420 C 1.0 + RCOND .EQ. 1.0 00000430 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING 00000440 C PRECISION. IN PARTICULAR, RCOND IS ZERO IF 00000450 C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE 00000460 C UNDERFLOWS. 00000470 C 00000480 C Z DOUBLE PRECISION(N) 00000490 C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. 00000500 C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS 00000510 C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT 00000520 C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . 00000530 C 00000540 C LINPACK. THIS VERSION DATED 08/14/78 . 00000550 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. 00000560 C 00000570 C SUBROUTINES AND FUNCTIONS 00000580 C 00000590 C LINPACK DGEFA 00000600 C BLAS DAXPY,DDOT,DSCAL,DASUM 00000610 C FORTRAN DABS,DMAX1,DSIGN 00000620 C 00000630 C INTERNAL VARIABLES 00000640 C 00000650 DOUBLE PRECISION DDOT,EK,T,WK,WKM 00000660 DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM 00000670 INTEGER INFO,J,K,KB,KP1,L 00000680 C 00000690 C 00000700 C COMPUTE 1-NORM OF A 00000710 C 00000720 ANORM = 0.0D0 00000730 DO 10 J = 1, N 00000740 ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1)) 00000750 10 CONTINUE 00000760 C 00000770 C FACTOR 00000780 C 00000790 CALL DGEFA(A,LDA,N,IPVT,INFO) 00000800 C 00000810 C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . 00000820 C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . 00000830 C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE 00000840 C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE 00000850 C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID 00000860 C OVERFLOW. 00000870 C 00000880 C SOLVE TRANS(U)*W = E 00000890 C 00000900 EK = 1.0D0 00000910 DO 20 J = 1, N 00000920 Z(J) = 0.0D0 00000930 20 CONTINUE 00000940 DO 100 K = 1, N 00000950 IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) 00000960 IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 00000970 S = DABS(A(K,K))/DABS(EK-Z(K)) 00000980 CALL DSCAL(N,S,Z,1) 00000990 EK = S*EK 00001000 30 CONTINUE 00001010 WK = EK - Z(K) 00001020 WKM = -EK - Z(K) 00001030 S = DABS(WK) 00001040 SM = DABS(WKM) 00001050 IF (A(K,K) .EQ. 0.0D0) GO TO 40 00001060 WK = WK/A(K,K) 00001070 WKM = WKM/A(K,K) 00001080 GO TO 50 00001090 40 CONTINUE 00001100 WK = 1.0D0 00001110 WKM = 1.0D0 00001120 50 CONTINUE 00001130 KP1 = K + 1 00001140 IF (KP1 .GT. N) GO TO 90 00001150 DO 60 J = KP1, N 00001160 SM = SM + DABS(Z(J)+WKM*A(K,J)) 00001170 Z(J) = Z(J) + WK*A(K,J) 00001180 S = S + DABS(Z(J)) 00001190 60 CONTINUE 00001200 IF (S .GE. SM) GO TO 80 00001210 T = WKM - WK 00001220 WK = WKM 00001230 DO 70 J = KP1, N 00001240 Z(J) = Z(J) + T*A(K,J) 00001250 70 CONTINUE 00001260 80 CONTINUE 00001270 90 CONTINUE 00001280 Z(K) = WK 00001290 100 CONTINUE 00001300 S = 1.0D0/DASUM(N,Z,1) 00001310 CALL DSCAL(N,S,Z,1) 00001320 C 00001330 C SOLVE TRANS(L)*Y = W 00001340 C 00001350 DO 120 KB = 1, N 00001360 K = N + 1 - KB 00001370 IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) 00001380 IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 00001390 S = 1.0D0/DABS(Z(K)) 00001400 CALL DSCAL(N,S,Z,1) 00001410 110 CONTINUE 00001420 L = IPVT(K) 00001430 T = Z(L) 00001440 Z(L) = Z(K) 00001450 Z(K) = T 00001460 120 CONTINUE 00001470 S = 1.0D0/DASUM(N,Z,1) 00001480 CALL DSCAL(N,S,Z,1) 00001490 C 00001500 YNORM = 1.0D0 00001510 C 00001520 C SOLVE L*V = Y 00001530 C 00001540 DO 140 K = 1, N 00001550 L = IPVT(K) 00001560 T = Z(L) 00001570 Z(L) = Z(K) 00001580 Z(K) = T 00001590 IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) 00001600 IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 00001610 S = 1.0D0/DABS(Z(K)) 00001620 CALL DSCAL(N,S,Z,1) 00001630 YNORM = S*YNORM 00001640 130 CONTINUE 00001650 140 CONTINUE 00001660 S = 1.0D0/DASUM(N,Z,1) 00001670 CALL DSCAL(N,S,Z,1) 00001680 YNORM = S*YNORM 00001690 C 00001700 C SOLVE U*Z = V 00001710 C 00001720 DO 160 KB = 1, N 00001730 K = N + 1 - KB 00001740 IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 00001750 S = DABS(A(K,K))/DABS(Z(K)) 00001760 CALL DSCAL(N,S,Z,1) 00001770 YNORM = S*YNORM 00001780 150 CONTINUE 00001790 IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) 00001800 IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 00001810 T = -Z(K) 00001820 CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) 00001830 160 CONTINUE 00001840 C MAKE ZNORM = 1.0 00001850 S = 1.0D0/DASUM(N,Z,1) 00001860 CALL DSCAL(N,S,Z,1) 00001870 YNORM = S*YNORM 00001880 C 00001890 IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM 00001900 IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 00001910 RETURN 00001920 END 00001930 SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) 00000010 C ---------------- C INTEGER LDA,N,IPVT(*),JOB 00000020 DOUBLE PRECISION A(LDA,*),B(*) 00000030 C 00000040 C DGESL SOLVES THE DOUBLE PRECISION SYSTEM 00000050 C A * X = B OR TRANS(A) * X = B 00000060 C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. 00000070 C 00000080 C ON ENTRY 00000090 C 00000100 C A DOUBLE PRECISION(LDA, N) 00000110 C THE OUTPUT FROM DGECO OR DGEFA. 00000120 C 00000130 C LDA INTEGER 00000140 C THE LEADING DIMENSION OF THE ARRAY A . 00000150 C 00000160 C N INTEGER 00000170 C THE ORDER OF THE MATRIX A . 00000180 C 00000190 C IPVT INTEGER(N) 00000200 C THE PIVOT VECTOR FROM DGECO OR DGEFA. 00000210 C 00000220 C B DOUBLE PRECISION(N) 00000230 C THE RIGHT HAND SIDE VECTOR. 00000240 C 00000250 C JOB INTEGER 00000260 C = 0 TO SOLVE A*X = B , 00000270 C = NONZERO TO SOLVE TRANS(A)*X = B WHERE 00000280 C TRANS(A) IS THE TRANSPOSE. 00000290 C 00000300 C ON RETURN 00000310 C 00000320 C B THE SOLUTION VECTOR X . 00000330 C 00000340 C ERROR CONDITION 00000350 C 00000360 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A 00000370 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY 00000380 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER 00000390 C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE 00000400 C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 00000410 C OR DGEFA HAS SET INFO .EQ. 0 . 00000420 C 00000430 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX 00000440 C WITH P COLUMNS 00000450 C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) 00000460 C IF (RCOND IS TOO SMALL) GO TO ... 00000470 C DO 10 J = 1, P 00000480 C CALL DGESL(A,LDA,N,IPVT,C(1,J),0) 00000490 C 10 CONTINUE 00000500 C 00000510 C LINPACK. THIS VERSION DATED 08/14/78 . 00000520 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. 00000530 C 00000540 C SUBROUTINES AND FUNCTIONS 00000550 C 00000560 C BLAS DAXPY,DDOT 00000570 C 00000580 C INTERNAL VARIABLES 00000590 C 00000600 DOUBLE PRECISION DDOT,T 00000610 INTEGER K,KB,L,NM1 00000620 C 00000630 NM1 = N - 1 00000640 IF (JOB .NE. 0) GO TO 50 00000650 C 00000660 C JOB = 0 , SOLVE A * X = B 00000670 C FIRST SOLVE L*Y = B 00000680 C 00000690 IF (NM1 .LT. 1) GO TO 30 00000700 DO 20 K = 1, NM1 00000710 L = IPVT(K) 00000720 T = B(L) 00000730 IF (L .EQ. K) GO TO 10 00000740 B(L) = B(K) 00000750 B(K) = T 00000760 10 CONTINUE 00000770 CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 00000780 20 CONTINUE 00000790 30 CONTINUE 00000800 C 00000810 C NOW SOLVE U*X = Y 00000820 C 00000830 DO 40 KB = 1, N 00000840 K = N + 1 - KB 00000850 B(K) = B(K)/A(K,K) 00000860 T = -B(K) 00000870 CALL DAXPY(K-1,T,A(1,K),1,B(1),1) 00000880 40 CONTINUE 00000890 GO TO 100 00000900 50 CONTINUE 00000910 C 00000920 C JOB = NONZERO, SOLVE TRANS(A) * X = B 00000930 C FIRST SOLVE TRANS(U)*Y = B 00000940 C 00000950 DO 60 K = 1, N 00000960 T = DDOT(K-1,A(1,K),1,B(1),1) 00000970 B(K) = (B(K) - T)/A(K,K) 00000980 60 CONTINUE 00000990 C 00001000 C NOW SOLVE TRANS(L)*X = Y 00001010 C 00001020 IF (NM1 .LT. 1) GO TO 90 00001030 DO 80 KB = 1, NM1 00001040 K = N - KB 00001050 B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) 00001060 L = IPVT(K) 00001070 IF (L .EQ. K) GO TO 70 00001080 T = B(L) 00001090 B(L) = B(K) 00001100 B(K) = T 00001110 70 CONTINUE 00001120 80 CONTINUE 00001130 90 CONTINUE 00001140 100 CONTINUE 00001150 RETURN 00001160 END 00001170 SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) 00000010 C ---------------- C INTEGER LDA,N,IPVT(*),INFO 00000020 DOUBLE PRECISION A(LDA,*) 00000030 C 00000040 C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. 00000050 C 00000060 C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED 00000070 C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. 00000080 C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) . 00000090 C 00000100 C ON ENTRY 00000110 C 00000120 C A DOUBLE PRECISION(LDA, N) 00000130 C THE MATRIX TO BE FACTORED. 00000140 C 00000150 C LDA INTEGER 00000160 C THE LEADING DIMENSION OF THE ARRAY A . 00000170 C 00000180 C N INTEGER 00000190 C THE ORDER OF THE MATRIX A . 00000200 C 00000210 C ON RETURN 00000220 C 00000230 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS 00000240 C WHICH WERE USED TO OBTAIN IT. 00000250 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE 00000260 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER 00000270 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. 00000280 C 00000290 C IPVT INTEGER(N) 00000300 C AN INTEGER VECTOR OF PIVOT INDICES. 00000310 C 00000320 C INFO INTEGER 00000330 C = 0 NORMAL VALUE. 00000340 C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR 00000350 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES 00000360 C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO 00000370 C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE 00000380 C INDICATION OF SINGULARITY. 00000390 C 00000400 C LINPACK. THIS VERSION DATED 08/14/78 . 00000410 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. 00000420 C 00000430 C SUBROUTINES AND FUNCTIONS 00000440 C 00000450 C BLAS DAXPY,DSCAL,IDAMAX 00000460 C 00000470 C INTERNAL VARIABLES 00000480 C 00000490 DOUBLE PRECISION T 00000500 INTEGER IDAMAX,J,K,KP1,L,NM1 00000510 C 00000520 C 00000530 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 00000540 C 00000550 INFO = 0 00000560 NM1 = N - 1 00000570 IF (NM1 .LT. 1) GO TO 70 00000580 DO 60 K = 1, NM1 00000590 KP1 = K + 1 00000600 C 00000610 C FIND L = PIVOT INDEX 00000620 C 00000630 L = IDAMAX(N-K+1,A(K,K),1) + K - 1 00000640 IPVT(K) = L 00000650 C 00000660 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED 00000670 C 00000680 IF (A(L,K) .EQ. 0.0D0) GO TO 40 00000690 C 00000700 C INTERCHANGE IF NECESSARY 00000710 C 00000720 IF (L .EQ. K) GO TO 10 00000730 T = A(L,K) 00000740 A(L,K) = A(K,K) 00000750 A(K,K) = T 00000760 10 CONTINUE 00000770 C 00000780 C COMPUTE MULTIPLIERS 00000790 C 00000800 T = -1.0D0/A(K,K) 00000810 CALL DSCAL(N-K,T,A(K+1,K),1) 00000820 C 00000830 C ROW ELIMINATION WITH COLUMN INDEXING 00000840 C 00000850 DO 30 J = KP1, N 00000860 T = A(L,J) 00000870 IF (L .EQ. K) GO TO 20 00000880 A(L,J) = A(K,J) 00000890 A(K,J) = T 00000900 20 CONTINUE 00000910 CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 00000920 30 CONTINUE 00000930 GO TO 50 00000940 40 CONTINUE 00000950 INFO = K 00000960 50 CONTINUE 00000970 60 CONTINUE 00000980 70 CONTINUE 00000990 IPVT(N) = N 00001000 IF (A(N,N) .EQ. 0.0D0) INFO = N 00001010 RETURN 00001020 END 00001030 SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) 00000010 C ---------------- C 00000020 C CONSTANT TIMES A VECTOR PLUS A VECTOR. 00000030 C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. 00000040 C JACK DONGARRA, LINPACK, 3/11/78. 00000050 C 00000060 DOUBLE PRECISION DX(*),DY(*),DA 00000070 INTEGER I,INCX,INCY,IXIY,M,MP1,N 00000080 C 00000090 IF(N.LE.0)RETURN 00000100 IF (DA .EQ. 0.0D0) RETURN 00000110 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 00000120 C 00000130 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 00000140 C NOT EQUAL TO 1 00000150 C 00000160 IX = 1 00000170 IY = 1 00000180 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 00000190 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 00000200 DO 10 I = 1,N 00000210 DY(IY) = DY(IY) + DA*DX(IX) 00000220 IX = IX + INCX 00000230 IY = IY + INCY 00000240 10 CONTINUE 00000250 RETURN 00000260 C 00000270 C CODE FOR BOTH INCREMENTS EQUAL TO 1 00000280 C 00000290 C 00000300 C CLEAN-UP LOOP 00000310 C 00000320 20 M = MOD(N,4) 00000330 IF( M .EQ. 0 ) GO TO 40 00000340 DO 30 I = 1,M 00000350 DY(I) = DY(I) + DA*DX(I) 00000360 30 CONTINUE 00000370 IF( N .LT. 4 ) RETURN 00000380 40 MP1 = M + 1 00000390 DO 50 I = MP1,N,4 00000400 DY(I) = DY(I) + DA*DX(I) 00000410 DY(I + 1) = DY(I + 1) + DA*DX(I + 1) 00000420 DY(I + 2) = DY(I + 2) + DA*DX(I + 2) 00000430 DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 00000440 50 CONTINUE 00000450 RETURN 00000460 END 00000470 DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) 00000010 C ------------------------------ C 00000020 C FORMS THE DOT PRODUCT OF TWO VECTORS. 00000030 C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. 00000040 C JACK DONGARRA, LINPACK, 3/11/78. 00000050 C 00000060 DOUBLE PRECISION DX(*),DY(*),DTEMP 00000070 INTEGER I,INCX,INCY,IX,IY,M,MP1,N 00000080 C 00000090 DDOT = 0.0D0 00000100 DTEMP = 0.0D0 00000110 IF(N.LE.0)RETURN 00000120 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 00000130 C 00000140 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 00000150 C NOT EQUAL TO 1 00000160 C 00000170 IX = 1 00000180 IY = 1 00000190 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 00000200 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 00000210 DO 10 I = 1,N 00000220 DTEMP = DTEMP + DX(IX)*DY(IY) 00000230 IX = IX + INCX 00000240 IY = IY + INCY 00000250 10 CONTINUE 00000260 DDOT = DTEMP 00000270 RETURN 00000280 C 00000290 C CODE FOR BOTH INCREMENTS EQUAL TO 1 00000300 C 00000310 C 00000320 C CLEAN-UP LOOP 00000330 C 00000340 20 M = MOD(N,5) 00000350 IF( M .EQ. 0 ) GO TO 40 00000360 DO 30 I = 1,M 00000370 DTEMP = DTEMP + DX(I)*DY(I) 00000380 30 CONTINUE 00000390 IF( N .LT. 5 ) GO TO 60 00000400 40 MP1 = M + 1 00000410 DO 50 I = MP1,N,5 00000420 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + 00000430 * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)00000440 50 CONTINUE 00000450 60 DDOT = DTEMP 00000460 RETURN 00000470 END 00000480 SUBROUTINE DSCAL(N,DA,DX,INCX) 00000010 C ----------------- C 00000020 C SCALES A VECTOR BY A CONSTANT. 00000030 C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. 00000040 C JACK DONGARRA, LINPACK, 3/11/78. 00000050 C 00000060 DOUBLE PRECISION DA,DX(*) 00000070 INTEGER I,INCX,M,MP1,N,NINCX 00000080 C 00000090 IF(N.LE.0)RETURN 00000100 IF(INCX.EQ.1)GO TO 20 00000110 C 00000120 C CODE FOR INCREMENT NOT EQUAL TO 1 00000130 C 00000140 NINCX = N*INCX 00000150 DO 10 I = 1,NINCX,INCX 00000160 DX(I) = DA*DX(I) 00000170 10 CONTINUE 00000180 RETURN 00000190 C 00000200 C CODE FOR INCREMENT EQUAL TO 1 00000210 C 00000220 C 00000230 C CLEAN-UP LOOP 00000240 C 00000250 20 M = MOD(N,5) 00000260 IF( M .EQ. 0 ) GO TO 40 00000270 DO 30 I = 1,M 00000280 DX(I) = DA*DX(I) 00000290 30 CONTINUE 00000300 IF( N .LT. 5 ) RETURN 00000310 40 MP1 = M + 1 00000320 DO 50 I = MP1,N,5 00000330 DX(I) = DA*DX(I) 00000340 DX(I + 1) = DA*DX(I + 1) 00000350 DX(I + 2) = DA*DX(I + 2) 00000360 DX(I + 3) = DA*DX(I + 3) 00000370 DX(I + 4) = DA*DX(I + 4) 00000380 50 CONTINUE 00000390 RETURN 00000400 END 00000410 DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) 00000010 C ------------------------------- C 00000020 C TAKES THE SUM OF THE ABSOLUTE VALUES. 00000030 C JACK DONGARRA, LINPACK, 3/11/78. 00000040 C 00000050 DOUBLE PRECISION DX(*),DTEMP 00000060 INTEGER I,INCX,M,MP1,N,NINCX 00000070 C 00000080 DASUM = 0.0D0 00000090 DTEMP = 0.0D0 00000100 IF(N.LE.0)RETURN 00000110 IF(INCX.EQ.1)GO TO 20 00000120 C 00000130 C CODE FOR INCREMENT NOT EQUAL TO 1 00000140 C 00000150 NINCX = N*INCX 00000160 DO 10 I = 1,NINCX,INCX 00000170 DTEMP = DTEMP + DABS(DX(I)) 00000180 10 CONTINUE 00000190 DASUM = DTEMP 00000200 RETURN 00000210 C 00000220 C CODE FOR INCREMENT EQUAL TO 1 00000230 C 00000240 C 00000250 C CLEAN-UP LOOP 00000260 C 00000270 20 M = MOD(N,6) 00000280 IF( M .EQ. 0 ) GO TO 40 00000290 DO 30 I = 1,M 00000300 DTEMP = DTEMP + DABS(DX(I)) 00000310 30 CONTINUE 00000320 IF( N .LT. 6 ) GO TO 60 00000330 40 MP1 = M + 1 00000340 DO 50 I = MP1,N,6 00000350 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) 00000360 * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 00000370 50 CONTINUE 00000380 60 DASUM = DTEMP 00000390 RETURN 00000400 END 00000410 INTEGER FUNCTION IDAMAX(N,DX,INCX) 00000010 C ----------------------- C 00000020 C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. 00000030 C JACK DONGARRA, LINPACK, 3/11/78. 00000040 C 00000050 DOUBLE PRECISION DX(*),DMAX 00000060 INTEGER I,INCX,IX,N 00000070 C 00000080 IDAMAX = 0 00000090 IF( N .LT. 1 ) RETURN 00000100 IDAMAX = 1 00000110 IF(N.EQ.1)RETURN 00000120 IF(INCX.EQ.1)GO TO 20 00000130 C 00000140 C CODE FOR INCREMENT NOT EQUAL TO 1 00000150 C 00000160 IX = 1 00000170 DMAX = DABS(DX(1)) 00000180 IX = IX + INCX 00000190 DO 10 I = 2,N 00000200 IF(DABS(DX(IX)).LE.DMAX) GO TO 5 00000210 IDAMAX = I 00000220 DMAX = DABS(DX(IX)) 00000230 5 IX = IX + INCX 00000240 10 CONTINUE 00000250 RETURN 00000260 C 00000270 C CODE FOR INCREMENT EQUAL TO 1 00000280 C 00000290 20 DMAX = DABS(DX(1)) 00000300 DO 30 I = 2,N 00000310 IF(DABS(DX(I)).LE.DMAX) GO TO 30 00000320 IDAMAX = I 00000330 DMAX = DABS(DX(I)) 00000340 30 CONTINUE 00000350 RETURN 00000360 END 00000370