[Fortran] 纯文本查看 复制代码
! All-choice principal element Gaussian elimination method
! Solve the system of linear algebraic equations AX = B.
module SUBROUTINES
contains ! Used to encapsulate functions or subroutines
!implicit none ! Used to encapsulate variables
Subroutine SUB_agaus(FmatrixCoeff_, FvectorRside_, Iorder_, x, l, js)
integer(kind=4) :: row = 0, col = 0
integer(kind=4) :: Iorder_, l, js(:)
real(kind=8) :: b, t
real(kind=8), allocatable :: FmatrixCoeff_(:, :), FvectorRside_(:), x(:)
allocate(FmatrixCoeff_(Iorder_, Iorder_), FvectorRside_(Iorder_), x(Iorder_), js(Iorder_))
l = 1
Do k = 1, Iorder_ - 1
d = 0.0
! Do i = k, Iorder_
! Do j = k, Iorder_
! If (abs(FmatrixCoeff_(i,j))>d) Then
! d = abs(FmatrixCoeff_(i,j))
! js(k) = j
! is = i
! End If
! End Do
! End Do
forall ( row = k:Iorder_, col = k:Iorder_, abs(FmatrixCoeff_(row, col))>d )
d = abs(FmatrixCoeff_(row, col))
js(k) = col
is = row
end forall
If (d+1.0==1.0) Then
l = 0
Else
If (js(k)/=k) Then
Do i = 1, Iorder_
t = FmatrixCoeff_(i, k)
FmatrixCoeff_(i, k) = FmatrixCoeff_(i, js(k))
FmatrixCoeff_(i, js(k)) = t
End Do
End If
If (is/=k) Then
Do j = k, Iorder_
t = FmatrixCoeff_(k, j)
FmatrixCoeff_(k, j) = FmatrixCoeff_(is, j)
FmatrixCoeff_(is, j) = t
End Do
t = FvectorRside_(k)
FvectorRside_(k) = FvectorRside_(is)
FvectorRside_(is) = t
End If
End If
!If l=0, it means that the coefficient matrix of the system of equations is singular, and the solution fails; if L=0, it means that the normal return
If (l==0) Then
Write (*, 100)
Return
End If
Do j = k + 1, Iorder_
FmatrixCoeff_(k, j) = FmatrixCoeff_(k, j)/FmatrixCoeff_(k, k)
End Do
FvectorRside_(k) = FvectorRside_(k)/FmatrixCoeff_(k, k)
Do i = k + 1, Iorder_
Do j = k + 1, Iorder_
FmatrixCoeff_(i, j) = FmatrixCoeff_(i, j) - FmatrixCoeff_(i, k)*FmatrixCoeff_(k, j)
End Do
FvectorRside_(i) = FvectorRside_(i) - FmatrixCoeff_(i, k)*FvectorRside_(k)
End Do
End Do
If (abs(FmatrixCoeff_(Iorder_, Iorder_))+1.0==1.0) Then
l = 0
Write (*, 100)
Return
End If
x(Iorder_) = FvectorRside_(Iorder_)/FmatrixCoeff_(Iorder_, Iorder_)
Do i = Iorder_ - 1, 1, -1
t = 0.0
Do j = i + 1, Iorder_
t = t + FmatrixCoeff_(i, j)*x(j)
End Do
x(i) = FvectorRside_(i) - t
End Do
js(Iorder_) = Iorder_
Do k = Iorder_, 1, -1
If (js(k)/=k) Then
t = x(k)
x(k) = x(js(k))
x(js(k)) = t
End If
End Do
Return
100 Format (1X, ' FAIL ')
End Subroutine SUB_agaus
end module SUBROUTINES
!// ForQuill v1.01 Beta www.fcode.cn
! NOTE：
!
program PROGRAM_NAME
use SUBROUTINES, only : SUB_agaus
implicit none
integer ExitKey
integer(kind=4) :: i, n = 4, l
integer(kind=4), allocatable :: js(:) !Declare a variable-size array
real(kind=8), allocatable :: a(:, :), b(:), x(:) !Declare a variable-size array
allocate(js(n), a(n, n), b(n), x(n)) !Configure memory space
Data a/0.2368, 0.1968, 0.1582, 1.1161, 0.2471, 0.2071, 1.1675, 0.1254 &
, 0.2568, 1.2168, 0.1768, 0.1397, 1.2671, 0.2271, 0.1871, 0.1490/
Data b/1.8471, 1.7471, 1.6471, 1.5471/
Call SUB_agaus(a, b, n, x, l, js)
If (l /= 0) write(unit=6, fmt="('X(', I2, ' )=', D15.6)") (i, x(i), i=1, 4)
! Avoid having the command window exit before the results are output
write ( *, *) 'Press input any words to exit...'
read ( *, *) ExitKey
stop
end program PROGRAM_NAME