algorithm - 将 VBA 代码转换为 Fortran

Question

我正在尝试通过部分旋转的高斯消去法求解这个方程。

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'。

谁能帮我找出我的算法有什么问题吗？？

Answer 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

做这些事情会消除很多错误。如果仍然有错误，您应该发布准确的错误消息，并指出它们引用了代码中的哪几行。