program Gauss_elimination
implicit real(8) (a-h,o-z)
allocatable:: a(:,:),b(:),x(:),a_bak(:,:),b1(:),xq(:),xz(:)
n=4
allocate(a(n,n),b(n),x(n),a_bak(n,n),b1(n),xq(n),xz(n))
a(1,1)=1.1348d0; a(1,2)=3.8326d0; a(1,3)=1.1651d0; a(1,4)=3.4017d0;
a(2,1)=0.5301d0; a(2,2)=1.7875d0; a(2,3)=2.5330d0; a(2,4)=1.5435d0;
a(3,1)=3.4129d0; a(3,2)=4.9317d0; a(3,3)=8.7643d0; a(3,4)=1.3142d0;
a(4,1)=1.2371d0; a(4,2)=4.9998d0; a(4,3)=10.6721d0; a(4,4)=0.0147d0;
b(1)=9.5342d0; b(2)=6.3941d0; b(3)=18.4231d0; b(4)=16.9237d0;
a_bak=a; b1=b
!-----消去过程-----
do i=1,n
xq(i)=i
enddo
do i=1,n
!-----找第i列绝对值最大的a(i_rec,i)----
tmpMax=0.d0
do ic=i,n
do jc=i,n
if(tmpMax<dabs(a(ic,jc))) then
tmpMax=dabs(a(ic,jc))
i_rec=ic
j_rec=jc
endif
enddo
enddo
!----把所找到的第i_rec列与第i列进行交换
if(i_rec.ne.i) then
do jca=i,n
tmp=a(i,jca)
a(i,jca)=a(i_rec,jca)
a(i_rec,jca)=tmp
enddo
tmp=b(i)
b(i)=b(i_rec)
b(i_rec)=tmp
endif
if(j_rec.ne.i) then
do ica=1,n
tmp=a(ica,i)
a(ica,i)=a(ica,j_rec)
a(ica,j_rec)=tmp
enddo
endif
do k=n,1,-1;
if(j_rec.ne.i) then
tmp=x(i)
x(i)=x(j_rec)
x(j_rec)=tmp
t=xq(i)
xq(i)=xq(j_rec)
xq(j_rec)=t
endif
enddo
if(dabs(a(i,i))<=1.d-6) then
write(*,*)'主元太小',i
stop
endif
do L=i+1,n
am=a(L,i)/a(i,i)
do k=i,n
a(L,k)=a(L,k)-am*a(i,k)
enddo
b(L)=b(L)-am*b(i)
enddo
enddo
!-----回代过程----
x(n)=b(n)/a(n,n)
do i=n-1,1,-1
tmp=0.d0
do j=n,i+1,-1
tmp=tmp+a(i,j)*x(j)
enddo
x(i)=(b(i)-tmp)/a(i,i)
enddo
!
do i=1,n
xz(xq(i))=x(i)
enddo
write(*,*),xz(1),xz(2),xz(3),xz(4)
!-----验算------
b=matmul(a_bak,xz)
do i=1,n
write(*,*)'b1-b',i,b(i)-b1(i)
enddo
end
