gauss7.f90

I suggest you either typed  this as it appears into the compiler, or as an easier way, just copy and paste it!

Note: I do not claim the program will 100% work, but I give it my best. If it doesn't, try tweaking it yourself, or email me about it, and I'll modify and post the improved version later. Good luck!

PROGRAM Linear_Systems
IMPLICIT NONE

INTEGER:: n,Status_allocation, i, j
double precision, ALLOCATABLE, DIMENSION(:,:) :: Linear_equation, Linear_equation2
double precision,ALLOCATABLE, DIMENSION(:) :: x,x2
LOGICAL :: singular
double precision, DIMENSION(:), ALLOCATABLE :: temp
double precision :: abspivot, multiplier
double precision, PARAMETER :: epsilon = 1E-14
INTEGER :: pivotrow, gaussjumper=1

PRINT*, " WELCOME!"
PRINT*, " Triangular Factorization and Gaussian Elimination Compared"
PRINT*, " By: ADCY Software Sdn. Bhd."
PRINT*, ""

WRITE (*, '(1X, A)', Advance="NO") "How many equations to be calculated? "
READ*, n
!To exit early
IF (n==0) THEN
PRINT*, "Exiting..."
PRINT*, ""
STOP
END IF

ALLOCATE(Linear_equation(n,n+1), X(n), X2(n),temp(n+1), Linear_equation2(n,n+1),stat=Status_allocation)
IF (Status_allocation/=0) STOP 

Singular = .false.

!READ coefficient and constants
DO i=1,n
PRINT*, "Coefficients and constants?",i,":"
READ*, (Linear_equation(i,j), j=1,n+1)
END DO

print*, ""
do i=1,n
Linear_equation2(i,n+1)=Linear_equation(i,n+1) !for triangular factorization
Linear_equation2(i,i)=1 !for triangular factorization
end do
!no problem!


!triangular factorization to obtain LU

100 DO i=1,n

!locate pivot element

abspivot = ABS(Linear_equation(i,i))
pivotrow=i
DO j=1+1,n
IF (ABS(Linear_equation(j,i)) >abspivot) THEN
abspivot = ABS(Linear_equation(j,i))
pivotrow=j
END IF
END DO

!check IF matrix is singular

IF (abspivot<epsilon) THEN
Singular = .true.
GOTO 8
END IF

!it isn't, so interchange row pivotrow and i IF necessary

IF (pivotrow/=i) THEN
temp(i:n+1) = Linear_equation(i,i:n+1)
Linear_equation(i,i:n+1)=Linear_equation(pivotrow, i:n+1)
Linear_equation(pivotrow, i:n+1) = temp(i:n+1)
END IF

!eliminate ith unknown from equation i+1,... n

DO j=i+1,n
multiplier=Linear_equation(j,i)/Linear_equation(i,i)
Linear_equation(j,i:n+1)=Linear_equation(j,i:n+1)-multiplier*Linear_equation(i,i:n+1)
Linear_equation(j,i)=0
Linear_equation2(j,i)=multiplier
!Linear_equation2 is used to find values of y
END DO
! Linear_equation2(i,n+1)=Linear_equation(i,n+1)
END DO

!using gaussian elimination
if (gaussjumper==2) then
goto 150
end if

!using forward substitution

X2(1)= Linear_equation2(1,n+1)/Linear_equation2(1,1)
DO j=2,n
X2(j)=(Linear_equation2(j,n+1)-SUM(Linear_equation2(j,1:j-1)*X2(1:j-1)))/Linear_equation2(j,j)
END DO
!no problem so far

do i=1,n
Linear_equation(i,n+1)=Linear_equation2(i,n+1)
end do

!find solution by back substitution

150 X(n)= Linear_equation(n,n+1)/Linear_equation(n,n)
DO j=n-1,1,-1
X(j)=(Linear_equation(j,n+1)-SUM(Linear_equation(j,j+1:n)*X(j+1:n)))/Linear_equation(j,j)
END DO

8 IF (.not. Singular) THEN
!using gaussian elimination
if (gaussjumper==2) then
goto 300
end if

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!format used to print out matrix for viewing
! print*, ""
! do i=1,n
! print*, (Linear_equation(i,j), j=1,n+1)
! end do
!
! do i=1,n
! print*, (x2(i))
! end do
!
! print*, ""
! do i=1,n
! print*, (Linear_equation2(i,j), j=1,n+1)
! end do
!
!find roots
PRINT*, "The roots are (using Triangular factorization):"
goto 313
300 Print*, ""
PRINT*, "The roots are (using Gsussian Elimination):"
313 DO i=1,n
PRINT '(1X, "X(", I2, ") =", F20.16)', i, X(i)
END DO
!using gaussian elimination
if (gaussjumper==2) then
goto 400
end if


gaussjumper=gaussjumper+1
goto 100

400 PRINT*, ""
PRINT*, "Thank you for using GEAC and have a nice day"
PRINT*, "Another quality product from Dang Chwuan Yew"
PRINT*, "To order the product, call 90217, TODAY!"
PRINT*, ""
ELSE
PRINT*, "This matrix is nearly singular."
END IF

DEALLOCATE(Linear_equation,X)
END PROGRAM Linear_Systems