Fortran程序求行列式

Fortran程序求行列式

先吐槽一波，Fortran真的不好用，很多简单的函数功能都要下载相应的函数库进行，没有或不想下载了只能自己进行编写！！！真的费时费力
当时程序主要考虑两种写法代数余子式法和等价转化法

1. 代数余子式法
   以三阶方阵为例：$$T_{3\times 3}=\left|\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right|$$以第一行展开$$T_{3\times 3}=a_{11}\left|\begin{array}{cc}{a}_{22}& a_{23}\\ a_{32}& a_{33}\end{array}\right|-a_{12}\left|\begin{array}{cc}{a}_{21}& a_{23}\\ a_{31}& a_{33}\end{array}\right|+a_{13}\left|\begin{array}{cc}{a}_{21}& a_{22}\\ a_{31}& a_{32}\end{array}\right|$$
   可以想象当行列式阶数更高时比如4阶，其计算机转换的迭代次数为（4-2）=2次，需要拆解共4*3个可供计算器求解的基础的2阶行列式。同理n(N>2)阶行列式需要拆解为共n!/2个可供计算器求解的基础的2阶行列式，所以计算量较大。

2. 等价转化法
   a.行列式的某一行（列）的各元素乘同一数然后加到另一行（列）对应的元素上去，行列式不变
   b. 行列式中某一行（列）的所有元素的公因子可以提到行列式记号的外面
   转化法的核心思想是将行列式转化成上三角行列式
   例：
   $$D=\left|\begin{array}{cccc}1& 2& 3& 7\\ 7& 9& 9& 6\\ 4& 9& 5& 6\\ 1& 2& 4& 5\end{array}\right|$$$$\Downarrow$$$$\left|\begin{array}{cccc}1& 2& 3& 7\\ 0& -5& -12& -43\\ 0& 1& -7& -22\\ 0& 0& 1& -2\end{array}\right|\Rightarrow \left|\begin{array}{cccc}1& 2& 3& 7\\ 0& 0& -47& -153\\ 0& 1& -7& -22\\ 0& 0& 1& -2\end{array}\right|\Rightarrow \left|\begin{array}{cccc}1& 2& 3& 7\\ 0& 0& 0& -247\\ 0& 1& -7& -22\\ 0& 0& 1& -2\end{array}\right|$$$$\Downarrow$$$$\left|\begin{array}{cccc}1& 2& 3& 7\\ 0& 1& -7& -22\\ 0& 0& 1& -2\\ 0& 0& 0& -247\end{array}\right|=-247$$
   可以看到其过程较为精简，故编程也采用该方法
   首先，我们采用matlab进行上述编程验证结果代码如下：

a=[1 2 3 7;7 9 9 6;4 9 5 6;1 2 4 5]
b=det(a)

结果为

a =

     1     2     3     7
     7     9     9     6
     4     9     5     6
     1     2     4     5
b =

   -247.0000

结果正确。
而后是利用fortran采用等价转化法编程，其主函数为

program main
    implicit none
    integer::N
    real*8::A(4,4),d
    N=4
    data A /1,7,4,1,2,9,9,2,3,9,5,4,7,6,6,5/
    call determinant(A,N,d)
    write(*,*) d
end program main

这里简单，fortran数组的批量赋值遵从的原则是先列后行，所以这边的话可能有点奇怪。
Fortran求解行列式的子函数为

subroutine determinant(A,N,d)
    !----------------------------------------------------------------------
    !--------输入-------
    !A——方阵A(N,N)
    !N——A的维度
    !--------输出--------
    !d——行列式
    !-----------------------------------------------------------------------
    implicit real*8(a-h,o-z), integer*4 (i-n)
    real*8 A(N,N),A_New(N,N)
    real*8 Paixu(N,1)
    !转化法
    xuhao=0.d0
    A_New=A
    d=1.d0
    do i=1,N
        Paixu(i,1)=0.d0
    enddo
    do j=1,N
        do i=1,N
        if((A_New(i,j)/=0.d0).and.(Paixu(i,1)==0.d0))then
            xuhao=xuhao+1.d0
            Paixu(i,1)=xuhao
            d=d*A_New(i,j)
            do m=1,N
                if(Paixu(m,1)==0.d0)then
                    A_New(m,:)=A_New(m,:)-A_New(i,:)*A_New(m,j)/A_New(i,j)
                endif
            enddo
            exit
        endif
        enddo
    enddo
    
    !write(*,*) Paixu
    !write(*,*)
    
    !逆序数
    Nixu_num=0.d0
    do i=1,N-1
        do j=i+1,N
            if(Paixu(i,1)>Paixu(j,1))then
                Nixu_num=Nixu_num+1.d0
            endif
        enddo
    enddo
    !write(*,*) Nixu_num
    !write(*,*)
    
    !d——行列式
    d=(-1.d0)**Nixu_num*d
    
 end

结果展示

主要思路是先进行等价转化，转化过程中记录缩放得因子，并记录序号，后续根据排序求解逆序数，最后得到结果