C Solution of Poisson's equation -△u = f in the rectangle
C [XMIN,XMAX]×[YMIN,YMAX]  by Finite Difference Method.
C 与えられた長方形の一辺を IB+1 等分し, 差分法を適用する
C 解くべき連立一次方程式の係数行列のサイズは IB*IB となる
C IB の上限値はサブルーチン INVVEC の中のパラメータ IBMAX の値で
C 制約される.
C 章末課題 12.2 に対応するためグラフィック表示をやめ, 分割数を
C 自動的に 10, 20, 40, 80 と取って Richardson 補外法により真値を
C 予測し, 誤差のオーダーを見るように変更したもの
C 値の比較は中央値だけで行っている
      PROGRAM POISN2
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IBMAX=159,NMAX=IBMAX*IBMAX,IRCT=5)
      PARAMETER(XMIN=0,XMAX=1,YMIN=0,YMAX=1,ZMIN=-0.5,ZMAX=1)
      REAL*8 A(0:IBMAX,NMAX),U(NMAX),V(NMAX),T(0:IRCT-1,IRCT)
C
      IB=10
      DO 100 KC=1,IRCT
          N=(IB-1)**2
          HX=(XMAX-XMIN)/IB
          HHX=HX*HX
          HX1=1/HHX
          HY=(YMAX-YMIN)/IB
          HHY=HY*HY
          HY1=1/HHY
          HCD=(HX1+HY1)*2
          HCX=-HX1
          HCY=-HY1
          X=XMIN
          DO 5 I=0,IB-1
             DO 3 J=1,N
                 A(I,J)=0
    3        CONTINUE
    5     CONTINUE
          DO 20 I=0,IB-2
             X=X+HX
             Y=YMIN
             DO 10 J=1,IB-1
                 Y=Y+HY
                 K=I*(IB-1)+J
                 U(K)=F(X,Y)
                 A(0,K)=HCD
                 A(IB-1,K)=HCY
                 IF (J.EQ.IB-1) THEN
                    A(1,K)=0
                 ELSE
                    A(1,K)=HCX
                 END IF
   10        CONTINUE
   20     CONTINUE
          CALL INVVEC(A,IB-1,U,V,N)
          T(0,KC)=V(N/2+1)
          IB=IB*2
  100 CONTINUE     
      M=2
      DO 300 I=1,IRCT-1
         DO 200 I2=1,IRCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
  200    CONTINUE
         M=M*2
  300 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 1st order:'
      DO 500 I=0,IRCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
  500 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IRCT-1,1)),I=1,IRCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IRCT-1,1))/
     +            (T(0,I)-T(IRCT-1,1)),I=1,IRCT-1)
      M=4
      DO 3000 I=1,IRCT-1
         DO 2000 I2=1,IRCT-I
            T(I,I2)=(M*T(I-1,I2+1)-T(I-1,I2))/(M-1)
 2000    CONTINUE
         M=M*4
 3000 CONTINUE
      WRITE(*,*)'Richardson extrapolation of 2nd order:'
      DO 5000 I=0,IRCT-1
          WRITE(*,'(5F18.15)')(T(I-I2+1,I2),I2=1,I+1)
 5000 CONTINUE
      WRITE(*,*)'Error:'
      WRITE(*,'(4F19.15)')((T(0,I)-T(IRCT-1,1)),I=1,IRCT)
      WRITE(*,*)'Ratio:'
      WRITE(*,'(4F19.15)')((T(0,I+1)-T(IRCT-1,1))/
     +              (T(0,I)-T(IRCT-1,1)),I=1,IRCT-1)
      END
C
C      -1
C V = A  U のコレスキー法による計算ルーチン 
C ここは Gauss の消去法でも, Gauss-Seidel 法でも何でもよい.
C コレスキー法については, 本日の講義で解説する
      SUBROUTINE INVVEC(A,IB,U,V,N)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(IBMAX=159,NMAX=IBMAX*IBMAX)
      REAL*8 A(0:IBMAX,NMAX),U(N),V(N),
     +         W0(NMAX),W1(NMAX),CHOL(0:IBMAX,NMAX)
      DO 5 I=1,N
          DO 3 K=0,MIN(N-I,IB)
             CHOL(K,I)=A(K,I)
             DO 2 J=MAX(1,K+I-IB),I-1
                CHOL(K,I)=CHOL(K,I)-CHOL(I-J,J)*CHOL(K+I-J,J)
    2        CONTINUE
             IF (K.EQ.0) THEN
                CHOL(0,I)=DSQRT(CHOL(0,I))
             ELSE
                CHOL(K,I)=CHOL(K,I)/CHOL(0,I)
             END IF
    3     CONTINUE
    5 CONTINUE
      DO 20 I=1,N
         W0(I)=U(I)
         DO 10 J=MAX(1,I-IB),I-1
            W0(I)=W0(I)-CHOL(I-J,J)*W0(J)
   10    CONTINUE
         W0(I)=W0(I)/CHOL(0,I)
   20 CONTINUE
      DO 40 I=N,1,-1
         W1(I)=W0(I)
         DO 30 J=I+1,MIN(I+IB,N)
            W1(I)=W1(I)-CHOL(J-I,I)*W1(J)
   30    CONTINUE
         W1(I)=W1(I)/CHOL(0,I)
   40 CONTINUE
      DO 300 J=1,N
         V(J)=W1(J)
  300 CONTINUE
      RETURN
      END
C
C GIVEN RIGHT-HAND SIDE
      DOUBLE PRECISION FUNCTION F(X,Y)
      REAL*8 F,X,Y
C      IF ((DABS(X-0.5D0).LT.0.2D0).AND.(DABS(Y-0.5D0).LT.0.2D0)) THEN
C       F=5.0D0
C      ELSE
C       F=0.0D0
C      END IF
      F=5.0D0
      RETURN
      END