PROJECT
(c) Write a program
for each of the above methods
using FORTRAN 90 as programming
language. - (7 Marks)
Fortran 90 code for solving
linear equations using Gauss
Elimination without pivoting
technique, gauss.f90
! ---------------------PROGRAM
GAUSSIAN
! gauss.f90
!
PROGRAM OF GAUSSIAN
ELIMINATION
PARAMETER (IN=20)
REAL:: A(IN,IN), B(IN)
PRINT *,'INPUT NUMBER
OF EQUATIONS N='
READ (5,*) N
PRINT *,'INPUT MATRIX
COEFFICIENTS A(I,J)='
READ (5,*) ((A(I,J),J=1,N),I=1,N)
PRINT *,'INPUT RIGHT-HAND
SIDE VECTOR B(I)='
READ (5,*) (B(I),I=1,N)
WRITE (6,*) ('****GAUSSIAN
ELIMINATION ****')
WRITE (6,*)
WRITE (6,*) ('COEFFICIENT
MATRIX you inputed:')
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) ('right hand
side vector you inputed:')
CALL PRINTV(B,N,6)
WRITE(6,*)
!
CONVERT TO UPPER TRIANGULAR
FORM
DO K = 1,N-1
IF (ABS(A(K,K)).GT.1.E-6)
THEN
DO I = K+1, N
X = A(I,K)/A(K,K)
DO J = K+1, N
A(I,J) = A(I,J) -A(K,J)*X
ENDDO
B(I) = B(I) - B(K)*X
ENDDO
ELSE
WRITE (6,*) 'ZERO PIVOT
FOUND IN LINE:'
WRITE (6,*) K
STOP
END IF
ENDDO
WRITE(6,*) 'MODIFIED
MATRIX'
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) 'MODIFIED
RIGHT HAND SIDE VECTOR'
CALL PRINTV (B,N,6)
WRITE(6,*)
!
BACK SUBSTITUTION
DO I = N,1,-1
SUM = B(I)
IF (I.LT.N) THEN
DO J= I+1,N
SUM = SUM - A(I,J)*B(J)
ENDDO
END IF
B(I) = SUM/A(I,I)
ENDDO
!
PRINT THE RESULTS
write(6,*) ('SOLUTION
VECTOR')
CALL PRINTV(B,N,6)
!
END
PROGRAM GAUSSIAN
!------------------------------------------
SUBROUTINE
PRINTA(A,IA,M,N,ICH)
!
!
WRITE A 2D ARRAY TO
OUTPUT CHANNEL 'ICH'
!
REAL A(IA,*)
DO I =1,M
WRITE(ICH,2) (A(I,J),J=1,N)
ENDDO
2
FORMAT(1X,6E12.4)
!
END
SUBROUTINE PRINTA
!-----------------------------------------
SUBROUTINE
PRINTV(VEC,N,ICH)
!
!
WRITE A COLUMN VECTOR
TO CHANNEL 'ICH'
!
REAL VEC(*)
WRITE(ICH,1) (VEC(I),I=1,N)
1
FORMAT(1X,6E12.4)
!
END
SUBROUTINE PRINTV
!-----------------------------------------
|
Compiling, input data, and
executing the objected code
%f90 -o run gauss.f90
%run
INPUT NUMBER OF EQUATIONS N=3
INPUT MATRIX COEFFICIENTS A(I,J)=
10.0 1.0 -5.0
-20.0 3.0 20.0
5.0 3.0 5.0
INPUT RIGHT-HAND SIDE VECTOR
B(I)=1.0 2.0 6.0
The sample of final results
can be displayed as
**** GAUSSIAN ELIMINATION
****
| Coefficient
matrix you inputed: |
| |
.1000E+02
.1000E+01 -.5000E+01
-.2000E+02 .3000E+01
.2000E+02
.5000E+01 .3000E+01
.5000E+01 |
| right
hand side vector
you inputed: |
| |
.1000E+01 .2000E+01
.6000E+01 |
Modified Matrix: |
| |
.1000E+02 .1000E+01
-.5000E+01
-.2000E+02 .5000E+01
.1000E+02
.5000E+01 .2500E+01
.2500E+01 |
| Modified
right hand side
vector: |
| |
.1000E+01 .4000E+01
.3500E+01 |
| Solution
vector: |
| |
.1000E+01 -.2000E+01
.1400E+01 |
|
If you want to input data from
a file rather than key strikes,
you can create a data file,
say, gauss.dat:
3
10.0 1.0 -5.0
-20.0 3.0 20.0
5.0 3.0 5.0
1.0 2.0 6.0
|
Then you can
use the following commands.
%f90 -o run gauss.f90
%run<gauss.dat
|
Fortran 90 code for
solving linear equations using
Gauss-Seidel iterative method
PROGRAM GAUSS_SEIDEL
! Gauss-Seidel method
for linear algebraic
euqations
! gauss_seidel.f90
!
PARAMETER (IN=50)
REAL A(IN,IN), B(IN),
X(IN), XNEW(IN), U(IN,IN)
PRINT *,'INPUT NUMBER
OF EQUATIONS N='
READ (5,*) N
PRINT *,'INPUT MATRIX
COEFFICIENTS A(I,J)='
READ (5,*) ((A(I,J),J=1,N),I=1,N)
PRINT *,'INPUT RIGHT-HAND
SIDE VECTOR B(I)='
READ (5,*) (B(I),I=1,N)
WRITE(*,*) 'Input
your gessed values:'
READ(5,*) (X(I),I=1,N)
PRINT *,'Input your
tolence and max iteration:'
READ(5,*) TOL,ITS
WRITE(6,*) ('******Gauss-Seidel
Method******')
WRITE(6,*)
WRITE(6,*) ('COEFFICIENT
MATRIX')
WRITE(6,*)
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) ('RIGHT
HAND VECTOR')
WRITE(6,*)
CALL PRINTV(B,N,6)
WRITE(6,*)
DO I=1,N
DIAG=A(I,I)
DO J=1,N
A(I,J)=A(I,J)/DIAG
ENDDO
B(I)=B(I)/DIAG
ENDDO
CALL NULL(U,IN,N,N)
DO I=1,N
DO J=I+1,N
U(I,J)=-A(I,J)
A(I,J)=0.0
ENDDO
ENDDO
WRITE(6,*) ('FIRST FEW
ITERATIONS')
ITERS=0
4 ITERS=ITERS+1
CALL MVMULT(U,IN,X,N,N,XNEW)
CALL VECADD(B,XNEW,XNEW,N)
CALL SUBFOR(A,IN,XNEW,N)
CALL CHECON(XNEW,X,N,TOL,ICON)
IF (ITERS.LE.5) CALL
PRINTV(X,N,6)
IF (ICON.EQ.0.AND.ITERS.LT.ITS)
GO TO 4
WRITE(6,*)
WRITE(6,*) ('ITERATIONS
TO CONVERGENCE')
WRITE(6,*) ITERS
WRITE(6,*)
WRITE(6,*) ('SOLUTION
VECTOR')
CALL PRINTV(X,N,6)
END PROGRAM GAUSS_SEIDEL
SUBROUTINE CHECON(LOADS,OLDLDS,N,TOL,ICON)
REAL LOADS(*),OLDLDS(*)
ICON=1
BIG=0.
DO I=1,N
IF (ABS(LOADS(I)).GT.BIG)
BIG=ABS(LOADS(I))
ENDDO
DO I=1,N
IF (ABS(LOADS(I)-OLDLDS(I))/BIG.GT.TOL)
ICON=0
OLDLDS(I) = LOADS(I)
ENDDO
END SUBROUTINE CHECON
SUBROUTINE NULL(A,IA,M,N)
REAL A(IA,*)
DO I=1,M
DO J=1,N
A(I,J)=0.0
ENDDO
ENDDO
END SUBROUTINE NULL
SUBROUTINE PRINTA(A,IA,M,N,ICH)
REAL A(IA,*)
DO I=1,M
WRITE (ICH,2) (A(I,J),J=1,N)
ENDDO
2 FORMAT (1X,6E12.4)
END SUBROUTINE PRINTA
SUBROUTINE PRINTV(VEC,N,ICH)
REAL VEC(*)
WRITE(ICH,1) (VEC(I),I=1,N)
1 FORMAT (1X,6E12.4)
END SUBROUTINE PRINTV
SUBROUTINE MVMULT(M,IM,V,K,L,Y)
REAL M(IM,*),V(*),Y(*)
DO I=1,K
X=0.
DO J=1,L
X=X+M(I,J)*V(J)
ENDDO
Y(I)=X
ENDDO
END SUBROUTINE MVMULT
SUBROUTINE VECADD(A,B,C,N)
REAL A(*),B(*),C(*)
DO I=1,N
C(I)=A(I)+B(I)
ENDDO
END SUBROUTINE VECADD
SUBROUTINE SUBFOR(A,IA,B,N)
REAL A(IA,*), B(*)
DO I=1,N
SUM=B(I)
IF (I.GT.1) THEN
DO J=1,I-1
SUM=SUM-A(I,J)*B(J)
ENDDO
END IF
B(I)=SUM/A(I,I)
ENDDO
END SUBROUTINE SUBFOR
SUBROUTINE MSMULT(A,IA,C,M,N)
REAL A(IA,*)
DO I=1,N
DO J=1,N
A(I,J)=A(I,J)*C
ENDDO
ENDDO
END SUBROUTINE
MSMULT
|
Compiling, input data,
and executing the objected code
%f90 -o run gauss_seidel.f90
%run
INPUT NUMBER OF EQUATIONS
N=
3
INPUT MATRIX COEFFICIENTS
A(I,J)=
16 4 8
4 5 -4
8 -4 22
INPUT RIGHT-HAND SIDE
VECTOR B(I)=
4 2 5
Input your gessed
values:
1.0 1.0 1.0
Input your tolence
and max iteration:
1.e-5
100
|
The sample of
final results can be displayed
as
******Gauss-Seidel Method******
COEFFICIENT MATRIX
.1600E+02 .4000E+01
.8000E+01
.4000E+01 .5000E+01
-.4000E+01
.8000E+01 -.4000E+01
.2200E+02
RIGHT HAND VECTOR
.4000E+01 .2000E+01
.5000E+01
FIRST FEW ITERATIONS
-.5000E+00 .1600E+01
.7000E+00
-.5000E+00 .1360E+01
.6564E+00
-.4182E+00 .1260E+01
.6084E+00
-.3691E+00 .1182E+01
.5764E+00
-.3337E+00 .1128E+01
.5537E+00
ITERATIONS TO CONVERGENCE
30
SOLUTION VECTOR
-.2500E+00 .1000E+01
.5000E+00
|
|
