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