! EUC-jp
! 文字コードはEUC-jp
! 自由形式のfortranで書いています。
! g77でコンパイルする場合は、
!g77 -ffree-form -o gau gau.f
! のようにコンパイルしてください。
! gfortran でコンパイルする場合は拡張子を.f90にして、
!日本語出力がある場合は、文字コードをutf-8にしてください。
! ガウスの消去法で対称行列a(i,j)の
!連立1次方程式a(i,j)*x(j)=b(i)を解くサブルーチン
! gausu の使用例。
implicit real*8(a-h, o-z)
dimension a(3,3),x(3),b(3)
! 例1答え:x(1)=0.,x(2)=1.,x(3)=2.
a(1,1)=1.;  a(1,2)=0.;  a(1,3)=0.
a(2,1)=0.;  a(2,2)=2.;  a(2,3)=-1.
a(3,1)=0.;  a(3,2)=-1.;  a(3,3)=1.
b(1)=0.
b(2)=0.
b(3)=1.
!
!
n=3
!
do i=2,n
do j=1,i-1
a(i,j)=9. !下三角に嘘の値を入れておく
end do
end do
!
!print*,'a='
do i=1,n
print*, (a(i,j),j=1,n)
end do
!
!print*,'b='
!print*, (b(i),i=1,n)
!
call gausu(n,a,x,b)
print*,'x='
print*, (x(i),i=1,n)
end
!
subroutine gausu(n,a,x,b)
implicit real*8(a-h, o-z)
dimension a(n,n),x(n),b(n)
!
!ガウスの消去法の参考としたのは、
!名取亮「すうがくぶっくす12 線形計算」(朝倉書店)p.10-15
!
do k=1,n-1  !a(k,k)を消去
do i=k+1,n  !k+1行からn行まで
do j=i,n    !i列からn列まで(対角項から右側)
!print*,'+kij=',k,i,j
! k<i, k<j, i<j
   a(i,j)=a(i,j)-a(k,j)*a(k,i)/a(k,k)
!  a(i,j)=a(i,j)-a(k,j)*a(i,k)/a(k,k)
  !消すのはa(k,k)の下のa(i,k)だけど、これは常に下三角
  !にあるから上三角のa(k,i)を使う
end do
   b(i)=b(i)-b(k)*a(k,i)/a(k,k)
   !b(i)=b(i)-b(k)*a(i,k)/a(k,k)
  !消すのはa(k,k)の下のa(i,k)だけど、これは常に下三角
  !にあるから上三角のa(k,i)を使う
end do
end do
!
!
! 後退代入
x(n)=b(n)/a(n,n)
do k=n,1,-1
   akjxj=0.d0
do j=k+1,n
   akjxj=akjxj+a(k,j)*x(j)
end do
   x(k)=(b(k)-akjxj)/a(k,k)
end do
!
return
end
!
!

