!###########################################################
! page.106-108 表5.2(例題5)
! ガウス・ザイデル法のプログラム
!###########################################################
PROGRAM prog20
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: KMAX=100,N=5
REAL(8),DIMENSION(N) :: bb,xx1,xx2,xx3
REAL(8),DIMENSION(N,N) :: aa,LL,UU,DD

!---------------
aa(1,1)=10.d0 ; aa(1,2)= 3.d0 ; aa(1,3)= 1.d0 ; aa(1,4)= 2.d0 ; aa(1,5)= 1.d0
aa(2,1)= 1.d0 ; aa(2,2)=19.d0 ; aa(2,3)= 2.d0 ; aa(2,4)=-1.d0 ; aa(2,5)= 1.d0
aa(3,1)=-1.d0 ; aa(3,2)= 1.d0 ; aa(3,3)=30.d0 ; aa(3,4)= 1.d0 ; aa(3,5)=10.d0
aa(4,1)=-2.d0 ; aa(4,2)= 0.d0 ; aa(4,3)= 1.d0 ; aa(4,4)=20.d0 ; aa(4,5)= 5.d0
aa(5,1)=-3.d0 ; aa(5,2)= 5.d0 ; aa(5,3)= 1.d0 ; aa(5,4)=-2.d0 ; aa(5,5)=25.d0

bb(1)=-22.d0
bb(2)= 27.d0
bb(3)= 89.d0
bb(4)=-73.d0
bb(5)= 22.d0

xx1=0.d0
xx2=0.d0
xx3=0.d0
DD=0.d0
LL=0.d0
UU=0.d0
delta=1.d-8
!---------------

DO i=1,N
 DO j=1,N
     IF(i>j)THEN
     LL(i,j)=aa(i,j)
     ELSEIF(i<j)THEN
     UU(i,j)=aa(i,j)
     ELSE
     DD(i,j)=aa(i,j)
     END IF
 END DO
END DO

aa=DD+LL
! (DD+LL) の逆行列を求める。
CALL gyaku(N,aa)


k=0
WRITE(*,"('k=',i2,' x=',1P5E14.6)") k,xx1

DO k=1,KMAX

xx3=bb-MATMUL(UU,xx1)
xx2=MATMUL(aa,xx3)

WRITE(*,"('k=',i2,' x=',1P5E14.6)") k,xx2
fsum=SUM((xx2-xx1)**2)
IF(fsum<delta)GOTO 100

xx1=xx2
xx2=0.d0
xx3=0.d0

END DO
WRITE(*,*) '解が求まりませんでした'
STOP

100 CONTINUE
WRITE(*,*)
WRITE(*,*) "ガウス・ザイデル法で求めた解は"
WRITE(*,"('x=',1P5E14.6)") xx2

STOP
END PROGRAM prog20
!###########################################################



!************************************
!         inverse matrix
!************************************
SUBROUTINE gyaku(N,aa)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,INTENT(IN) :: N
REAL(8),DIMENSION(N,N),INTENT(INOUT) :: aa
REAL(8),DIMENSION(N,N):: aa1,b

aa1=aa
b=0.d0

DO i=1,N
b(i,i)=1.d0
END DO

DO i=1,N

      aii=aa1(i,i)
      IF(aii==0.d0)THEN
      WRITE(*,*) 'ERROR'
      STOP
      END IF

      DO j=1,N
      b(i,j)=b(i,j)/aii
      IF(j<i) CYCLE
      aa1(i,j)=aa1(i,j)/aii
      END DO
IF(i==N)CYCLE

   DO ii=i+1,N
   p=aa1(ii,i)
      DO j=1,N
      b(ii,j)=b(ii,j)-b(i,j)*p
      IF(j<i) CYCLE
      aa1(ii,j)=aa1(ii,j)-aa1(i,j)*p
      END DO
   END DO

END DO


DO i=N,2,-1
   DO ii=i-1,1,-1
   p=aa1(ii,i)
      DO j=1,N
      b(ii,j)=b(ii,j)-b(i,j)*p
      END DO
   END DO
END DO

aa=b

RETURN
END SUBROUTINE gyaku
!************************************
!
!************************************