***GAUSS-JORDAN METHOD************************************************
* THE PROGRAM FOR SOLVING SYSTEM OF LINEAR EQUATION BY GAUSS- *
* JORDAN METHOD WITH PARTIAL AND COMPLETE PIVOTING OPTION *
**********************************************************************
*
program GAUSSJORDN
implicit doubleprecision(a-h,o-z)
parameter (nd = 40)
doubleprecision A(nd,nd),INV(nd,nd),CM(nd),X(nd),Y(nd),ORD(nd)
$,ORDC(nd),SM(nd)
integer option
parameter (epsil = 1d-9)
*
print*, 'Select Type of Pivoting Required < 1 - PARTIAL PIVOT / 2
$- COMPLETE PIVOT > :'
read(*,*) option
write(6,*)
print*, 'Enter the number of equation, n:'
read(*,*) n
print*, 'Enter the coefficient matrix row-wise:'
read(*,*) ((A(i,j), j = 1,n), i = 1,n)
write(6,*)
*
m = n + 1
nn = 2*n
l = n - 1
c
c partial pivoting
if (option == 1) then
CALL PPIVOT(A,n,ORD)
endif
c
c complete pivoting
if (option == 2) then
CALL CPIVOT(A,n,ORD,ORDC)
endif
c
c augmenting the unit matrix
do i = 1, n
do j = m, nn
if(i == (j - n)) then
A(i,j) = 1
else
A(i,j) = 0
endif
enddo
enddo
c
c transforming the matrix upper triangular
do k = 1, l
kk = k + 1
do i = kk, n
qt = A(i,k)/A(k,k)
do j = k, nn
A(i,j) = A(i,j) - qt *A(k,j)
enddo
enddo
enddo
c
c checking for singularity of matrix
do i = 1, n
if (dabs(A(i,i)) < epsil) then
write(6,*)'----------------------------------------------'
write(6,*)'The coeff. matrix singular. No solution exist.'
write(6,*)'----------------------------------------------'
stop
endif
enddo
c
c making the diagonal elements unity
do i = 1, n
q = A(i,i)
do j = i, nn
A(i,j) = A(i,j) /q
enddo
enddo
c
c transformation of original matrix to unit matrix
do k = n, 1, -1
kk = k - 1
do i = kk, 1, -1
prod = A(i,k)
do j = k, nn
A(i,j) = A(i,j) - prod *A(k,j)
enddo
enddo
enddo
c
c storing the inverse
do i = 1, n
do j = 1, n
jj = j + n
INV(i,j) = A(i,jj)
enddo
enddo
*
write(6,*) '-------------------------------------------------'
write(6,*) ' The inverse of the pivoted matrix is:'
write(6,*)
CALL pr(INV,n,n)
write(6,*) '-------------------------------------------------'
*
icontinue = 1
10 if (icontinue == 1) then
write(6,*)
print*, 'Enter the right-hand side constant vector:'
read(*,*) (CM(i), i = 1,n)
c
c rearranging elements of constant vector
do i = 1, n
j = ORD(i)
SM(i) = CM(j)
enddo
*
CALL MATMUL(INV,SM,X,n,n,1)
c
c rearranging solution vector for complete pivoting
do i = 1, n
if (option == 2) then
j = ORDC(i)
Y(j) = X(i)
else
Y(i) = X(i)
endif
enddo
*
write(6,*) '-------------------------------------------------'
write(6,*) ' The solution vector is:'
write(6,*)
write(6,101) (Y(i), i = 1,n)
101 format('',29x,f12.4)
write(6,*)'--------------------------------------------------'
write(6,*) ' More right-hand side vector ? < 0 for NO / 1 for
$YES >:'
read(*,*) icontinue
go to 10
endif
write(6,*) '----------------------------------------------------'
*
stop
end
*
***MATRIX MULTIPLICATION**********************************************
*
subroutine MATMUL(A,B,PRO,m,l,n)
implicit doubleprecision(a-h,o-z)
parameter (nd = 40)
doubleprecision A(nd,nd),B(nd,nd),PRO(nd,nd)
*
do i = 1, m
do j = 1, n
PRO(i,j) = 0
do index = 1,l
PRO(i,j) = PRO(i,j) + A(i,index) *B(index,j)
enddo
enddo
enddo
*
return
end
*
***PARTIAL PIVOTING***************************************************
*
subroutine PPIVOT (A,n,ORD)
implicit doubleprecision(a-h,o-z)
parameter (nd = 40)
doubleprecision A(nd,nd),ORD(nd)
integer row,p
c
c establish initial order in order vector
do i = 1, n
ORD(i) = i
enddo
c
c partial pivoting
do i = 1, n
row = i
do p = i, n
if (dabs(A(row,i)) < dabs(A(p,i))) then
row = p
endif
enddo
*
if (row.ne.i) then
do j = 1, n
temp = A(i,j)
A(i,j) = A(row,j)
A(row,j) = temp
enddo
temp = ORD(i)
ORD(i) = ORD(row)
ORD(row) = temp
endif
enddo
*
return
end
*
***COMPLETE PIVOTING**************************************************
*
subroutine CPIVOT (A,n,ORD,ORDC)
implicit doubleprecision(a-h,o-z)
parameter (nd = 40)
doubleprecision A(nd,nd),ORD(nd),ORDC(nd)
integer col,row,p,r
c
c establish initial order in the order vector
do i = 1, n
ORD(i) = i
ORDC(i) = i
enddo
c
c complete pivoting
do i = 1, n
row = i
col = i
do p = i, n
do r = i, n
if (dabs(A(row,col)) < dabs(A(p,r))) then
row = p
col = r
endif
enddo
enddo
*
if (col.ne.i) then
do ii = 1, n
temp = A(ii,i)
A(ii,i) = A(ii,col)
A(ii,col) = temp
enddo
t = ORDC(i)
ORDC(i) = ORDC(col)
ORDC(col) = t
endif
*
m = n + 1
if (row.ne.i) then
do jj = 1, m
temp = A(i,jj)
A(i,jj) = A(row,jj)
A(row,jj) = temp
enddo
t = ORD(i)
ORD(i) = ORD(row)
ORD(row) = t
endif
enddo
*
return
end
*
***SUBROUTINE FOR PRINTING THE MATRIX*********************************
*
subroutine pr(C,m,n)
implicit doubleprecision(a-h,o-z)
parameter (nd = 40)
doubleprecision C(nd,nd)
*
if(n == 1) then
do i = 1,m
write(6,10)C(i,1)
10 format('',f10.4)
enddo
elseif(n == 2) then
do i = 1,m
write(6,20)C(i,1),C(i,2)
20 format('',2f10.4)
enddo
elseif(n == 3) then
do i = 1,m
write(6,30)C(i,1),C(i,2),C(i,3)
30 format('',3f10.4)
enddo
elseif(n == 4) then
do i = 1,m
write(6,40)C(i,1),C(i,2),C(i,3),C(i,4)
40 format('',4f10.4)
enddo
elseif(n == 5) then
do i = 1,m
write(6,50)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5)
50 format('',5f10.4)
enddo
elseif(n == 6) then
do i = 1,m
write(6,60)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6)
60 format('',6f10.4)
enddo
elseif(n == 7) then
do i = 1,m
write(6,70)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7)
70 format('',7f10.4)
enddo
elseif(n == 8) then
do i = 1,m
write(6,80)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7),
$C(i,8)
80 format('',8f10.4)
enddo
elseif(n == 9) then
do i = 1,m
write(6,90)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7),
$C(i,8),C(i,9)
90 format('',9f10.4)
enddo
elseif(n == 10) then
do i = 1,m
write(6,100)C(i,1),C(i,2),C(i,3),C(i,4),C(i,5),C(i,6),C(i,7),
$C(i,8),C(i,9),C(i,10)
100 format('',10f10.4)
enddo
endif
*
return
end
|
