使用矩阵求逆程序出现NaN

Escaliborc

程序出现NaN

[Fortran] 纯文本查看 复制代码

module inv_mat
!----------------------------------------module coment
!  Version     :  V1.0    
!  Coded by    :  syz 
!  Date        :  2010-4-8
!-----------------------------------------------------
!  Description : 计算逆矩阵
!    
!-----------------------------------------------------
!  Contains    :
!      1.   driver  驱动函数
!      2.   solve   方法函数 
!-----------------------------------------------------

contains

subroutine drive(N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  
!-----------------------------------------------------
!  Purpose   :  驱动程序
!               用以读取数据，计算结果并把结果保存于文件中
!
!----------------------------------------------------
!  Post Script :
!       1.     fin   输入文件
!       2.     fout  输出文件
!----------------------------------------------------
implicit real*8(a-z)

integer::N,i,j
real*8::A(N,N),invA(N,N)

read(11,*)((A(i,j),j=1,N),i=1,N)

write(12,101)
101 format(/,T12,'逆矩阵为',/)

call solve(A,invA,N)


write(12,102)((invA(i,j),j=1,n),i=1,n)

!变量输出格式只针对 IVF编译器，在CVF中不支持
102 format(<N>(<N>F20.5/))

end subroutine drive




subroutine solve(A,invA,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  
!-----------------------------------------------------
!  Purpose   :  计算逆矩阵
!    
!-----------------------------------------------------
!  Input  parameters  :
!       1.
!       2. 
!  Output parameters  :
!       1.
!       2.
!  Common parameters  :
!
!----------------------------------------------------
!  Post Script :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)

integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)

E=0

!设置E为单位矩阵
do i=1,n
   E(i,i)=1
end do

call mateq(A,E,invA,N,N)


end subroutine solve


subroutine mateq(A,B,X,N,M)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-4-8
!-----------------------------------------------------
!  Purpose   :  高斯列主元消去法计算矩阵方程组
!                 AX=B
!-----------------------------------------------------
!  Input  parameters  :
!       1.   A(N,N)系数矩阵
!       2.   B(N,M)右矩阵
!       3.   N方程维数
!       4.   M右矩阵的列数
!  Output parameters  :
!       1.  X 方程的解矩阵
!       2.
!  Common parameters  :
!
!----------------------------------------------------
implicit real*8(a-z)

integer::i,k,N,M
integer::id_max  !主元素标号

real*8::A(N,N),B(N,M),X(N,M)

real*8::Aup(N,N),Bup(N,M)

!Ab为增广矩阵 [AB]
real*8::AB(N,N+M)

real*8::vtemp1(N+M),vtemp2(N+M)
real*8::vtmp(N),xtmp(N)

AB(1:N,1:N)=A

AB(:,N+1:N+M)=B


!##########################################################
!  这段是列主元消去法的核心部分
do k=1,N-1

    elmax=dabs(Ab(k,k))
    id_max=k

    !这段为查找主元素        
    !这段程序的主要目的不是为了赋值最大元素给elmax，而是为了找出最大元素对应的标号

        
        do i=k+1,n
      if (dabs(Ab(i,k))>elmax) then
         elmax=Ab(i,k)

         id_max=i
      end if          
    end do


!至此，已经完成查找最大元素，查找完成以后与 第k行交换
!交换两行元素，其他不变
    vtemp1=Ab(k,:)
    vtemp2=Ab(id_max,:)


    Ab(k,:)=vtemp2
    Ab(id_max,:)=vtemp1   
!
!以上一大段是为交换两行元素，交换完成以后即按照消元法进行
!#########################################################

   do i=k+1,N

     temp=Ab(i,k)/Ab(k,k)

     Ab(i,:)=Ab(i,:)-temp*Ab(k,:)

   end do

end do

!-----------------------------
! 经过上一步，Ab已经化为如下形式的矩阵
!            | *  *  *  *  #  #|
!     [A b]= | 0  *  *  *  #  #|
!            | 0  0  *  *  #  #|
!            | 0  0  0  *  #  #|
!
Aup(:,:)=AB(1:N,1:N)




do i=1,m
!调用用上三角方程组的回带方法
vtmp=AB(:,N+i)
call uptri(Aup,vtmp,xtmp,n)
!把计算结果赋值给X
X(:,i)=xtmp
end do

end subroutine mateq


subroutine uptri(A,b,x,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-4-8
!-----------------------------------------------------
!  Purpose   :  上三角方程组的回带方法
!                 Ax=b
!-----------------------------------------------------
!  Input  parameters  :
!       1.   A(N,N)系数矩阵
!       2.   b(N)右向量
!       3.   N方程维数
!  Output parameters  :
!       1.  x  方程的根
!       2.
!  Common parameters  :
!
!----------------------------------------------------

implicit real*8(a-z)

integer::i,j,N

real*8::A(N,N),b(N),x(N)

x(N)=b(N)/A(N,N)

!回带部分
do i=n-1,1,-1

    x(i)=b(i)
   do j=i+1,N
    x(i)=x(i)-a(i,j)*x(j)
   end do
    x(i)=x(i)/A(i,i)

end do

end subroutine uptri

end module inv_mat




program main
!--------------------------------------program comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2010-4-8
!-----------------------------------------------------
!  Purpose   :  逆矩阵的计算
!    
!-----------------------------------------------------
!  In put data  files :
!       1.  fin.txt  输入方程系数
!       2.
!  Output data files  :
!       1. fout.txt  计算结果
!       2.
!-----------------------------------------------------
!  Post Script :
!       1.    需要准备输入数据
!       2.    由于驱动函数调用方法函数
!-----------------------------------------------------
use inv_mat

integer::N

!N表示矩阵的维数，由文件读入

open(unit=11,file='fin.txt')
open(unit=12,file='fout.txt')

read(11,*)
!读说明文字

!读入方程维数系数
read(11,*)N


call drive(N)


end program main

Escaliborc

            输入矩阵
     6        
            1.095445            0.000000            0.000000            0.000000            0.000000            0.000000
           45.643546          176.776692            0.000000            0.000000            0.000000            0.000000
           -1.095445            0.000000            1.095445            0.000000            0.000000            0.000000
           45.643546          -58.925564           45.643545          242.956325            0.000000            0.000000
            0.000000            0.000000           -1.095445            0.000000            0.000048            0.000000
            0.000000            0.000000           45.643545          -42.874645            0.000000          171.498582

这是输入文档

vvt

mateq函数中，第147，148 行。
记住，列主元是找绝对值大的，而不是数值大的。

本例中，你遇到了 -234，-5 和 0。然后你就找到 0 的那列了。然后算比例因子 temp=Ab(i,k)/Ab(k,k)，得，除数为 0 ，结果为 NaN

147      if (dabs(Ab(i,k))>elmax) then
148         elmax=Ab(i,k)
改为
147      if (dabs(Ab(i,k))>elmax) then
148         elmax=abs(Ab(i,k))