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

.