      SUBROUTINE IESLV (KERNEL,RHFCN,A,B,EP,IFLAG,T,X,NT,NUPPER,
     *                     MUPPER,NF,MF,NORM,W,IER)
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
C   GAUSSIAN QUADRATURE, WITH AN ITERATIVE TECHNIQUE OF SOLUTION
C   FOR THE RESULTING LINEAR SYSTEM.
C
C   KERNEL     THESE ARE REAL FUNCTIONS OF TWO AND ONE
C   RHFCN      VARIABLES, RESPECTIVELY. THEY MUST BE DECLARED IN AN
C              EXTERNAL STATEMENT IN THE PROGRAM CALLING IESLV.
C   EP         THE DESIRED ERROR. THE VARIABLE EP IS CHANGED ON
C              COMPLETION OF THE PROGRAM. SEE THE DISCUSSION OF IER
C              AND IFLAG FOR MORE INFORMATION.
C   IFLAG =0   EP IS INTERPRETED AS AN ABSOLUTE ERROR TOLERANCE.
C         =1   EP IS INTERPRETED AS A RELATIVE ERROR TOLERANCE.
C   T          CONTAINS THE NODE POINTS AT WHICH THE SOLUTION OF THE
C              INTEGRAL EQUATION IS DESIRED. SEE THE VARIABLE NT FOR
C              MORE INFORMATION.
C   X          THE COMPUTED APPROXIMATE SOLUTION OF THE INTEGRAL
C              EQUATION, EVALUATED AT THE NODE POINTS IN T, IS
C              STORED IN X ON COMPLETION OF THE ROUTINE. THIS IS
C              TRUE IRREGARDLESS OF WHETHER OR NOT THE DESIRED ERROR
C              TOLERANCE WAS ATTAINED.
C   NT         IF NT=0, THEN T AND X WILL BE SET EQUAL TO THE FINAL
C              GAUSSIAN NODES AND THE CORRESPONDING SOLUTION VALUES,
C              AND NT WILL BE SET TO THE NUMBER OF THE SOLUTION
C              VALUES STORED IN X AND T. THE ARRAYS T AND X SHOULD
C              HAVE DIMENSION AT LEAST MUPPER, ASSIGNED IN THE
C              CALLING PROGRAM.
C              IF NT .GT. 0, THEN T CONTAINS NT USER SUPPLIED NODES
C              AT WHICH THE SOLUTION X IS DESIRED.
C   NUPPER     AN UPPER LIMIT ON THE VARIABLE N IN THIS PROGRAM.
C              N IS THE ORDER OF A LINEAR SYSTEM WHICH IS BEING
C              USED TO ITERATIVELY SOLVE A LARGER LINEAR SYSTEM OF
C              ORDER M WHICH APPROXIMATES THE INTEGRAL EQUATION.
C              (FOR FURTHER DETAILS CONCERNING THE MAXIMUM VALUE
C              THAT N CAN TAKE, SEE THE DESCRIPTION OF NMAX BELOW.)
C   MUPPER     AN UPPER LIMIT ON THE VARIABLE M IN THE PROGRAM.
C              N AND M ARE ALWAYS POWERS OF TWO.
C   NF         SAME AS NFINAL (SEE BELOW)
C   MF         SAME AS MFINAL (SEE BELOW)
C   NORM       SAME AS NORMK (SEE BELOW)
C   W          TEMPORARY WORKING STORAGE FOR THE PROGRAM. IT MUST
C              CONTAIN AT LEAST 5*NU*NU+9*(NU+MU) POSITIONS, WITH
C              NU=NUPPER, MU=MUPPER.
C   IER =0     THIS ERROR COMPLETION CODE MEANS THE ROUTINE WAS
C              COMPLETED SATISFACTORILY. EP CONTAINS THE PREDICTED
C              ERROR.
C       =1     THE ERROR TEST WAS NOT SATISFIED. EP CONTAINS THE
C              PREDICTED ERROR.
C       =2     THE ERROR TEST WAS NOT SATISFIED. EP HAS BEEN SET
C              TO ZERO.
C       =3     THE ORIGINAL VALUE OF EP WAS TOO SMALL, DUE TO
C              POSSIBLE ILL-CONDITIONING PROBLEMS IN THE INTEGRAL
C              EQUATION. THE VALUE OF EP WAS RESET TO A MORE
C              REALISTIC VALUE, AND THAT TOLERANCE WAS ATTAINED.
C       =4     THE ERROR WAS SATISFACTORY AT THE GAUSSIAN NODE
C              POINTS (IER=0), BUT THE INTERPOLATION PROCESS(DUE TO
C              NT .GT. 0) MAY NOT PRESERVE THIS ACCURACY. CHECK THE
C              VALUE OF NORM(K) FOR POSSIBLE INDICATIONS THAT THE
C              INTEGRAL EQUATION MAY BE ALMOST FIRST KIND. SUCH
C              EQUATIONS ARE QUITE ILL-CONDITIONED. THE ERROR IN EP
C              IS THE PREDICTED ERROR FOR THE SOLUTION AT THE
C              GAUSSIAN NODE POINTS OF ORDER MFINAL.
C       =5     THE ANALOGUE OF IER=4, BUT WITH IER=1 AT THE
C              GAUSSIAN NODE POINTS.
C       =6     THE ANALOGUE OF IER=4, BUT WITH IER=3 AT THE
C              GAUSSIAN NODE POINTS.
C
C
C   *** REFERENCES ***
C   (1) AN AUTOMATIC PROGRAM FOR FREDHOLM INTEGRAL EQUATIONS OF THE
C       SECOND KIND, ACM TRANS. MATH SOFTWARE 2(1976), PP.154-171.
C   (2) A SURVEY OF NUMERICAL METHODS FOR THE SOLUTION OF FREDHOLM
C       INTEGRAL EQUATIONS OF THE SECOND KIND, SIAM PUB., 1976,
C       PART II, CHAP. 5.
C
      REAL KERNEL,NORM,NORMK
      DIMENSION X(1),T(1),W(1)
      EXTERNAL KERNEL,RHFCN
