CODEWITHSUNDEEP
matrixfortrangaussiansolver
tackoo
tackoo

Reputation: 33

Basic Gauss Elimination solver yields wrong result

I'm a newbie to Fortran and I need to write a Gauss Elimination code to solve a 4x4 matrix. My code below returns a wrong result and I couldn't debug the problem. I'd much appreciate if you can help me.

      common /grid/ A(100,100), NEQ, C(100), X(100)

      open(10, file="NEQ.txt", status='unknown')
      read(10,*) NEQ
      close (10)

      open(12, file="C1.txt", status='unknown')
        do i=1,NEQ
        read(12,*) C(i)
        enddo
      close (12)

      open(11, file="A1.txt", status='unknown')
        do i=1,NEQ
        read(11,*) (A(i,k), k=1,NEQ)
        enddo
      close (11)

      call SOL

      open(13, file="X.txt", status='unknown')
        do i=1,NEQ
        write(13,*) X(i)
        enddo
      close (13)

      stop
      end

      subroutine SOL
      common /grid/ A(100,100), NEQ, C(100), X(100)


c     Forward Reduction Phase:

      do 10 K=2,NEQ
      do 10 I=K,NEQ
      R=A(I,K-1)/A(K-1,K-1)

      C(I)=C(I)-R*C(K-1)

      do 10 J=K-1,NEQ
10    A(I,J)=A(I,J)-R*A(K-1,J)

c     Back Substitution Phase:

      X(NEQ)=C(NEQ)/A(NEQ,NEQ)
      do 30 K=NEQ-1,1,-1
      X(K)=C(K)
      do 20 J=K+1,NEQ

20    X(K)=X(K)-A(K,J)*X(J)
30    X(K)=X(K)/A(K,K)

      return
      end

For my case NEQ is read from a text file as 4 And my A1.txt is:

18, -6, -6, 0
-6, 12, 0, -6
-6, 0, 12, -6
0, -6, -6, 18

And C1.txt is:

60
0
20
0

The resultant X matrix cames out to be:

8.3333330    
6.6666665    
8.3333321    
4.9999995

Instead of:

13.13
14.17
15.83
15.00

Upvotes: 2

Views: 1017

Answers (1)

sigma
sigma

Reputation: 2963

As Vladimir F has commented, it is extremely helpful to use at least Fortran 90 if you have to write any new code at all. It offers far better options than Fortran 77 for keeping your code structured and organised (and readable). I would like to add that implicit none is also an invaluable statement for reducing error.

That said, here is an example of what your algorithm might look like in Fortran 90 (I have written it as a function):

function gaussian_elimination(A, C) result(X)
  implicit none
  real, intent(inout) :: C(:), A(size(C), size(C))
  real :: X(size(C))

  real    :: R(size(C))
  integer :: i, j, neq

  neq = size(C)

  ! Forward reduction, only two loops since reduction is now row by row
  do i = 1, neq
    R = A(:,i)/A(i,i)

    do j = i+1, neq
      A(j,:) = A(j,:) - R(j)*A(i,:)
      C(j) = C(j) - R(j)*C(i)
    enddo
  enddo

  ! Back substitution, only one loop
  do i = neq, 1, -1
    x(i) = (C(i) - sum(A(i, i+1:) * x(i+1:))) / A(i,i)
  enddo
end function gaussian_elimination

This only leaves me to wonder why you think that your result X = [8.33, 6.67, 8.33, 5.00] is incorrect? It is the right solution, which you can verify by multiplying the matrix A with it: matmul(A, X) should be (nearly) equal to C.

Upvotes: 4

Related Questions