
C     TITLE: SOLVE AN INTEGRAL EQUATION USING THE NYSTROM METHOD.
C            THE PROGRAM ILLUSTRATES THE RESIDUAL CORRECTION
C            METHOD FOR SOLVING A LINEAR SYSTEM.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LAMBDA,KERMAT
      PARAMETER (MAXN=50)
      DIMENSION A(MAXN,MAXN), B(MAXN), WORK(MAXN), T(MAXN), IPVT(MAXN),
     *     KERMAT(MAXN,MAXN), RESID(MAXN), X0(MAXN)
      COMMON/BLOCKA/A
      COMMON/BLOCKK/KERMAT
      DATA IWRITE/8/
C
      OPEN(IWRITE,FILE='iteration.output')
C     INPUT PROBLEM PARAMETERS.
10    PRINT *, 'GIVE LAMBDA AND N.'
      READ *, LAMBDA,N
      IF (N .EQ. 0)  THEN
          CLOSE(5)
          STOP
      END IF
C
      WRITE(IWRITE,*) 'LAMBDA, N=',LAMBDA,N
C
C     CREATE THE LINEAR SYSTEM.
      DO 20 I=1,N
20      T(I)=(I-.5D0)/N
      DO 40 I=1,N
        DO 30 J=1,N
          KERMAT(I,J)=FCNK(T(I)-T(J))/N
30        A(I,J)=-KERMAT(I,J)
        B(I)=LAMBDA - (ATAN(1.0D0-T(I)) + ATAN(T(I)))
        X0(I)=B(I)
40      A(I,I)=A(I,I)+LAMBDA
C
C     SOLVE THE LINEAR SYSTEM.
      CALL DGECO(A,MAXN,N,IPVT,RCOND,WORK)
      WRITE(IWRITE,1003) RCOND
1003  FORMAT(/,' RCOND=',1PD8.2,/)
      CALL DGESL(A,MAXN,N,IPVT,X0,0)
C
C     CHECK ON THE ERROR WHEN COMPARED TO THE TRUE ANSWER.
90    ERRMAX=0.0D0
      DO 100 I=1,N
        TRUE=1.0D0
        ERROR=TRUE-X0(I)
        ERRMAX=MAX(ERRMAX,ABS(ERROR))
100     CONTINUE
      WRITE(IWRITE,1002) ERRMAX
1002  FORMAT(/,' MAX ABS[F(T(I))-X(I)] = ',1PD9.3,/)
C
C     BEGIN ITERATION SECTION OF PROGRAM.
110   PRINT *, 'GIVE LAMBDA.'
      READ *, LAMBDA
      IF (LAMBDA .EQ. 0.0D0) GO TO 10
      PRINT *, 'GIVE ITERATION ERROR TOLERANCE'
      READ *, EPSLON
      WRITE(IWRITE,1004) LAMBDA,EPSLON
1004  FORMAT(' LAMBDA=',F8.4,5X,'ITERATION EPSILON = ',1PD8.2,/)
C
C     INITIALIZE.
      ITCNT=0
      DO 130 I=1,N
        X0(I)=0.0D0
130     B(I)=LAMBDA-ATAN(1.0D0-T(I))-ATAN(T(I))
      WRITE(IWRITE,1000)
1000  FORMAT(//,' ITERATION RESULTS.',//,3X,'M',6X,'DIFF(M)',4X,
     *       'RATIO',/)
C
C     ITERATION.
140   DO 160 I=1,N
        RESID(I)=0.0D0
        DO 150 J=1,N
150       RESID(I)=RESID(I)+KERMAT(I,J)*X0(J)
160     RESID(I)=B(I)-LAMBDA*X0(I)+RESID(I)
C
      CALL DGESL(A,MAXN,N,IPVT,RESID,0)
C
C     CALCULATE THE SIZE OF THE CORRECTION TERM.
      ERRMAX=0.0D0
      DO 170 J=1,N
170     ERRMAX=MAX(ERRMAX,ABS(RESID(J)))
      ITCNT=ITCNT+1
      IF(ITCNT .EQ. 1) THEN
          WRITE(IWRITE,1001) ITCNT,ERRMAX
      ELSE
          WRITE(IWRITE,1001)ITCNT,ERRMAX,ERRMAX/PASTEM
      END IF
1001  FORMAT(I4,D14.4,F10.5)
      PASTEM=ERRMAX
C
      DO 180 I=1,N
180     X0(I)=X0(I)+RESID(I)
      IF(ERRMAX .GT. EPSLON) GO TO 140
      GO TO 90
      END

      FUNCTION FCNK(T)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      FCNK=1.0D0/(1.0D0+T*T)
      RETURN
      END
