ガウス・ジョルダン法(部分ピボット選択あり)

 部分ピボット選択ありの場合、


module subprogs
implicit none
contains
subroutine gauss_jordan_pv(a0,x,b,n)!---部分ピボット選択あり---

integer ,intent(in)::n
real(8),intent(in)::a0(n,n),b(n)
real(8),intent(out)::x(n)
integer i,k,m
real(8) ar,am,t,a(n,n),w(n)
a(:,:)=a0(:,:)
x(:)=b(:)
do k=1,n
!---部分ピボット選択
m=k
am=abs(a(k,k))
do i=k+1,n
if(abs(a(i,k))>am)then
am=abs(a(i,k))
m=i
endif
enddo
if(am==0.0d0)stop 'A is singular'!Aが特異なら停止
if(k/=m)then
w( k:n)=a(k,k:n)
a(k,k:n)=a(m,k:n)
a(m,k:n)=w( k:n)
t=x(k)
x(k)=x(m)
x(m)=t
endif
if(a(k,k)==0.0d0)stop 'pivot = 0'
ar=1.0d0/a(k,k)
a(k,k)=1.0d0
a(k,k+1:n)=ar*a(k,k+1:n)
x(k)=ar*x(k)
do i=1,n
if(i/=k)then
a(i,k+1:n)=a(i,k+1:n)-a(i,k)*a(k,k+1:n)
x(i)=x(i)-a(i,k)*x(k)
a(i,k)=0.0d0
endif
enddo
enddo
end subroutine gauss_jordan_pv


subroutine set_random_ab(a,b,x,n)
implicit none
real(8),allocatable,intent(out)::a(:,:),b(:),x(:)
integer n
write(*,'(a)',advance='no')' input n: '
read (*,*)n
if(n<1.or.n>100)stop 'n must be 0 < n <101 '
allocate (a(n,n))
call random_number(a)
allocate (b(n))
call random_number(b)
allocate (x(n))
end subroutine set_random_ab

subroutine print_mat(a,b,n)
implicit none
real(8) a(n,n),b(n)
integer i,n
do i =1,n
write(*,'(101e12.4)') a(i,1:n),b(i)
enddo


end subroutine print_mat
end module subprogs


program main
use subprogs
implicit none
real(8),allocatable::a(:,:),b(:),x(:),r(:)
integer n,i
call set_random_ab(a,b,x,n)
call print_mat(a,b,n)
call gauss_jordan_pv(a,x,b,n)
do i=1,n
write(*,'(100e12.4)')x(i)
enddo
allocate (r(n))
r(:)=b(:)-matmul(a,x)
write(*,*) 'Gauss-Jordan error = ',dot_product(r,r)
deallocate(a,b,x)
end program main

ブログ気持玉

クリックして気持ちを伝えよう!

ログインしてクリックすれば、自分のブログへのリンクが付きます。

→ログインへ

なるほど(納得、参考になった、ヘー)
驚いた
面白い
ナイス
ガッツ(がんばれ!)
かわいい

気持玉数 : 0

この記事へのコメント

この記事へのトラックバック