mimp.f90

求解对角线元素很小的矩阵的逆矩阵,避免一般求逆的溢出问题,用Fortran编写

program main
real::m(4,4)
integer::row
real::returnvalue(4,4)
m(1,1)=1.0;m(1,2)=3.0;m(1,3)=0.0;m(1,4)=1.0
m(2,1)=2.0;m(2,2)=2.0;m(2,3)=1.0;m(2,4)=3.0
m(3,1)=3.0;m(3,2)=4.0;m(3,3)=5.0;m(3,4)=6.0
m(4,1)=0;m(4,2)=7.0;m(4,3)=3.0;m(4,4)=1.0
row=4
call mimp(m,returnvalue,row)
do i=1,4
  print*,(returnvalue(i,j),j=1,4)
enddo
end program main
subroutine mimp(m,returnvalue,row)
!这个程序是用来求解矩阵的逆矩阵,而m是要求逆矩阵的矩阵,也就是原来的矩阵
!而returnvalue是所求出的逆矩阵,row是矩阵的维数.
implicit none
integer::row
real::m(row,row),returnvalue(row,row)
real::temp(row)
logical::nocompare
integer::i,j,n,k,r(row),c(row)
real::pivot
!---------------------
print*,3==5
do i=1,row
    returnvalue(i,i)=1
enddo
do i=1,row
    pivot=0
    nocompare=.False.
    do j=1,row
       do n=1,i-1
         nocompare=(j==c(n))
         if(nocompare)exit 
       enddo
       if(.not.nocompare)then
           do k=1,row
              do n=1,i-1
                nocompare=(k==r(n))
                if(nocompare) exit
              enddo
            if(.not.nocompare)then
               if(abs(m(k,j))>=Pivot)then
                pivot=abs(m(k,j))
                r(i)=k
			    c(i)=j
               endif
            endif
           enddo
        endif
     enddo
     pivot=m(r(i),c(i))
     if(abs(Pivot-1)>0)then
        do k=1,row
          m(r(i),k)=m(r(i),k)/pivot
          returnValue(r(i),k)=returnvalue(r(i),k)/pivot
        enddo
     endif
     do j=1,row
        pivot=m(j, c(i))
        if(abs(j-r(i))>0.and. abs(pivot)>0) then
           do k=1,row
              m(j,k)=m(j,k)-m(r(i),k)*pivot
              returnValue(j,k)=returnvalue(j,k)-returnvalue(r(i),k)*pivot
           enddo
        endif
     enddo
enddo
do i=1,row-1
  if(abs(r(i)-c(i))>0) then
       do j=1,row
         temp(j)=returnvalue(r(i),j)
         returnvalue(r(i),j)=returnvalue(c(i),j)
         returnvalue(c(i),j)=temp(j)
       enddo
	   do j=i+1,row
         if(r(j)==c(i))then 
			r(j)=r(i)
			exit 
	     endif
       enddo
       r(i)=c(i)
   endif
enddo
end subroutine mimp