C
C*******************************************************************
C                                                                  *
      COMMON/XXINFO/R1,R2,FINLEP,NORMK,NFINAL,MFINAL
C                                                                  *
C   THE NUMBERS IN XXINFO GIVE ADDITIONAL INFORMATION ABOUT THE    *
C   FUNCTIONING OF IESLV. R1 IS THE ITERATIVE RATE OF CONVERGENCE  *
C   IN THE MOST RECENTLY COMPUTED LINEAR SYSTEM. R2 IS THE RATE OF *
C   CONVERGENCE OF THE GAUSSIAN QUADRATURE VARIANT OF THE NYSTROM  *
C   METHOD. FINLEP IS THE FINAL VALUE OF EP USED AS THE DESIRED    *
C   ERROR TOLERANCE. USUALLY FINLEP WILL EQUAL THE INPUT VALUE OF  *
C   EP, UNLESS EP WAS MUCH TOO SMALL. NORMK IS AN APPROXIMATE      *
C   VALUE FOR THE NORM OF THE INTEGRAL OPERATOR K, AND IT IS       *
C   CALCULATED ONLY IF NT .GT. 0.                                  *
C   NFINAL AND MFINAL ARE THE FINAL VALUES OF N AND M USED IN      *
C   IESLV.  IF NFINAL=MFINAL, THEN ITERATION WAS NOT INVOKED       *
C   SUCCESSFULLY.                                                  *
C                                                                  *
C*******************************************************************
C                                                                  *
C   NMAX IS THE MAXIMUM VALUE FOR N THAT IS PERMITTED BY IESLV.    *
C   THUS MIN0(NUPPER,NMAX) IS THE MAXIMUM VALUE FOR N THAT CAN BE  *
C   USED. THERE IS ALSO AN UPPER LIMIT OF 128 ON N IMPOSED BY THE  *
C   SUBROUTINE LNSYS.                                              *
C                                                                  *
              DATA NMAX /64/
C                                                                  *
C*******************************************************************
C                                                                  *
C     UNITRD IS A MACHINE DEPENDENT PARAMETER. ASSIGN UNITRD THE   *
C     VALUE U WHERE U IS THE SMALLEST FLOATING POINT NUMBER SUCH   *
C     THAT 1.0 + U .GT. 1.0.                                       *
C                                                                  *
              UNITRD = SPMPAR(1)
C                                                                  *
C*******************************************************************
C
      CUTOFF = 0.5
      ROOTRT = 0.1
      NUP = MIN0(NUPPER,NMAX)
C
C   SET UP THE RELATIVE BASE ADDRESSES FOR THE VARIOUS ARRAYS INTO
C   WHICH W IS TO BE SPLIT.
      N=NUP
      M=MUPPER
      NSQ=N*N
      I1=1
      I2=I1+NSQ
      I3=I2+NSQ
      I4=I3+NSQ/2
      I5=I4+NSQ/2
      I6=I5+M
      I7=I6+M
      I8=I7+N
      I9=I8+N
      I10=I9+M
      I11=I10+N
      I12=I11+M
      I13=I12+M
      I14=I13+M
      I15=I14+N
      I16=I15+M
      I17=I16+M
      I18=I17+N
      I19=I18+M
      I20=I19+4*N
      I21=I20+NSQ
      NHALF=N/2
      CALL IEGS(KERNEL,RHFCN,A,B,EP,IFLAG,X,T,NT,NUP,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))
      NORM=NORMK
      NF=NFINAL
      MF=MFINAL
      RETURN
      END
