我正在尝试通过部分旋转的高斯消去法求解这个方程。
x-2y-z=2
5x+2y+2z=9
-3x+5y-z=1
所以我放
1 2 -1
5 2 2
-3 5 -1
到 INPUT1.DAT 和
2
9
1
到 INPUT2.DAT。
这是运行良好的 VBA 代码，
Option Explicit
Sub GaussElim()
Dim n As Integer, er As Integer, i As Integer
Dim a(10, 10) As Double, b(10) As Double, x(10) As Double
n = 3
a(1, 1) = 1: a(1, 2) = 2: a(1, 3) = -1
a(2, 1) = 5: a(2, 2) = 2: a(2, 3) = 2
a(3, 1) = -3: a(3, 2) = 5: a(3, 3) = -1
b(1) = 2: b(2) = 9: b(3) = 1
Call Gauss(a, b, n, x, er)
If er = 0 Then
For i = 1 To n
MsgBox "x(" & i & ") = " & x(i)
Next i
Else
MsgBox "ill-conditioned system"
End If
End Sub
Sub Gauss(a, b, n, x, er)
Dim i As Integer, j As Integer
Dim s(10) As Double
Const tol As Double = 0.000001
er = 0
For i = 1 To n
s(i) = Abs(a(i, 1))
For j = 2 To n
If Abs(a(i, j)) > s(i) Then s(i) = Abs(a(i, j))
Next j
Next i
Call Eliminate(a, s, n, b, tol, er)
If er <> -1 Then
Call Substitute(a, n, b, x)
End If
End Sub
Sub Pivot(a, b, s, n, k)
Dim p As Integer, ii As Integer, jj As Integer
Dim factor As Double, big As Double, dummy As Double
p = k
big = Abs(a(k, k) / s(k))
For ii = k + 1 To n
dummy = Abs(a(ii, k) / s(ii))
If dummy > big Then
big = dummy
p = ii
End If
Next ii
If p <> k Then
For jj = k To n
dummy = a(p, jj)
a(p, jj) = a(k, jj)
a(k, jj) = dummy
Next jj
dummy = b(p)
b(p) = b(k)
b(k) = dummy
dummy = s(p)
s(p) = s(k)
s(k) = dummy
End If
End Sub
Sub Substitute(a, n, b, x)
Dim i As Integer, j As Integer
Dim sum As Double
x(n) = b(n) / a(n, n)
For i = n - 1 To 1 Step -1
sum = 0
For j = i + 1 To n
sum = sum + a(i, j) * x(j)
Next j
x(i) = (b(i) - sum) / a(i, i)
Next i
End Sub
Sub Eliminate(a, s, n, b, tol, er)
Dim i As Integer, j As Integer, k As Integer
Dim factor As Double
For k = 1 To n - 1
Call Pivot(a, b, s, n, k)
If Abs(a(k, k) / s(k)) < tol Then
er = -1
Exit For
End If
For i = k + 1 To n
factor = a(i, k) / a(k, k)
For j = k + 1 To n
a(i, j) = a(i, j) - factor * a(k, j)
Next j
b(i) = b(i) - factor * b(k)
Next i
Next k
If Abs(a(k, k) / s(k)) < tol Then er = -1
End Sub
我尝试将此代码转换为 Fortran，如下所示，
program Gauss_Emlimination !with partial pivoting
implicit none
INTEGER n, i, j
REAL A(3,3), B(3), X(3), ER, tol
tol = 0.000001
n = 3
OPEN(UNIT=2, FILE='INPUT1.DAT')
OPEN(UNIT=3, FILE='INPUT2.DAT')
OPEN(UNIT=4, FILE='RESULT.DAT')
READ(2,*) ((A(I,J),J=1,3),I=1,3)
READ(3,*) (B(I), I=1,3)
CALL Gauss(a, b, n, x, er)
IF (er .EQ. 0) THEN
DO i =1, N
WRITE(4,*) X(i)
END DO
ELSE
WRITE(4,*) "ill-conditioned system"
END IF
contains
SUBROUTINE Gauss(a, b, n, x, er)
real, intent(inout) :: a(3,3)
real, intent(inout) :: b(3)
integer, intent(in) :: n
real, intent(out) :: x(N)
REAL, intent(out) :: er
real, dimension(10) :: S(10)
INTEGER I, J
ER=0
DO I= 1, N
s(i) = ABS(A(i,1))
DO j = 2, n
IF (ABS(A(i,j)) .GT. s(i)) THEN
s(i) = ABS(A(i,j))
END IF
END DO
END DO
CALL Eliminate(a, s, n, b, tol, er)
IF (er .ne. -1) THEN
CALL Substitute(a, n, b, x)
END IF
END SUBROUTINE Gauss
SUBROUTINE Pivot(a, b, s, n, k)
INTEGER ii, jj
real, intent(inout) :: a(3,3)
real, intent(inout) :: b(3)
integer, intent(in) :: n, K
integer p
real, dimension(10) :: S(10)
DOUBLE PRECISION big, dummy, factor
p = k
big = ABS(A(k,k)/S(k))
DO ii = k+1, n
dummy = ABS(A(ii, k)/S(ii))
IF (dummy .GT. big) THEN
big = dummy
p = ii
END IF
END DO
IF (p .ne. k) THEN
DO jj = k, n
dummy = A(p, jj)
A(p, jj) = A(k, jj)
A(k, jj) = dummy
END DO
dummy = B(p)
B(p) = B(k)
B(k) = dummy
dummy = S(p)
S(p) = S(k)
S(k) = dummy
END IF
END SUBROUTINE Pivot
SUBROUTINE Substitute(a, n, b, x)
INTEGER i, j
real, intent(inout) :: a(3,3)
real, intent(inout) :: b(3)
integer, intent(in) :: n
real, intent(out) :: x(N)
DOUBLE PRECISION sum
X(n) = B(n)/A(n, n)
DO i = n-1, 1, -1
sum = 0
DO j = i+1, n
sum = sum +A(i, j)*X(j)
END DO
X(n) = (B(n)-sum)/A(n,n)
END DO
END SUBROUTINE Substitute
SUBROUTINE Eliminate(a, s, n, b, tol, er)
real, intent(in) :: tol
real, intent(inout) :: a(3,3)
real, intent(inout) :: b(3)
integer, intent(in) :: n
real, dimension(10) :: S(10)
real, intent(INout) :: er
INTEGER i, j, k
DOUBLE PRECISION factor
DO K = 1, N-1
CALL Pivot (a, b, s, n, k)
IF (ABS(A(K,K)/S(K)) .LT. tol) THEN
er=-1
EXIT
END IF
DO i = k+1, n
factor = A(i,k)/A(k,k)
DO j= k+1, n
A(i,j) = A(i,j) - factor*B(k)
END DO
B(i) = B(i) - factor * B(k)
END DO
END DO
IF (ABS(A(n,n)/S(n)) .LT. tol) THEN
er= -1
END IF
END SUBROUTINE Eliminate
end program Gauss_Emlimination
而且这段代码没有错误。
但问题是我得到了 ' 0.0000000E+00 0.0000000E+00 -7.1424372E-02' 结果。
它应该是 'x(1)=1, x(2)=1, x(3)=1'。
谁能帮我找出我的算法有什么问题吗？？
最佳答案
首先，您应该确保将所有这些子例程都放在一个模块中。这样你就不需要在每个子例程中声明 External GaussElim，因为编译器会知道模块中的所有子例程，以及它们期望的参数。为此，只需将所有这些子例程放在一个文件中并将它们放在中间:
module gauss_mod
implicit none
contains
! your code here
end module gauss_mod
然后在你的主程序中，只需将use gauss_mod 放在顶部，你就可以访问模块中的所有子程序。
implicit none 告诉您的编译器您将声明所有变量，并且它不应该猜测您没有告诉它的任何类型。例如，这会捕获很多由打字错误引起的错误。
其次，您需要声明子程序的参数。这是导致大部分错误的原因。在 GaussElim 之外，其他子程序都不知道像 A 这样的变量是什么。结果，当你的编译器看到
s(i) = ABS(A(i,1))
它认为 A(i,1) 是一个函数，并给出与此相关的错误。您可以通过简单地将以下行添加到您的子例程来修复它:
double precision, dimension(:,:) :: A
这告诉子例程 A 必须有两个维度，但可以是任意大小。
您还应该为输入参数添加 intent(in)，为输出参数添加 intent(out)，为更改的参数添加 intent(inout)通过你的子程序。
此外，不使用 double ，而是使用real 并设置kind 参数:
module gauss_mod
implicit none
integer, parameter :: dp = selected_real_kind(15)
contains
! your code here
! As an example:
subroutine gauss(a, b, n, x, er)
! dummy arguments
real(kind=dp), dimension(:,:), intent(in) :: a, b
integer, intent(in) :: n
real(kind=dp), dimension(:), intent(out) :: x
real(kind=dp), intent(out) :: er
real(kind=dp), dimension(10) :: S(10)
real(kind=dp), parameter :: tol = 0.000001
! rest of subroutine
end subroutine gauss
end module gauss_mod
做这些事情会消除很多错误。如果仍然有错误，您应该发布准确的错误消息，并指出它们引用了代码中的哪几行。
关于algorithm - 将 VBA 代码转换为 Fortran，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/23006816/