
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


