Converting VBA code to Fortran

Question

I am trying to solve this equation with Gauss elimination with partial pivoting.

x-2y-z=2
5x+2y+2z=9
-3x+5y-z=1

so I put

 1    2   -1

 5    2    2

-3    5   -1

to INPUT1.DAT and

2
9
1

to INPUT2.DAT.

This is VBA code running well,

   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

and I tried to convert this code to Fortran like below,

      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

and I got no error for this code.

but the problem is that I got ' 0.0000000E+00 0.0000000E+00 -7.1424372E-02' as a result.

it supposed to be 'x(1)=1, x(2)=1, x(3)=1'.

Could anyone help me to find what is wrong in my algorithm please??

I suggest you highlight specific lines that you get errors on and what those error messages were... — Dan, Apr 11 '14 at 08:27
Your first code looks like VBA, not Matlab, to me. If it is, please edit the question and title. — Fortranner, Apr 11 '14 at 12:14
@kol This is my school subject.. Professor made me to do it. — leaving_traces, Apr 11 '14 at 13:02
You now have lots of reasonably helpful compiler errors :) The number in brackets after the file name is the line number where the error is. At the end of some of the error messages, in square brackets, is the variable which is the cause of that error. Read each error carefully and compare it to the line that it refers to. Try and fix one error at a time and then recompile. — Yossarian, Apr 11 '14 at 15:31
@Yossarian Hoorrayy!! I got 0 error with revised code as you can see above. but the problem is that the answer is wrong. it supposed to be like "x(1)=1, x(2)=1, x(3)=1" but I got " 0.0000000E+00 0.0000000E+00 -7.1424372E-02"... Can you find what's wrong in my algorithm? — leaving_traces, Apr 12 '14 at 13:12
@leaving_traces That's great. If you feel my answer helped, you can click accept. I'm afraid I'm not familiar with Gauss elimination. You may find it helpful to write out variables during intermediate parts of your code using `print*, variable_name`. This way you can see whether things have sensible values or not. — Yossarian, Apr 14 '14 at 07:55

score 3 · Accepted Answer · answered Apr 11 '14 at 08:52

Firstly, you should make sure you put all these subroutines in a module. That way you don't need to declare External GaussElim, for example, in each subroutine as the compiler will know about all the subroutines in the module, and what arguments they expect. To do that, just put all these subroutines in one file and put them inbetween:

module gauss_mod

implicit none

contains

! your code here

end module gauss_mod

Then in your main program, you can just put use gauss_mod at the top, and you will have access to all the subroutines in the module. The implicit none tells your compiler that you are going to declare all variables, and it shouldn't guess the type of any you haven't told it about. This will catch a lot of errors caused by typos, for example.

Secondly, you need to declare the arguments to your subroutines. This is what is causing most of your errors. Outside of GaussElim, none of the other subroutines know what variables like A are. As a result, when your compiler sees

s(i) = ABS(A(i,1))

It thinks A(i,1) is a function and gives you errors related to that. You can fix it by simply adding the following line to your subroutines:

double precision, dimension(:,:) :: A

This tells the subroutine that A must have two dimensions, but can be any size. You should also add intent(in) for input arguments, intent(out) for output arguments and intent(inout) for arguments which are changed by your subroutines.

Additionally, instead of using double precision, use real and set the kind parameter:

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

Doing these things will get rid of lots of the errors. If you still have errors, you should post the exact error message, and indicate which lines in the code they refer to.

Thank you for helping me a lot.. but I still getting same amount of error messages. Could you help me a little bit more? I have edited my original post to advised code. — leaving_traces, Apr 11 '14 at 13:05
If you are calling these subroutines from a main program, do you have a `use` statement to cite the module? — M. S. B., Apr 11 '14 at 13:08
@M.S.B. As you can see in the code, I just wanted to get solution for 3 unknowns polynomial equations by inputting coefficients by INPUT1.DAT, INPUT2.DAT using these subroutines only. I don't need main program.. — leaving_traces, Apr 11 '14 at 14:01
@leaving_traces In fortran, you always need a `program`. In your case, you can easily turn the subroutine `GaussElim` into it. You may find it helpful to read a fortran tutorial. — Yossarian, Apr 11 '14 at 15:34
Yes, you must have a main program unless the subroutines are being called from another language. See http://stackoverflow.com/questions/6511711/computing-the-cross-product-of-two-vectors-in-fortran-90 for an example of how to organize your code. I think you need `end subroutine` rather than end. That should clear up some of the errors. Please edit after that with the revised code and the errors that remain (if any). — M. S. B., Apr 11 '14 at 16:04

Converting VBA code to Fortran

1 Answers1