# Programming Example 11

Program report11
implicit none
character(len=50):: matrixfile
integer:: status=0, n
real(8):: r
real(8), dimension(:,:), allocatable:: Matrix, Inv_Matrix, Identity
integer:: i, j, row

write(*,*) "Enter the file name of matrix"
read(*,*) matrixfile
!===== Read file =====
open(1, file=matrixfile)
n=0
do
read(1, *, iostat=status) r
if (status/=0) exit
n=n+1
end do
rewind(1)
allocate(Matrix(n,n))
do i=1, n
read(1, *) Matrix(i, :)
end do
close(1)
!===== Process data =====
allocate(Inv_Matrix(n,n), Identity(n,n))
Inv_Matrix = Minverse(Matrix)
Identity = matmul(Matrix, Inv_Matrix)

!Check
print*, "The matrix before inversion."
do i=1, n
write(*,'(a)',advance="no") "|"
do j=1, n
write(*,'(E12.3)',advance="no") Matrix(i,j)
end do
write(*,*) "|"
end do

!===== Write file =====
open(2, file="invmat.dat")

print*,"The matrix after inversion"
do i=1, n
write(*,'(a)',advance="no") "|"
do j=1, n
write(*,'(E12.3)',advance="no") Inv_Matrix(i,j)
end do
write(*,*) "|"
end do

print*, "Check the identity matrix"
do i=1, n
write(*,'(a)',advance="no") "|"
do j=1, n
write(2,'(E12.3)',advance="no") Inv_Matrix(i,j)
write(*,'(E12.3)',advance="no") Identity(i,j)
end do
write(2,*) " "
write(*,*) "|"
end do
close(2)

Contains
Function Minverse(A)
real(8), dimension(:,:):: A
real(8), dimension(:,:), allocatable:: Minverse, B, V
real(8), dimension(:), allocatable:: temp
real(8), parameter:: TOL=1D-6
real(8):: pivot
integer:: i, j, jmax, k, n,l
n=size(A)**0.5
allocate(Minverse(n,n),B(n,n),V(n,n),temp(n))
B=A
V(:,:)=0.0
do i=1, n
V(i,i)=1.0
end do

do i=1, n
jmax=i
do j=i+1, n
if (abs(B(j,i))>abs(B(jmax,i))) then
jmax=j
temp=B(i,:);B(i,:)=B(jmax,:);B(jmax,:)=temp
temp=V(i,:);V(i,:)=V(jmax,:);V(jmax,:)=temp
end if
end do
pivot=B(i,i)
if (abs(pivot)<TOL) then
print *, "The matrix is singular!"; stop
end if

B(i,:)=B(i,:)/pivot !All pivots on B(i,:) are divide by diagnole in order to set diagnole 1.
V(i,:)=V(i,:)/pivot !Apply the above process for V which leads to Inverse of the matrix.
!This is the wrong code I wrote.
! do j=1,n
!   B(i,j)=B(i,j)/B(i,i)
!   V(i,j)=V(i,j)/B(i,i)
! end do
!If we failed to set pivot=B(i,i), the value is changeable so we can't get correct value.

do k =1,n
if (i .ne. k) then
pivot=B(k,i) !Set pivot B(k,i) which means the number located on the column the above process works.
do l=1,n
B(k,l)=B(k,l)-B(i,l)*pivot !B(k,i) is changed to new row in order to change one column into 0 .
V(k,l)=V(k,l)-V(i,l)*pivot !Apply above process for V.
end do
end if
enddo
end do
!This is the wrong code I wrote.
!  do l=1,n
!   B(k,l)=B(k,l)-B(i,l)*B(k,i)
!   V(k,l)=V(k,l)-V(i,l)*B(k,i)
! end do
!If we failed to set pivot=B(k,i),the same thing as above will happen.
Minverse=V
end Function Minverse

End Program report11