!###########################################################
! page.104-106 表5.1(例題4)
! ヤコビ法のプログラム
!###########################################################
PROGRAM prog19
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,DD
REAL(8),DIMENSION(N,N) :: aa,LU

!---------------
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
LU=0.d0
delta=1.d-8
!---------------

LU=aa
DO i=1,N
LU(i,i)=0.d0
DD(i)=aa(i,i)
END DO

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

DO k=1,KMAX

xx3=MATMUL(LU,xx1)
xx3=bb-xx3
DO i=1,N
xx2=xx3/DD
END DO

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

xx1=xx2
xx2=0.d0

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

100 CONTINUE
WRITE(*,*)
WRITE(*,*) "ヤコビ法で求めた解は"
WRITE(*,"('x=',1P5E14.6)") xx2

STOP
END PROGRAM prog19
!###########################################################