C TITLE: TEST OF DIEGAU: C --------------------- C C THIS IS A TEST PROGRAM FOR THE INTEGRAL EQUATION PROGRAM DIEGAU. C 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 C FUNCTIONS 'FCNY' WHICH CAN BE CALLED UPON BY SETTING 'NUMKER' C AND 'NUMFNY', RESPECTIVELY. THE FUNCTIONS FCNK AND FCNY USUALLY C INVOLVE SOME PARAMETERS, CALLED 'CONST' AND 'AUXCON' IN THE C BELOW INPUT STATEMENTS. THERE IS ALSO A PARAMETER 'LAMBDA', C WHICH IS MULTIPLIED TIMES THE INTEGRAL 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). 'NT' CONTROLS THE CHOICE OF OUTPUT NODES FROM C THE INTEGRAL EQUATION. IF NT=0, THEN THE OUTPUT VALUES ARE C EVALUATED AT APPROPRIATE GAUSSIAN QUADRATURE NODE POINTS. IF C NT IS GREATER THAN 0, THEN THE INTERVAL OF INTEGRATION IS C DIVIDED INTO AN EVEN MESH T(I)=A+(I-1)*(B-A)/(NT-1), I=1,...,NT. C SEE THE INTRODUCTORY COMMENTS IN DIEGAU FOR MORE INFORMATION C ON THE MEANING OF THESE AND OTHER VARIABLES. C THE PROGRAM 'DIEGAU' USES THE LINPACK SUBROUTINE 'DGECO' AND C 'DGESL'. THEIR USE IS EMBEDDED IN THE SUBROUTINE 'LINSYS', THAT C IS USED IN SUBROUTINES 'DIEGAU' AND 'DITERT'. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER TMLOUT, TML_IN DOUBLE PRECISION LAMBDA,NORMX,INFO CHARACTER Q PARAMETER(MAXN=32, MAXM=256, MAXWRK=5*MAXN**2+10*MAXN+9*MAXM) DIMENSION X(MAXM),T(MAXM),WORK(MAXWRK),IWORK(MAXN),INFO(6) PARAMETER (ZERO=0.0D0) EXTERNAL FCNK,FCNY COMMON/PARAM/LAMBDA,CONST,AUXCON,NUMKER,NUMFNY C C********************************************************************** C * C FOLLOWING IS THE OUTPUT UNIT NUMBER TO BE USED IN THE OUTPUT * C OF THE RESULTS TO THE FILE 'IEGAUS.OUT' * C * IWRITE = 8 C * C FOLLOWING ARE THE UNIT NUMBERS FOR THE STANDARD INPUT AND * C OUTPUT UNITS. * C * TML_IN = 5 TMLOUT = 6 C * C********************************************************************** C C READ OUTPUT DIRECTIONS PRINT *, 'DO YOU WANT THE OUTPUT ON THE TERMINAL? (Y/N)' READ(*,'(A1)') Q IF((Q .EQ. 'Y') .OR. (Q .EQ. 'y')) THEN IOUT = TMLOUT ELSE OPEN(IWRITE,FILE='IEGAUS.OUT') IOUT = IWRITE END IF C C READ INPUT PARAMETERS AND INITIALIZE TEST. 10 PRINT *, ' GIVE NUMKER AND NUMFNY.' READ *, NUMKER,NUMFNY PRINT *, ' GIVE THE INTERVAL [A,B].' READ *, A, B PRINT *, ' GIVE LAMBDA, CONST, AND EPS.' READ *, LAMBDA, CONST, EPS PRINT *, ' GIVE NT, AUXCON, AND IFLAG.' READ *, NT, AUXCON, IFLAG C C PRINT THE INPUT VARIABLES. WRITE(IOUT,1001) NUMKER,NUMFNY,LAMBDA,CONST,AUXCON,A,B, * EPS,NT,IFLAG C C IF NECESSARY, PRODUCE NODES FOR EVALUATING UNKNOWN FUNCTION. IF(NT .GT. 0) THEN C INITIALIZE T. H=(B-A)/(NT-1) DO 20 I=1,NT 20 T(I)=A+(I-1)*H END IF C CALL IEGAUS TO SOLVE INTEGRAL EQUATION. CALL DIEGAU(FCNK,FCNY,A,B,EPS,IFLAG,X,T,NT,MAXN,MAXM, * WORK,IWORK,IER,INFO) C NFINAL = INT(INFO(5)) MFINAL = INT(INFO(6)) C PRINT THE OUTPUT VALUES. WRITE(IOUT,1004) EPS,IER,NFINAL,MFINAL,INFO(1),INFO(2),INFO(3), * INFO(4) IF(IFLAG .EQ. 0) THEN WRITE(IOUT,1002) ELSE C IFLAG=1. WRITE(IOUT,1007) END IF C NORMX=ZERO DO 30 I=1,NT 30 NORMX=MAX(NORMX,ABS(X(I))) C C PRINT COMPUTED RESULTS, TRUE ANSWERS, AND TRUE ERRORS, C INTERPRETED ACCORDING TO IFLAG. ERRMAX=ZERO DO 40 I=1,NT ANS=TRUE(T(I)) ERROR=ANS-X(I) ERRMAX=MAX(ERRMAX,ABS(ERROR)) IF(IFLAG .EQ. 1) ERROR=ERROR/NORMX 40 WRITE(IOUT,1003) T(I),X(I),ANS,ERROR C WRITE(IOUT,1005) ERRMAX WRITE(IOUT,1006) ERRMAX/NORMX C C CHECK WHETHER TO INPUT MORE DATA. PRINT *, ' DO YOU WANT TO SOLVE ANOTHER EQUATION? (Y/N)' READ (*,'(A1)') Q IF((Q .EQ. 'Y') .OR. (Q .EQ. 'y')) THEN GO TO 10 ELSE IF(IOUT .EQ. IWRITE) CLOSE(IWRITE) STOP END IF C C FORMAT STATEMENTS FOR OUTPUT. 1001 FORMAT('1KERNEL NO.',I2,5X,'RHFCN NO.',I2,/,' LAMBDA=',1PD17.10, * /,5X,'CONSTANT=',D17.10,5X,'AUX. CONSTANT=',D17.10,/, * 5X,'[A,B] = ',D17.10,D20.10,/,5X,' DESIRED ERROR=',D8.2, * 5X,'NT=',I3,5X,' IFLAG=',I1,//) 1002 FORMAT(//,12X,'T',14X,'APPROX(T)',14X,'TRUE(T)',9X,'ERROR',/) 1003 FORMAT(1PD20.10,2D20.10,D12.2) 1004 FORMAT(' PREDICTED ERROR = ',1PD8.2,5X,'IER = ',I1,/, * ' FINAL N = ',I3,5X,'FINAL M = ',I4,/,' FINAL R1 = ',D8.2,5X, * 'FINAL R2 = ',D8.2,/,' FINAL DESIRED ERROR = ',D8.2, * 5X,'NORM(K) = ',D8.2,//) 1005 FORMAT(//,' ACTUAL ERROR = ',1PD10.2) 1006 FORMAT(' ACTUAL RELATIVE ERROR = ',1PD8.2) 1007 FORMAT(//,12X,'T',14X,'APPROX(T)',14X,'TRUE(T)',8X,'REL ERR',/) END DOUBLE PRECISION FUNCTION FCNK(S,T) C ------------------------------ C C DEFINE THE KERNEL FUNCTION OF THE INTEGRAL EQUATION. C THE PARAMETERS C AND ALPHA MAY BE USED IN DEFINING K(S,T). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LAMBDA PARAMETER (ONE=1.0D0, TWO=2.0D0, FIVE=5.0D0) COMMON/PARAM/LAMBDA,C,ALPHA,NUMKER,NUMFNY DATA PI/3.14159265358979D0/,CI/2202.5465794806717D0/ C GO TO (10,20,30,40,50,60,70,80) NUMKER C C USE WITH NUMFNY=1. 10 FCNK = LAMBDA*COS(PI*S*T) RETURN C C USE WITH NUMFNY=5. 20 FCNK = C/(C*C+(S-T)*(S-T))*LAMBDA RETURN C C USE WITH NUMFNY=2, 3, OR 10. 30 C2=C*C FCNK = (ONE-C2)/(ONE+C2-TWO*C*COS(TWO*PI*(S+T)))*LAMBDA RETURN C C USE WITH NUMFNY=6. 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=7 OR 8. 50 FCNK = LAMBDA/(S+S+T+C) RETURN C C USE WITH NUMFNY=9. 60 N = C FCNK = LAMBDA*(T-S)**N RETURN C C USE WITH NUMFNY=4. 70 BETA = C FCNK = LAMBDA*EXP(BETA*S*T) RETURN C C USE WITH NUMFNY=12. 80 FCNK = LAMBDA*EXP(FIVE*(S+T))/CI RETURN END DOUBLE PRECISION FUNCTION FCNY(S) C ------------------------------ C C DEFINE THE RIGHT SIDE OF THE INTEGRAL EQUATION. C THE PARAMETERS C AND ALPHA MAY BE USED IN DEFINING Y(S). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LAMBDA PARAMETER (ZERO=0.0D0, P06=.06D0, P2=.2D0, P5=.5D0, * P8=.8D0, ONE=1.0D0, TWO=2.0D0, THREE=3.0D0, * SIX=6.0D0, SEVEN=7.0D0) COMMON/PARAM/LAMBDA,C,ALPHA,NUMKER,NUMFNY DATA PI/3.14159265358979D0/ DATA C1/.4472135955D0/,C2/.8944271910D0/ C GO TO (10,20,30,40,50,60,70,80,90,100,110),NUMFNY C C USE WITH NUMKER=1. 10 T1 = SEVEN+PI*S T2 = SEVEN-PI*S EB = EXP(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=7. 40 BETA = C FCNY = EXP(ALPHA*S)-LAMBDA*(EXP(ALPHA+BETA*S)-ONE)/(ALPHA * +BETA*S) RETURN C C USE WITH NUMKER=2. 50 A = -P8 B = P06 CS = C*C T1 = LOG((CS+(ONE-S)**2)/(CS+S*S)) T2 = ATAN((ONE-S)/C)+ATAN(S/C) FCNY = C+C*(S+A/TWO)*T1+(B-CS+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 C C USE WITH NUMKER=5. 70 TEMP = SQRT(C+S+S) TEMP = TWO*(ONE-TEMP*ATAN(ONE/TEMP)) FCNY = SQRT(S)-LAMBDA*TEMP RETURN C C USE WITH NUMKER=5. 80 TEMP = SQRT(C+S+S+P2) T2 = P5*TEMP*LOG((TEMP+C1)/(TEMP-C1)) T3 = TEMP*ATAN(C2/TEMP) TEMP = TWO*(C2-C1+T2-T3) FCNY = SQRT(ABS(S-P2))-LAMBDA*TEMP RETURN C C USE WITH NUMKER=6. 90 N = C SUM = ONE/(C+TWO)**2 PROD = 1 DO 91 I=1,N FI = I PROD = -PROD*S*((C-FI+ONE)/FI) 91 SUM = SUM+PROD/(C+TWO-FI)**2 SUM = -SUM IF(S .NE. ZERO) THEN FCNY = S*LOG(S)-LAMBDA*SUM ELSE FCNY = -LAMBDA*SUM END IF RETURN C C USE WITH NUMKER=3. 100 FCNY = COS(TWO*PI*S) RETURN C C USE WITH NUMKER =8. 110 FCNY = S RETURN END DOUBLE PRECISION FUNCTION TRUE(S) C ------------------------------ C C THIS DEFINES THE TRUE SOLUTION X(S) OF THE INTEGRAL EQUATION. C THE PARAMETERS C AND ALPHA MAY BE USED IN DEFINING X(S). C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION LAMBDA PARAMETER (ZERO=0.0D0, P06=.06D0, P2=.2D0, * P8=.8D0, ONE=1.0D0, TWO=2.0D0, FIVE=5.0D0, * SIX=6.0D0, SEVEN=7.0D0) COMMON/PARAM/LAMBDA,C,ALPHA,NUMKER,NUMFNY DATA PI/3.14159265358979D0/ DATA DP/23.7861054564D0/,CI/2202.5465794806717D0/ C GO TO (10,20,30,40,50,60,70,80,90,100,110),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=7. 40 TRUE = EXP(ALPHA*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 C C USE WITH NUMKER=5. 70 TRUE = SQRT(S) RETURN C C USE WITH NUMKER=5. 80 TRUE = SQRT(ABS(S-P2)) RETURN C C USE WITH NUMKER=6. 90 IF(S .NE. ZERO) THEN TRUE = S*LOG(S) ELSE TRUE = ZERO END IF RETURN C C USE WITH NUMKER=3. 100 TRUE = COS(TWO*PI*S)/(ONE-LAMBDA*C) RETURN C C USE WITH NUMKER=8. 110 TEMP = LAMBDA*DP/(CI*(ONE-LAMBDA)) TRUE = S+TEMP*EXP(FIVE*S) RETURN END SUBROUTINE DIEGAU (KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, * NUPPER, MUPPER, W, IW, IER, INFO) C***BEGIN PROLOGUE DIEGAU 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 (DIEGAU-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 IFLAG, NT, NUPPER, MUPPER, IW(NUPPER), IER C DOUBLE PRECISION KERNEL, RHFCN, A, B, EP, X(MUP), T(MUP), W(*), C INFO(6) C EXTERNAL KERNEL, RHFCN C PARAMETER (P1 = .1D0, P5 = .5D0) C C CALL DIEGAU (KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, NUPPER, C MUPPER, W, IW, IER, INFO) 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 DIEGAU. C C RHFCN :EXT THIS IS A DOUBLE PRECISION FUNCTION OF ONE VARIABLE. IT C MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE C PROGRAM CALLING DIEGAU. C C A :IN IS THE LOWER INTEGRATION LIMIT. C C B :IN IS THE UPPER INTEGRATION LIMIT. C C EP :INOUT 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 C IFLAG :IN IFLAG = 0: C EP IS INTERPRETED AS AN ABSOLUTE ERROR TOLERANCE. C IFLAG = 1: C EP IS INTERPRETED AS A RELATIVE ERROR TOLERANCE. C C X :OUT THE COMPUTED APPROXIMATE SOLUTION OF THE INTEGRAL C EQUATION, EVALUATED AT THE NODE POINTS IN T, IS STORED C IN X ON COMPLETION OF THE ROUTINE. THIS IS TRUE C IRREGARDLESS OF WHETHER OR NOT THE DESIRED ERROR C TOLERANCE WAS ATTAINED. C C T :IN CONTAINS THE NODE POINTS AT WHICH THE SOLUTION OF THE C INTEGRAL EQUATION IS DESIRED. SEE THE VARIABLE NT FOR C MORE INFORMATION. C C NT :IN NT = 0: C T AND X WILL BE SET EQUAL TO THE FINAL GAUSSIAN NODES C AND THE CORRESPONDING SOLUTION VALUES, AND NT WILL BE C SET TO THE NUMBER OF THE SOLUTION VALUES STORED IN X C AND T. THE ARRAYS T AND X SHOULD HAVE DIMENSION AT C LEAST MUPPER, ASSIGNED IN THE CALLING PROGRAM. C NT .GT. 0: C T CONTAINS NT USER SUPPLIED NODES AT WHICH THE SOLUTION C X IS DESIRED. C C NUPPER :IN AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM. N IS C THE ORDER OF A LINEAR SYSTEM WHICH IS BEING USED TO C ITERATIVELY SOLVE A LARGER LINEAR SYSTEM OF ORDER M C WHICH APPROXIMATES THE ABOVE INTEGRAL EQUATION. THE C VALUE OF NUPPER SHOULD BE AN EVEN INTEGER. THIS IS DUE C TO THE GAUSSIAN QUADRATURE PROGRAM BEING USED. C C MUPPER :IN AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM. THE C VALUE OF MUPPER SHOULD BE AN EVEN INTEGER. N AND M ARE C ALWAYS POWERS OF TWO. C C W : OUT 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 C IW :WORK INTEGER TEMPORARY WORKING STORAGE OF LENGTH NUPPER. C C IER :OUT IS AN ERROR COMPLETION CODE WITH THE FOLLOWING POSSIBLE C VALUES: C IER = 0: C THE ROUTINE WAS COMPLETED SATISFACTORILY. EP CONTAINS C THE PREDICTED ERROR. C IER = 1: C THE ERROR TEST WAS NOT SATISFIED. EP CONTAINS THE C PREDICTED ERROR. C IER = 2: C THE ERROR TEST WAS NOT SATISFIED. EP HAS BEEN SET TO C ZERO, BECAUSE THERE IS NOT A SATISFACTORY PREDICTED C ERROR. C IER = 3: C THE ORIGINAL VALUE OF EP WAS TOO SMALL, DUE TO POSSIBLE C ILL-CONDITIONING PROBLEMS IN THE INTEGRAL EQUATION. THE C VALUE OF EP WAS RESET TO A MORE REALISTIC VALUE, AND C THAT TOLERANCE WAS ATTAINED. C IER = 4: C THE ERROR WAS SATISFACTORY AT THE GAUSSIAN NODE POINTS C (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 IS C THE PREDICTED ERROR FOR THE SOLUTION AT THE GAUSSIAN C NODE POINTS OF ORDER MFINAL(INFO(6)). C IER = 5: C THE ANALOGUE OF IER=4, BUT WITH IER=1 AT THE GAUSSIAN C NODE POINTS. C IER = 6: C THE ANALOGUE OF IER=4, BUT WITH IER=3 AT THE GAUSSIAN C NODE POINTS. C IER = 7: C THE MATRIX I-KNN IS SINGULAR. ORDINARILY, THIS CANNOT C OCCUR. C C INFO :OUT THE NUMBERS IN INFO GIVE ADDITIONAL INFORMATION ABOUT C THE FUNCTIONING OF DIEGAU. INFO(1) IS THE ITERATIVE C RATE OF CONVERGENCE IN THE MOST RECENTLY COMPUTED C LINEAR SYSTEM. INFO(2) IS THE RATE OF CONVERGENCE OF C THE GAUSSIAN QUADRATURE VARIANT OF THE NYSTROM METHOD. C INFO(3) IS THE FINAL VALUE OF EP USED AS THE DESIRED C ERROR TOLERANCE. USUALLY INFO(3) WILL EQUAL THE INPUT C VALUE OF EP, UNLESS EP WAS MUCH TOO SMALL. INFO(4) IS C AN APPROXIMATE VALUE FOR THE NORM OF THE INTEGRAL C OPERATOR K, AND IT IS CALCULATED ONLY IF NT .GT. 0. C INFO(5) AND INFO(6) ARE THE FINAL VALUES OF N AND M C USED IN DIEGAU. IF INFO(5)=INFO(6), THEN ITERATION WAS C NOT INVOKED SUCCESSFULLY. 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 GAUSSIAN C QUADRATURE, WITH AN ITERATIVE TECHNIQUE OF SOLUTION FOR THE C RESULTING 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 D1MACH, DAVERR, DGAUS3, DIEGS, DINTRP, DITERT, C DLEAVE, DNORM, DRMIN, LINSYS 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 DIEGAU C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL,INFO DIMENSION X(*),T(*),W(*),IW(*),INFO(6) PARAMETER (P1=.1D0, P5=.5D0) C EXTERNAL KERNEL,RHFCN DATA CUTOFF/P5/, ROOTRT/P1/ C C********************************************************************* C * UNITRD = D1MACH(3) 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 DGAUS3; THEY ARE ACCURATE TO 15 DIGITS. THE * C ROUTINE 'D1MACH' IS THE STANDARD SLATEC AND BELL LABS ROUTINE * C FOR GIVING DOUBLE PRECISION MACHINE CONSTANTS. IT IS ALSO * C USED BELOW IN FUNCTION 'DAVERR'. C * C********************************************************************* C C SET UP THE RELATIVE BASE ADDRESSES FOR THE VARIOUS ARRAYS INTO C WHICH W IS TO BE SPLIT. C C***FIRST EXECUTABLE STATEMENT DIEGAU C 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 I22=I21+NSQ NHALF=N/2 CALL DIEGS(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), * W(I22),IW,INFO) RETURN END SUBROUTINE DIEGS (KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, 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, TMPVEC, IPIVOT, INFO) C***BEGIN PROLOGUE DIEGS C***SUBSIDIARY C***PURPOSE THIS ROUTINE CONTROLS THE SOLUTION OF THE INTEGRAL EQUATION. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DIEGS-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 IFLAG, NT, NUP, MUP, NHALF, IPIVOT(NUP) C DOUBLE PRECISION KERNEL, RHFCN, A, B, EP, X(MUP), T(MUP), CUTOFF, C ROOTRT, UNITRD, 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 SAVE(MUP), XN(NUP), SAVE2(MUP), ASIDE(NUP,4), C ASIDE3(NUP,NUP), IMKNN(NUP,NUP), TMPVEC(NUP) C EXTERNAL KERNEL, RHFCN C PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, P5 = .5D0, C P1 = .1D0, P01 = 0.01D0, P0001 = .0001D0) C C CALL DIEGS (KERNEL, RHFCN, A, B, EP, IFLAG, X, T, NT, NUP, MUP, C IER, CUTOFF, ROOTRT, UNITRD, NHALF, LUFACT, KMM, KMN, C KNM, RHS, R, RH, DELN, TM, TN, XM, XMZ, WM, WN, OLDX, C SAVE, XN, SAVE2, ASIDE, ASIDE3, IMKNN, TMPVEC, IPIVOT, C INFO) 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 DIEGAU. 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 DIEGAU. C C A :IN IS THE LOWER INTEGRATION LIMIT. C C B :IN IS THE UPPER INTEGRATION LIMIT. C C EP :INOUT 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 C IFLAG :IN IFLAG = 0: C EP IS INTRPRETED AS AN ABSOLUTE ERROR TOLERANCE. C IFLAG = 1: C EP IS INTRPRETED AS A RELATIVE ERROR TOLERANCE. C C X :OUT THE COMPUTED APPROXIMATE SOLUTION OF THE INTEGRAL C EQUATION, EVALUATED AT THE NODE POINTS IN T, IS STORED C IN X ON COMPLETION OF THE ROUTINE. THIS IS TRUE C IRREGARDLESS OF WHETHER OR NOT THE DESIRED ERROR C TOLERANCE WAS ATTAINED. C C T :IN CONTAINS THE NODE POINTS AT WHICH THE SOLUTION OF THE C INTEGRAL EQUATION IS DESIRED. SEE THE VARIABLE NT FOR C MORE INFORMATION. C C NT :IN NT = 0: C T AND X WILL BE SET EQUAL TO THE FINAL GAUSSIAN NODES C AND THE CORRESPONDING SOLUTION VALUES, AND NT WILL BE C SET TO THE NUMBER OF THE SOLUTION VALUES STORED IN X C AND T. THE ARRAYS T AND X SHOULD HAVE DIMENSION AT C LEAST MUPPER, ASSIGNED IN THE CALLING PROGRAM. C NT .GT. 0: C T CONTAINS NT USER SUPPLIED NODES AT WHICH THE SOLUTION C X IS DESIRED. C C NUP :IN AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM. N IS C THE ORDER OF A LINEAR SYSTEM WHICH IS BEING USED TO C ITERATIVELY SOLVE A LARGER LINEAR SYSTEM OF ORDER M C WHICH APPROXIMATES THE ABOVE INTEGRAL EQUATION. THE C VALUE OF NUPPER SHOULD BE AN EVEN INTEGER. THIS IS DUE C TO THE GAUSSIAN QUADRATURE PROGRAM BEING USED. C C MUP :IN AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM. THE C VALUE OF MUPPER SHOULD BE AN EVEN INTEGER. N AND M ARE C ALWAYS POWERS OF TWO. C C IER :OUT IS AN ERROR COMPLETION CODE WITH THE FOLLOWING POSSIBLE C VALUES: C IER = 0: C THE ROUTINE WAS COMPLETED SATISFACTORILY. EP CONTAINS C THE PREDICTED ERROR. C IER = 1: C THE ERROR TEST WAS NOT SATISFIED. EP CONTAINS THE C PREDICTED ERROR. C IER = 2: C THE ERROR TEST WAS NOT SATISFIED. EP HAS BEEN SET TO C ZERO, BECAUSE THERE IS NOT A SATISFACTORY PREDICTED C ERROR. C IER = 3: C THE ORIGINAL VALUE OF EP WAS TOO SMALL, DUE TO POSSIBLE C ILL-CONDITIONING PROBLEMS IN THE INTEGRAL EQUATION. THE C VALUE OF EP WAS RESET TO A MORE REALISTIC VALUE, AND C THAT TOLERANCE WAS ATTAINED. C IER = 4: C THE ERROR WAS SATISFACTORY AT THE GAUSSIAN NODE POINTS C (IER=0), BUT THE INTERPOLATION PROCESS(DUE TO C NT .GT. 0) MAY NOT PRESERVE THIS ACCURACY. CHECK THE C VALUE OF DNORM(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 IS C THE PREDICTED ERROR FOR THE SOLUTION AT THE GAUSSIAN C NODE POINTS OF ORDER MFINAL(INFO(6)). C IER = 5: C THE ANALOGUE OF IER=4, BUT WITH IER=1 AT THE GAUSSIAN C NODE POINTS. C IER = 6: C THE ANALOGUE OF IER=4, BUT WITH IER=3 AT THE GAUSSIAN C NODE POINTS. C IER = 7: C THE MATRIX I-KNN IS SINGULAR. ORDINARILY, THIS CANNOT C OCCUR. 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 ROOTRT :IN DEPENDS ON THE NUMERICAL INTEGRATION RULE BEING USED. C SEE THE REFERENCES OF DIEGAU IN ORDER TO SET THEM C PROPERLY WHEN CHANGING NUMERICAL INTEGRATION RULES. C C UNITRD :IN UNITRD IS THE SMALLEST NUMBER U FOR WHICH 1+U .GT. 1 . C UNITRD VARIES WITH THE COMPUTER AND WITH THE ARITHMETIC C BEING USED. TO CHANGE TO ANOTHER ARITHMETIC, UNITRD C MUST BE CHANGED. BUT UNITRD ALSO REFLECTS THE ACCURACY C OF THE CONSTANTS USED IN SUBROUTINE DGAUS3; THEY ARE C ACCURATE TO 15 DIGITS. THE ROUTINE 'D1MACH' IS THE C STANDARD SLATEC AND BELL LABS ROUTINE FOR GIVING DOUBLE C PRECISION MACHINE CONSTANTS. IT IS ALSO USED BELOW IN C FUNCTION 'DAVERR'. C C NHALF :IN NHALF = (NUP + 1)/ 2. C C LUFACT :WORK IS AN NUP BY NUP ARRAY, IS THE SAME AS W(I1) IN DIEGAU. 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 DIEGAU. 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 DIEGAU. 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 DIEGAU. 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 DIEGAU. C C R :WORK IS AN MUP ARRAY, IS THE SAME AS W(I6) IN DIEGAU. C C RH :WORK IS AN NUP ARRAY, IS THE SAME AS W(I7) IN DIEGAU. C C DELN :WORK IS AN NUP ARRAY, IS THE SAME AS W(I8) IN DIEGAU. C C TM :OUT THE M NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B). IT C IS AN MUP ARRAY, IS THE SAME AS W(I9) IN DIEGAU. C C TN :OUT THE N NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B). IT C IS AN NUP ARRAY, IS THE SAME AS W(I10) IN DIEGAU. C C XM :OUT THE NYSTROM INTERPOLATES. IT IS AN MUP ARRAY, IS THE C SAME AS W(I11) IN DIEGAU. 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 DIEGAU. C C WM :OUT THE WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B), WITH M C NODES. IT IS AN MUP ARRAY, IS THE SAME AS W(I13) IN C DIEGAU. C C WN :OUT THE WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B), WITH N C NODES. IT IS AN NUP ARRAY, IS THE SAME AS W(I14) IN C DIEGAU. C C OLDX :INOUT IS AN MUP ARRAY, IS THE SAME AS W(I15) IN DIEGAU. C C SAVE :WORK IS AN MUP ARRAY, IS THE SAME AS W(I16) IN DIEGAU. C C XN :OUT THE NYSTROM INTERPOLATES. IT IS AN NUP ARRAY, IS THE C SAME AS W(I17) IN DIEGAU. C C SAVE2 :WORK IS AN MUP ARRAY, IS THE SAME AS W(I18) IN DIEGAU. C C ASIDE :WORK IS AN NUP BY 4 ARRAY, IS THE SAME AS W(19) IN DIEGAU. C C ASIDE3 :WORK IS AN NUP BY NUP ARRAY, IS THE SAME AS W(I20) IN C DIEGAU. C C IMKNN :OUT THE LINEAR SYSTEM OF THE INTEGRAL EQUATION. IT IS AN C NUP BY NUP ARRAY, IS THE SAME AS W(I21) IN DIEGAU. C C TMPVEC :WORK IS AN NUP ARRAY, IS THE SAME AS W(22) IN DIEGAU. 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 DIEGAU. C C INFO :OUT THE NUMBERS IN INFO GIVE ADDITIONAL INFORMATION ABOUT C THE FUNCTIONING OF DIEGAU. INFO(1) IS THE ITERATIVE C RATE OF CONVERGENCE IN THE MOST RECENTLY COMPUTED C LINEAR SYSTEM. INFO(2) IS THE RATE OF CONVERGENCE OF C THE GAUSSIAN QUADRATURE VARIANT OF THE NYSTROM METHOD. C INFO(3) IS THE FINAL VALUE OF EP USED AS THE DESIRED C ERROR TOLERANCE. USUALLY INFO(3) WILL EQUAL THE INPUT C VALUE OF EP, UNLESS EP WAS MUCH TOO SMALL. INFO(4) IS C AN APPROXIMATE VALUE FOR THE NORM OF THE INTEGRAL C OPERATOR K, AND IT IS CALCULATED ONLY IF NT .GT. 0. C INFO(5) AND INFO(6) ARE THE FINAL VALUES OF N AND M C USED IN DIEGAU. IF INFO(5)=INFO(6), THEN ITERATION WAS C NOT INVOKED SUCCESSFULLY. C C***SEE ALSO DIEGAU C***ROUTINES CALLED DAVERR, DGAUS3, DINTRP, DITERT, DLEAVE, DNORM, C DRMIN, LINSYS 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 DIEGS C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL,LUFACT,KMM,KMN,KNM,IMKNN,DNORM, * NUMR1,NUMR2,INFO INTEGER FLAG,OLDM DIMENSION X(*),T(*),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),IPIVOT(NUP),TMPVEC(NUP),INFO(6) PARAMETER (ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0, P5=.5D0, * P1=.1D0, P01=0.01D0, P0001=.0001D0) C EXTERNAL KERNEL,RHFCN C C***FIRST EXECUTABLE STATEMENT DIEGS C C INITIALIZATION C LOOP=1 N=2 M=2*N EPS=EP C 'EPS' IS THE ERROR TOLERANCE BEING REQUESTED OF THE PROGRAM. C IT MAY BE INCREASED DURING THE COMPUTATION, FOR VARIOUS C POSSIBLE REASONS. 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 ITERATIVE METHOD OF SOLUTION. C C CREATE THE NODES AND WEIGHTS TN(I) AND WN(I), I=1,...,N CALL DGAUS3(A,B,N,TN,WN) C C SET UP MATRIX FOR (I-KN)*XN=RHFCN DO 20 J=1,N DO 10 I=1,N 10 IMKNN(I,J)=-WN(J)*KERNEL(TN(I),TN(J)) XMZ(J)=RHFCN(TN(J)) 20 IMKNN(J,J)=IMKNN(J,J)+ONE C C THIS IS THE ENTRANCE WHEN ITERATION FAILS AND WE NEED TO C INCREASE N TO OBTAIN A BETTER ITERATION RATE. 30 CONTINUE C 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********************************************************************* C * CALL LINSYS(IMKNN,LUFACT,N,XMZ,XN,2,RCOND,IPIVOT,NUP, * TMPVEC,ELINSY) C * C LINSYS IS A GENERAL LINEAR EQUATION SOLVER, BASED ON THE * C LINPACK PROGRAMS DGECO AND DGESL. LINSYS HAS SPECIAL * C OPTIONS WHICH ARE USED IN THIS PROGRAM, AND THUS CARE * C SHOULD BE TAKEN IN REPLACING IT. * C * C LINSYS IS ALSO USED IN THE SUBROUTINE DITERT. * C * C********************************************************************* C IF(RCOND .EQ. ZERO) THEN IER=7 RETURN END IF C C UPDATE CONDITION NUMBER AND OTHER QUANTITIES FOR ERROR ANALYSIS. COND=ONE/RCOND PASTER = D1MACH(3) CALL DAVERR(ELINSY,AVERR,PASTER) RELMIN=DRMIN(N,N,COND,UNITRD,AVERR) C C SET UP APPROXIMATE RATE OF CONVERGENCE OF SOLUTIONS XN TO TRUE C SOLUTION X. ALSO SET UP DESIRED RATIO FOR ITERATIVE METHOD. IF(LOOP .EQ. 1) THEN INFO(2)=P5 R1RAT=ROOTRT GO TO 40 ELSEIF (LOOP .EQ. 2) THEN NUMR2=DNORM(XN,OLDX,N,1) DENR2=ZERO INFO(2)=P5 R1RAT=ROOTRT ELSE C LOOP .GT. 2 NUMR2=DNORM(XN,OLDX,N,1) INFO(2)=MIN(P5,MAX(NUMR2/DENR2,P0001)) R1RAT=MIN(ROOTRT,SQRT(INFO(2))) END IF C C CHECK FOR ERROR IN XN USING TEST INVOLVING INFO(2) AND OLDX, C ACCORDING TO THE THEORY FOR ASYMPTOTIC ERROR BOUNDS. C MODIFY ERROR IF IT IS OUTSIDE PRECISION RANGE OF COMPUTER, C POSSIBLY DUE TO ILL-CONDITIONING. ERROR=(INFO(2)/(ONE-INFO(2)))*NUMR2 XNORM=DNORM(XN,XN,N,0) RELERR=ERROR/XNORM IF(IFLAG .EQ. 0) THEN EPS=MAX(EP,XNORM*RELMIN) ERROR=XNORM*MAX(RELERR,RELMIN) ELSE C IFLAG=1. EPS = MAX(EP,RELMIN) ERROR = MAX(RELERR,RELMIN) END IF C C IF(ERROR .LE. EPS) THEN C EXIT FOR SUCCESSFUL RETURN. ITERATION WAS NOT NECESSARY. CALL DLEAVE(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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN ELSE C THE ERROR IN XN IS NOT SUFFICIENTLY SMALL. INITIALIZE C FOR ANOTHER LOOP ON N (OR M, IF ITERATION IS SUCCESSFUL). DENR2=NUMR2 END IF C C ATTEMPT TO SOLVE (I-KM)*XM=RHFCN ITERATIVELY, CHECKING TO C SEE IF THE RATE OF CONVERGENCE IS SUFFICIENTLY FAST SO AS C TO ENTER STAGE B. C C CALCULATE TM(I) AND WM(I), I=1,...,M. 40 CALL DGAUS3(A,B,M,TM,WM) FLAG=0 C CALCULATE INITIAL GUESS XMZ FOR ITERATION METHOD. CALL DINTRP(TM,XMZ,M,TN,WN,XN,N,KERNEL,RHFCN,RHS) DO 50 I=1,M 50 OLDX(I)=XMZ(I) C C CALCULATE FIRST ITERATE. CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM,RHS, * IMKNN,LUFACT,R,RH,DELN,IPIVOT,TMPVEC,NUP,NHALF,ELINSY,FLAG) CALL DAVERR(ELINSY,AVERR,PASTER) DENR1=DNORM(XM,XMZ,M,1) FLAG=1 DO 60 I=1,M 60 XMZ(I)=XM(I) C CALCULATE SECOND ITERATE. CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM,RHS, * IMKNN,LUFACT,R,RH,DELN,IPIVOT,TMPVEC,NUP,NHALF,ELINSY,FLAG) CALL DAVERR(ELINSY,AVERR,PASTER) NUMR1=DNORM(XM,XMZ,M,1) C C CALCULATE THE RATE OF CONVERGENCE OF THE ITERATION METHOD. INFO(1)=MAX(NUMR1/DENR1,P0001) C C THE FOLLOWING VARIABLE 'RATE' IS USED IN COMPUTING THE SPEED C OF CONVERGENCE INFO(2) OF THE SOLUTIONS XM TO X. IT IS USED TO C CAPTURE THE SPEED OF CONVERGENCE OF THE ITERATION METHOD AS C IT AFFECTS THE ACCURACY OF THE FINAL ITERATE USED IN APPROXI- C MATING THE DESIRED SOLUTION XM. IT IS PROBABLY OF LITTLE C REAL IMPORTANCE, BUT IS AN EXTRA SAFEGUARD. RATE=INFO(1) C C CHECK THE SIZE OF M. IF(M .GT. NUP) THEN C N CANNOT BE INCREASED ANY FURTHER. SEE IF ITERATION WILL C WORK, EVEN IF IT CONVERGES SLOWER THAN DESIRABLE. IF(INFO(1) .GT. CUTOFF) THEN C THE ITERATES ARE CONVERGING VERY SLOWLY OR NOT AT C ALL. THUS RETURN WITHOUT FURTHER ATTEMPTS TO LESSEN C THE ERROR. IF(LOOP .EQ. 1) THEN CALL DLEAVE(2,N,N,XN,TN,WN,ZERO,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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN ELSE CALL DLEAVE(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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN END IF ELSE C THE ITERATES ARE CONVERGING RAPIDLY ENOUGH TO CONTINUE C WITH ITERATION IN STAGE B. OLDM=N ASIDE2(2)=M GO TO 200 END IF END IF 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. IF(INFO(1) .GT. R1RAT) THEN C THE ITERATION DID NOT WORK WELL ENOUGH, AND STAGE A IS C TO BE REPEATED. RE-INITIALIZE FOR SOLVING (I-KN)*XN=RHFCN C AGAIN WITH A LARGER N. N=M LOOP=LOOP+1 M=2*N DO 80 J=1,N DO 70 I=1,N 70 IMKNN(I,J)=-KMM(I,J) WN(J)=WM(J) TN(J)=TM(J) XMZ(J)=RHS(J) 80 IMKNN(J,J)=IMKNN(J,J)+ONE GO TO 30 END IF C 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 CONTROL HAS TO RETURN TO STAGE A. OLDM=M DO 90 I=1,M ASIDE(I,1)=OLDX(I) ASIDE(I,2)=WM(I) ASIDE(I,3)=TM(I) 90 ASIDE(I,4)=RHS(I) ASIDE2(1)=LOOP ASIDE2(2)=M ASIDE2(3)=INFO(2) ASIDE2(4)=DENR2 ASIDE2(5)=R1RAT DO 100 J=1,M DO 100 I=1,M 100 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 C ACCURATE COMPARED TO THE TRUE XM. 200 RATE=INFO(1)*RATE TEMP=DNORM(XM,OLDX,M,1) C TEMP IS AN ESTIMATE OF NUMR2. C TEMP2 IS AN ESTIMATE OF INFO(2). IF(LOOP .EQ. 1) THEN TEMP2=P5 ELSE TEMP2=TEMP/DENR2 END IF C C RT IS AN "ESTIMATE" FOR THE RATE OF CONVERGENCE OF {XM} TO X. C THE USE OF THIS QUANTITY WAS DETEDRMINED EXPERIMENTALLY OVER C A WIDE CLASS OF PROBLEMS. FOR A YET MORE CONSERVATIVE C ITERATION TEST, FORCE RT TO BE LARGER. AT PRESENT, RT .LE. 0.1. RT= MIN(P01,MAX(TEMP2,P0001))/TWO XNORM=DNORM(XM,XM,M,0) C ESTERR IS AN ESTIMATE OF THE ERROR DNORM(X-XM). ESTERR=XNORM*MAX(RELMIN,(RT/(ONE-RT))*TEMP/XNORM) C TEST IS AN ESTIMATE OF THE NEEDED SIZE FOR THE DIFFERENCE C DNORM(XM-XMZ) IN THE MOST CURRENT VALUES OF THE ITERATES C XM AND XMZ. TEST=((ONE-INFO(1))/INFO(1))*ESTERR IF(NUMR1 .LE. TEST) GO TO 260 C C ITERATE NOT SUFFICIENTLY ACCURATE. INITIALIZE FOR C COMPUTATION OF ANOTHER ITERATE. 210 DENR1=NUMR1 DO 220 I=1,M 220 XMZ(I)=XM(I) CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM,RHS, * IMKNN,LUFACT,R,RH,DELN,IPIVOT,TMPVEC,NUP,NHALF,ELINSY,FLAG) CALL DAVERR(ELINSY,AVERR,PASTER) NUMR1=DNORM(XM,XMZ,M,1) C INFO(1)=MAX(NUMR1/DENR1,P0001) IF(INFO(1) .LE. CUTOFF) GO TO 200 C 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. C 230 MNEW=ASIDE2(2) IF(MNEW .GT. NUP) THEN C ABORTIVE EXIT FROM STAGE B. N CANNOT BE INCREASED FURTHER, C AND INFO(1) IS NOT SUFFICIENTLY SMALL. IF(LOOP .EQ. 1) THEN CALL DLEAVE(2,N,N,XN,TN,WN,ZERO,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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN ELSE CALL DGAUS3(A,B,OLDM,TM,WM) CALL DLEAVE(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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN END IF ELSE C INITIALIZE FOR RETURN TO STAGE A. N=MNEW DO 250 J=1,N DO 240 I=1,N 240 IMKNN(I,J)=-ASIDE3(I,J) OLDX(J)=ASIDE(J,1) WN(J)=ASIDE(J,2) TN(J)=ASIDE(J,3) XMZ(J)=ASIDE(J,4) 250 IMKNN(J,J)=IMKNN(J,J)+ONE M=2*N LOOP=ASIDE2(1)+ONE INFO(2)=ASIDE2(3) DENR2=ASIDE2(4) R1RAT=ASIDE2(5) GO TO 30 END IF C C C AN ACCURATE VALUE OF XM HAS BEEN OBTAINED. INFO(2) 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) THEN DENR2 = TEMP LOOP = 2 ELSE NUMR2 = TEMP INFO(2) = MAX(P0001,RATE,MIN(NUMR2/DENR2,P5)) DENR2 = NUMR2 LOOP = LOOP+1 END IF C ERROR=(INFO(2)/(ONE-INFO(2)))*TEMP XNORM=DNORM(XM,XM,M,0) RELERR=ERROR/XNORM RELMIN=DRMIN(N,M,COND,UNITRD,AVERR) IF(IFLAG .EQ. 0) THEN EPS=MAX(EP,XNORM*RELMIN) ERROR=XNORM*MAX(RELERR,RELMIN) ELSE C IFLAG=1. EPS=MAX(EP,RELMIN) ERROR=MAX(RELERR,RELMIN) END IF C IF(ERROR .LE. EPS) THEN CALL DLEAVE(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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN END IF C C THE ERROR IN XM IN NOT SUFFICIENTLY SMALL. INCREASE M, IF C POSSIBLE, AND CALCULATE A NEW MORE ACCURATE XM. MNEW=2*M IF(MNEW .GT. MUP) THEN C M CANNOT BE INCREASED ANY FURTHER. CALL DLEAVE(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,IPIVOT,TMPVEC, * NUP,NHALF,XNORM,INFO) RETURN END IF C C M IS INCREASED AND TWO MORE ITERATES ARE COMPUTED WITH C THE NEW M. OLDM=M M=MNEW DO 270 I=1,OLDM SAVE2(I)=WM(I) 270 SAVE(I)=TM(I) CALL DGAUS3(A,B,M,TM,WM) FLAG=0 CALL DINTRP(TM,XMZ,M,SAVE,SAVE2,XM,OLDM,KERNEL,RHFCN,RHS) DO 280 I=1,OLDM 280 SAVE(I)=XM(I) DO 290 I=1,M 290 OLDX(I)=XMZ(I) C CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM,RHS, * IMKNN,LUFACT,R,RH,DELN,IPIVOT,TMPVEC,NUP,NHALF,ELINSY,FLAG) CALL DAVERR(ELINSY,AVERR,PASTER) DENR1=DNORM(XM,XMZ,M,1) FLAG=1 DO 300 I=1,M 300 XMZ(I)=XM(I) CALL DITERT(KERNEL,RHFCN,N,TN,WN,M,TM,WM,XM,XMZ,KMM,KMN,KNM,RHS, * IMKNN,LUFACT,R,RH,DELN,IPIVOT,TMPVEC,NUP,NHALF,ELINSY,FLAG) CALL DAVERR(ELINSY,AVERR,PASTER) NUMR1=DNORM(XM,XMZ,M,1) C INFO(1)=NUMR1/DENR1 RATE=INFO(1) IF(INFO(1) .LE. CUTOFF) THEN C THE RATE OF CONVERGENCE IS ACCEPTABLE. GO TO 200 ELSE C THE RATE OF CONVERGENCE IS NOT ACCEPTABLE. GO TO 230 END IF END SUBROUTINE DITERT (KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XMZ, * KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, IPIVOT, * TMPVEC, NUP, NHALF, ELINSY, IFLG) C***BEGIN PROLOGUE DITERT C***SUBSIDIARY C***PURPOSE THIS ROUTINE CALCULATES ONE ITERATE XM GIVEN THE INITIAL C GUESS XM_Z. 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, IPIVOT(NUP), NUP, NHALF, 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), IMKNN(NUP,NUP), C LUFACT(NUP,NUP), R(M), RH(N), DELN(N), C TMPVEC(NUP), ELINSY C PARAMETER (ZERO = 0.0D0) C C CALL DITERT (KERNEL, RHFCN, N, TN, WN, M, TM, WM, XM, XM_Z, KMM, C KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, IPIVOT, C TMPVEC, NUP, NHALF, ELINSY, 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 N NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B). C C WN :IN THE N WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B). C C M :IN THE NUMBER OF SUBINTERVALS ON (A, B), USUALLY, M IS C GREATER THEN N. C C TM :IN THE M NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B). C C WM :IN THE M WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B). 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 AT C (TM(I),TN(J)). 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 IMKNN :INOUT THE LINEAR SYSTEM OF THE INTEGRAL EQUATION. C C LUFACT :OUT THIS CONTAINS (OR WILL CONTAIN) THE LU FACTORIZATION C OF THE MATRIX IN IMKNN. ALSO SEE LINSYS. 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 = NUP. C C TMPVEC :WORK IS AN NUP ARRAY. ALSO SEE LINSYS. 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 ELINSY :OUT IS THE ERRMAX IN LINSYS. C C IFLG :IN IFLG = 0: C THE MATRICES KMM AND KNM MUST BE COMPUTED AND STORED. C C***ROUTINES CALL LINSYS 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,IMKNN,LUFACT 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),R(M), * LUFACT(NUP,NUP),RH(N),DELN(N),IPIVOT(NUP),TMPVEC(NUP) PARAMETER (ZERO=0.0D0) 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 100 C C IF(IFLG .EQ. 0) THEN C THE MATRICES KMM, KMN, AND KNM MUST BE COMPUTED C 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*XMZ)(TM(I)) DO 40 I=1,M SUM=ZERO DO 30 J=1,M 30 SUM=SUM+KMM(I,J)*XMZ(J) 40 R(I)=RHS(I)-(XMZ(I)-SUM) C C COMPUTE RH=KM*R AT ALL TN(I). DO 50 I=1,N RH(I)=ZERO DO 50 J=1,M 50 RH(I)=RH(I)+KNM(I,J)*R(J) C CALCULATE DELN=((I-KN)**(-1))*KM*R AT ALL TN(I). C C********************************************************************* C * CALL LINSYS(IMKNN,LUFACT,N,RH,DELN,4,RCOND,IPIVOT,NUP, * TMPVEC,ELINSY) C * C SEE THE ORIGINAL REFERENCE IN DIEGS. * C********************************************************************* C C CALCULATE NEW XM. DO 80 I=1,M SUM=ZERO DO 60 J=1,M 60 SUM=SUM+KMM(I,J)*R(J) DO 70 J=1,N 70 SUM=SUM+KMN(I,J)*DELN(J) 80 XM(I)=SUM+R(I)+XMZ(I) RETURN C C ENTRANCE WHEN M .GT. NUP. C CALCULATE RESIDUALS. 100 DO 120 I=1,M SUM=ZERO DO 110 J=1,M 110 SUM=SUM+WM(J)*KERNEL(TM(I),TM(J))*XMZ(J) 120 R(I)=RHS(I)-(XMZ(I)-SUM) C CALCULATE RH=KM*R. DO 130 I=1,N RH(I)=ZERO DO 130 J=1,M 130 RH(I)=RH(I)+WM(J)*KERNEL(TN(I),TM(J))*R(J) C C********************************************************************* C * CALL LINSYS(IMKNN,LUFACT,N,RH,DELN,4,RCOND,IPIVOT,NUP, * TMPVEC,ELINSY) C * C SEE THE ORIGINAL REFERENCE IN DIEGS. * C********************************************************************* C C CALCULATE XM. DO 160 I=1,M SUM=ZERO DO 140 J=1,M 140 SUM=SUM+WM(J)*KERNEL(TM(I),TM(J))*R(J) DO 150 J=1,N 150 SUM=SUM+WN(J)*KERNEL(TM(I),TN(J))*DELN(J) 160 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 TM(M), XM(M), TN(N), WN(N), XN(N), KERNEL, RHFCN, C 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 GAUSSIAN QUADRATURE ON (A, B), WITH C M 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 GAUSSIAN QUADRATURE ON (A, B), WITH C N NODES. C C WN :IN THE WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B), WITH N C N 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 DIEGAU. 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 DIEGAU. 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 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 SAVE RHS(I) FOR LATER USE IN ITERT. C CALCULATE NYSTROM INTERPOLATING FORMULA. C C***FIRST EXECUTABLE STATEMENT DINTRP C DO 10 I=1,M RHS(I)=RHFCN(TM(I)) XM(I)=RHS(I) DO 10 J=1,N 10 XM(I)=XM(I)+WN(J)*KERNEL(TM(I),TN(J))*XN(J) RETURN END C SUBROUTINE DLEAVE (IERSET, NF, MF, XV, TV, WV, 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, * IPIVOT, TMPVEC, NUP, NHALF, XNORM, INFO) C***BEGIN PROLOGUE DELEAVE C***SUBSIDIARY C***PURPOSE THIS ROUTINE SETS ALL NECESSARY PARAMETERS FOR LEAVING C DIEGAU. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DLEAVE-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 IERSET, NF, MF, IFLAG, NT, IER, IPIVOT, NUP, NHALF C DOUBLE PRECISION XV(MF), TV(MF), WV(MF), ERROR, KERNEL, RHFCN, EP, C X(1), T(1), EPS, ELINSY, TN(NF), WN(NF), TM(MF), C WM(MF), XM(MF), XMZ(MF), KMM(NUP,NUP), C KMN(NUP,NHALF), KNM(NHALF,NUP), RHS(MF), C IMKNN(NUP,NUP), LUFACT(NUP,NUP), R(MF), RH(NF), C DELN(NF), TMPVEC(NUP), XNORM, INFO(6) C EXTERNAL KERNEL, RHFCN C PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C C CALL DLEAVE (IERSET, NF, MF, XV, TV, WV, ERROR, KERNEL, RHFCN, EP, C IFLAG, X, T, NT, IER, EPS, ELINSY, TN, WN, TM, WM, XM, C XMZ, KMM, KMN, KNM, RHS, IMKNN, LUFACT, R, RH, DELN, C IPIVOT, TMPVEC, NUP, NHALF, XNORM, INFO) C C *ARGUMENTS: C C IERSET :IN SEE DIEGS. C C NF :IN THE NUMBER OF NODE POINTS, SEE DIEGIS. C C MF :IN THE NUMBER OF NODE POINTS, SEE DIEGIS. C C XV :IN THE NYSTROM INTERPOLATES. IT IS AN MF ARRAY. C C TV :IN THE NF NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B). IT C IS AN MF ARRAY. C C WV :IN THE WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B), WITH MF C NODES. IT IS AN MF ARRAY. C C ERROR :IN SEE DIEGS. C C KERNEL :EXT THIS IS A DOUBLE PRECISION FUNCTION OF TWO VARIABLES. IT C MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM C CALLING DIEGAU. C C RHFCN :EXT THIS IS A DOUBLE PRECISION FUNCTION OF ONE VARIABLE. IT C MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM C CALLING DIEGAU. C C EP :INOUT THE DESIRED ERROR. SEE THE DISCUSSION OF EP IN DIEGS FOR C MORE INFORMATION. C C IFLAG :IN IFLAG = 0: C EP IS INTRPRETED AS AN ABSOLUTE ERROR TOLERANCE. C IFLAG = 1: C EP IS INTRPRETED AS A RELATIVE ERROR TOLERANCE. C C X :OUT NT = 0: C NO NYSTROM INTERPOLATION IS DESIRED. RETURN THE VALUES C AT THE GAUSSIAN NODE POINTS(X(I) = XV(I)). C NT .GT. 0: C COMPUTE NYSTROM INTERPOLATES AT THE NODES IN T. C C T :INOUT NT = 0: C NO NYSTROM INTERPOLATION IS DESIRED. RETURN THE VALUES C AT THE GAUSSIAN NODE POINTS(T(I) = TV(I)). C NT .GT. 0: C T CONTAINS NT USER SUPPLIED NODES AT WHICH THE SOLUTION C X IS DESIRED. C C NT :IN NT = 0: C T AND X WILL BE SET EQUAL TO TV AND XV RESPECTIVELY. C NT .GT. 0: C T CONTAINS NT USER SUPPLIED NODES AT WHICH THE SOLUTION C X IS DESIRED. C C IER :OUT IS AN ERROR COMPLETION CODE WITH THE FOLLOWING POSSIBLE C VALUES: C IER = 0: C THE ROUTINE WAS COMPLETED SATISFACTORILY. EP CONTAINS C THE PREDICTED ERROR. C IER = 1: C THE ERROR TEST WAS NOT SATISFIED. EP CONTAINS THE C PREDICTED ERROR. C IER = 2: C THE ERROR TEST WAS NOT SATISFIED. EP HAS BEEN SET TO C ZERO, BECAUSE THERE IS NOT A SATISFACTORY PREDICTED C ERROR. C IER = 3: C THE ORIGINAL VALUE OF EP WAS TOO SMALL, DUE TO POSSIBLE C ILL-CONDITIONING PROBLEMS IN THE INTEGRAL EQUATION. THE C VALUE OF EP WAS RESET TO A MORE REALISTIC VALUE, AND C THAT TOLERANCE WAS ATTAINED. C IER = 4: C THE ERROR WAS SATISFACTORY AT THE GAUSSIAN NODE POINTS C (IER=0), BUT THE INTERPOLATION PROCESS(DUE TO C NT .GT. 0) MAY NOT PRESERVE THIS ACCURACY. CHECK THE C VALUE OF DNORM(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 IS C THE PREDICTED ERROR FOR THE SOLUTION AT THE GAUSSIAN C NODE POINTS OF ORDER INFO(6). C IER = 5: C THE ANALOGUE OF IER=4, BUT WITH IER=1 AT THE GAUSSIAN C NODE POINTS. C IER = 6: C THE ANALOGUE OF IER=4, BUT WITH IER=3 AT THE GAUSSIAN C NODE POINTS. C IER = 7: C THE MATRIX I-KNN IS SINGULAR. ORDINARILY, THIS CANNOT C OCCUR. C EPS :IN IS THE ERROR TOLERANCE BEING REQUESTED OF THE PROGRAM C DIEGS. C C ELINSY :IN IS THE ERRMAX FROM LINSYS. C C TN :IN THE NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B), C WITH NF NODES. C C C WN :IN THE WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B), C WITH NF NODES. C C TM :IN THE NODE POINTS FOR GAUSSIAN QUADRATURE ON (A, B), C WITH MF NODES. C C C WM :IN THE WEIGHTS FOR GAUSSIAN QUADRATURE ON (A, B), C WITH MF NODES. C C XM :OUT THE NYSTROM INTERPOLATES. IT IS AN MF ARRAY. C C XMZ :OUT THE NYSTROM INTERPOLATES. IT IS AN MF ARRAY. 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 IMKNN :IN THE LINEAR SYSTEM OF THE INTEGRAL EQUATION. C C LUFACT :OUT THIS CONTAINS (OR WILL CONTAIN) THE LU FACTORIZATION C OF THE MATRIX IN IMKNN. ALSO SEE LINSYS. 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 = NF. C C IPIVOT :WORK THIS IS TEMPORARY WORKING STORAGE FOR INTEGER C VARIABLES. IT SHOULD HAVE DIMENSION = NUP. C C TMPVEC :WORK IS AN NUP ARRAY. ALSO SEE LINSYS. C C NUP :IN AN UPPER LIMIT ON THE VARIABLE N IN THE PROGRAM( C DIEGAU). N IS THE ORDER OF APPROXIMATE INVERSE WHICH C IS BEING USED TO ITERATIVELY SOLVE A LARGER SYSTEM OF C ORDER M. C C NHALF :IN NHALF = (NUP + 1)/2. C C XNORM :IN THE MAXIMUM NORM OF THE VECTOR XN. C C INFO :OUT THE NUMBERS IN INFO GIVE ADDITIONAL INFORMATION ABOUT C THE FUNCTIONING OF DIEGAU. INFO(1) IS THE ITERATIVE C RATE OF CONVERGENCE IN THE MOST RECENTLY COMPUTED C LINEAR SYSTEM. INFO(2) IS THE RATE OF CONVERGENCE OF C THE GAUSSIAN QUADRATURE VARIANT OF THE NYSTROM METHOD. C INFO(3) IS THE FINAL VALUE OF EP USED AS THE DESIRED C ERROR TOLERANCE. USUALLY INFO(3) WILL EQUAL THE INPUT C VALUE OF EP, UNLESS EP WAS MUCH TOO SMALL. INFO(4) IS C AN APPROXIMATE VALUE FOR THE NORM OF THE INTEGRAL C OPERATOR K, AND IT IS CALCULATED ONLY IF NT .GT. 0. C INFO(5) AND INFO(6) ARE THE FINAL VALUES OF N AND M C USED IN DIEGAU. IF INFO(5)=INFO(6), THEN ITERATION WAS C NOT INVOKED SUCCESSFULLY. C C *DESCRIPTION: C THIS ROUTINE SETS ALL NECESSARY PARAMETERS FOR LEAVING DIEGAU. IF C NT .GT. 0, IT ALSO PERFORMS THE NECESSARY NYSTROM INTERPOLATION AT C THE NODES GIVEN IN T. C C***ROUTINES CALLED DAVERR, DITERT, DNORM, LINSYS 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 DLEAVE C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KERNEL,KMM,KMN,KNM,IMKNN,LUFACT,NUMR1,NORM,XNORM * ,INFO DIMENSION X(1),T(1),XV(MF),TV(MF),WV(MF),TN(NF),WN(NF),WM(MF), * TM(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),DELN(NF),IPIVOT(NUP),TMPVEC(NUP),INFO(6) PARAMETER (ZERO=0.0D0, ONE=1.0D0) EXTERNAL KERNEL,RHFCN C C***FIRST EXECUTABLE STATEMENT DLEAVE C C SET ERROR PARAMETERS FOR RETURN. INFO(5)=NF INFO(6)=MF INFO(3)=EPS C CALCULATE NORM(K) DO 10 I=1,NF 10 IMKNN(I,I)=IMKNN(I,I)-ONE INFO(4)=ZERO DO 30 I=1,NF SUM=ZERO DO 20 J=1,NF 20 SUM=SUM+ABS(IMKNN(I,J)) 30 INFO(4)=MAX(INFO(4),SUM) DO 40 I=1,NF 40 IMKNN(I,I)=IMKNN(I,I)+ONE C IF((EPS .GT. EP) .AND. (ERROR .LE. EPS)) THEN IER=3 C SINCE EPS IS THE SMALLEST ERROR POSSIBLE, SET EP=EPS C FOR THE FINAL ERROR ESTIMATE. EP=EPS ELSE IER=IERSET EP=ERROR END IF C IF(NT .EQ. 0) THEN C NO NYSTROM INTERPOLATION IS DESIRED. RETURN THE C VALUES AT THE GAUSSIAN NODE POINTS. DO 50 I=1,MF X(I)=XV(I) 50 T(I)=TV(I) NT=MF RETURN END IF C SAVEP=EP IF(NF .LT. MF) THEN C ITERATE TO DECREASE THE NOISE LEVEL IN X. THIS SHOULD C REDUCE POSSIBLE ERRORS IN NYSTROM INTERPOLATION. C THE FOLLOWING IS THE ITERATION DIFFERENCE BEING C SOUGHT IN THE ITERATION. IF(IFLAG .EQ. 0) THEN DERROR=((ONE-INFO(1))/INFO(1))*EPS/INFO(4) ELSE DERROR=((ONE-INFO(1))/INFO(1))*EPS*XNORM/INFO(4) END IF C ITLOOP=0 PASTER=ELINSY DO 60 I=1,MF 60 XM(I)=XV(I) C 70 DO 80 I=1,MF 80 XMZ(I)=XM(I) CALL DITERT(KERNEL,RHFCN,NF,TN,WN,MF,TM,WM,XM,XMZ,KMM,KMN, * KNM,RHS,IMKNN,LUFACT,R,RH,DELN,IPIVOT,TMPVEC,NUP, * NHALF,ELINSY,1) C CALL DAVERR(ELINSY,AVERR,PASTER) NUMR1=DNORM(XM,XMZ,MF,1) ITLOOP=ITLOOP+1 IF((NUMR1 .GT. DERROR) .AND. (ITLOOP .LT. 5)) GO TO 70 C C THE ITERATION IS COMPLETE. DO 90 I=1,MF 90 XV(I)=XM(I) C ESTIMATE NEW ERROR BOUND FOR NYSTROM INTERPOLATES. IF(IFLAG .EQ. 0) THEN TEMP=INFO(4)*(INFO(1)/(ONE-INFO(1)))*NUMR1 ELSE TEMP=INFO(4)*(INFO(1)/(ONE-INFO(1)))*NUMR1/XNORM END IF EP=MAX(EP,TEMP) ELSE C NO ITERATION USED IN COMPUTING X. JUST COMPUTE ERROR C ESTIMATE IN NYSTROM INTERPOLATE. IF(IFLAG .EQ. 0) THEN TEMP=INFO(4)*ELINSY*XNORM ELSE TEMP=NORK*ELINSY END IF IF(IER .NE. 2) EP=MAX(EP,TEMP) END IF C C COMPUTE NYSTROM INTERPOLATES AT THE NODES IN T. DO 110 I=1,NT SUM=ZERO DO 100 J=1,MF 100 SUM=SUM+WV(J)*KERNEL(T(I),TV(J))*XV(J) 110 X(I)=RHFCN(T(I))+SUM C C SET FINAL ERROR INDICATORS. 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 DGAUS3 (A, B, N, TV, WV) C***BEGIN PROLOGUE DGAUS3 C***SUBSIDIARY C***PURPOSE THIS ROUTINE CALCULATES THE WEIGHTS AND NODE POINTS FOR C GAUSSIAN QUADRATURE ON (A,B). C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DGAUS3-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 TV(N), WV(N) C PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) C C CALL DGAUS3 (A, B, N, TV, WV) C C *ARGUMENTS: C C A :IN IS THE LOWER INTEGRATION LIMIT. C C B :IN IS THE UPPER INTEGRATION LIMIT. C C N :IN IS THE NUMBER OF NODES. C C TV :OUT INTEGRATION NODES. C C WV :OUT INTEGRATION WEIGHTS. C C *DESCRIPTION: C INTEGRATION WEIGHTS AND NODES ARE TO BE CALCULATED AND STORED IN WV C AND TV, RESPECTIVELY. N IS ASSUMED TO BE A POWER OF TWO. IF C 1 .LE. N .LE. 256, THEN GAUSSIAN QUADRATURE IS USED. IF N .GT. 256, C THEN THE INTERVAL (A,B) IS DIVIDED N/256 TIMES AND THE 256 POINT C FORMULA IS APPLIED TO EACH SUBINTERVAL. 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 DGAUS3 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WV(N),TV(N) DIMENSION W(255),T(255) PARAMETER (ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0) 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/ C C***FIRST EXECUTABLE STATEMENT DGAUS3 C IF(N .EQ. 1) THEN TV(1) = (B+A)/TWO WV(1) = B-A RETURN END IF C LOOP = MAX(1,N/256) FLOOP = LOOP H = (B-A)/FLOOP SCALE = H/TWO M = MIN0(128,N/2) MT = 2*M NPLACE = M-1 DO 10 L=1,LOOP FL = L AL = A+(FL-ONE)*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*(ONE+S)+(ONE-S)*BL)/TWO TV(I2) = (AL*(ONE-S)+(ONE+S)*BL)/TWO WV(I1) = R 10 WV(I2) = R RETURN END C SUBROUTINE DAVERR (ELINSY, AVERR, PASTER) C***BEGIN PROLOGUE DAVERR C***SUBSIDIARY C***PURPOSE UPDATING THE VALUE OF THE CONDITION NUMBER IN DIEGS. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DAVERR-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 DOUBLE PRECISION ELINSY, AVERR, PASTER C PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C C CALL DAVERR (ELINSY, AVERR, PASTER) C C *ARGUMENTS: C C ELINSY :IN IS THE ERRMAX FROM LINSYS. C C AVERR :OUT IS THE CONDITION NUMBER IN DIEGS. C C PASTER :INOUT IS USED IN DAVERR ONLY. C C***ROUTINES CALLED D1MACH 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 DAVERR C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (ZERO=0.0D0, ONE=1.0D0) C C***FIRST EXECUTABLE STATEMENT DAVERR C IF (ELINSY .NE. ZERO) THEN AVERR = SQRT(ELINSY*PASTER) PASTER = ELINSY ELSE AVERR = D1MACH(3) PASTER = AVERR END IF 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 DOUBLE PRECISION FUNCTION DRMIN (N, M, COND, UNITRD, AVERR) C***BEGIN PROLOGUE DRMIN C***SUBSIDIARY C***PURPOSE FIND THE VALUE, RELMIN, USED IN IEGS FOR A LINEAR SYSTEM C (I-KMM)*XM=RHFCN OF ORDER M. C***LIBRARY SLATEC C***CATEGORY I3 C***TYPE DOUBLE PRECISION (DRMIN-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 C DOUBLE PRECISION COND, UNTIRD,AVERR C PARAMETER (R1P5 = 1.5D0) C C VALUE = DRMIN(N, M, COND, UNITRD, AVERR) C C *ARGUMENTS: C C N :IN IS THE ORDER OF THE LU DECOMPOSITION CURRENTLY BEING USED C IN SOLVING THE ORDER M SYSTEM. ALWAYS, M .GE. N; AND C M = N MEANS THAT ITERATION IS NOT BEING USED. C C M :IN THE ORDER M SYSTEM IS BEING SOLVED. C C COND :IN COMES FROM THE SUBROUTINE LINSYS (IMPLICITLY). C C UNTIRD :IN IS DEFINED IN DIEGAU. C C AVERR :IN COMES FROM THE SUBROUTINE DAVERR. C C *FUNCTION RETURN VALUES: C C VALUE : RELMIN. 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 DRMIN C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (R1P5=1.5D0) C C***FIRST EXECUTABLE STATEMENT DRMIN C FLOAT1=M FLOAT2=M/N DRMIN=MAX((FLOAT1**R1P5)*COND*UNITRD, * (FLOAT2**R1P5)*AVERR) 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