Fortran Gaussian Elimination in Fortran

[50Cent]

Hi,

I am trying to recreate the naive gauss elimination method in fotran 95 but am having a few problems with it. The idea is to read in a nxn matrix of equations, so you can type in any number when u start the program and then the program will ask you to enter the relavant amount of coefficients.

This is my program so far:

PROGRAM Gauss1
! Solution of a system of N linear equations
! A.X = C
IMPLICIT NONE
PARAMETER (IN=20)
REAL:: A(IN,IN), X(IN), C(IN), F, SUM
INTEGER:: K, I, J, N, IN

! Read or input A and C here

PRINT *,'Input the number of equations (N):'
READ (5,*) N

PRINT *,'Input the matrix coefficients [A(I,J)]:'
READ (5,*) ((A(I,J),J=1,N),I=1,N)

PRINT *,'Input the solution vector[C(I)]:'
READ (5,*) (C(IN),N=1,N)

! Decomposition (Elimination)
DO K = 1, N-1
DO I = K+1, N
F = A(I,K) / A(K,K)
DO J = K+1, N
A(I,J) = A(I,J) - F * A(K,J)
ENDDO
C(I) = C(I) - F * C(K)
ENDDO
ENDDO

! Back-substitution
X(IN) = C(IN) / A(IN,IN)
DO I = N-1, 1, -1
SUM = 0.0
DO J = I+1, N
SUM = SUM + A(I,J) * X(J)
ENDDO
X(I) = (C(I) - SUM) / A(I,I)
ENDDO

! Print out results X here
WRITE (6,*) 'x(1)= ',X(1), ' x(2)= ',X(2), ' x(3)= ',X(3)

STOP

END PROGRAM Gauss1

I know the last print command is incorrect because its only specific for 3 unknowns. I didnt know how to write it to print N unknowns.

At the moment the program compiles and runs, asks for the the coefficients but then stops there. Can anyone see the problem?

Any help is much appreciated
Thanks

 

[jtbell]

Does it ask for the solution vector?

If it does, then I would add some extra PRINT statements so as to see exactly how far the program is going. For example, first add one immediately before the back-substitution section. If you don't see that output, then add a PRINT statement inside the loops of the decomposition section, that displays the values of I, J and K. Otherwise add a similar PRINT statement inside the loops of the back-substitution section. Etc.

The general idea is to make the program itself tell you what it is actually doing, instead of trying to puzzle it out simply by reading the code.

 

[50Cent]

Hi, Thanks for the reply. Yes i do need to find the solution vector X(IN).

I think the problem is actually the last print statement that i have. The one i have used is for a set 3x3 matrix. Do you know how to code it so that it will print "N" solutions?

i tried this:

WRITE(6,*) (X(IN),IN=1,N)

but get the error:

Compiling file: Gauss1.f95
H:\ENG3068\Gauss1\Gauss1.F95(43) : error 227 - DO control variables must be INTEGER or REAL scalar variables only
H:\ENG3068\Gauss1\Gauss1.F95(43) : error 52 - Compilation abandoned
Compilation failed.

In the meanwhile i will add the extra print statements to check how far the program is going

Thanks

 

[Borek]

Hi, Thanks for the reply. Yes i do need to find the solution vector X(IN).

Hardly surprising That's not what jtbell asked about. He asked IF the program asks you to enter solution vector - you wrote

asks for the the coefficients but then stops there.

So we don't know if problems started BEFORE elimination loop, or if the program was able to get INTO elimination loop.

 

[jtbell]

Exactly.