Fortran90克莱姆法则求解线性系统行列式计算异常排查

问题现象 #

编写大学课程作业的线性系统求解程序时，基于拉普拉斯第一定理实现了计算方阵行列式的递归函数，所有单独测试用例运行结果均正确。后续在程序主体中通过do循环实现克莱姆法则求解线性系统，需要依次计算n个将第i列替换为方程组常数项的矩阵的行列式，此时自定义行列式函数出现稳定复现的特定错误：

• 以3阶矩阵测试时，替换第一列的矩阵行列式计算结果正确
• 第二、第三个替换矩阵的行列式计算结果为理论正确值与前序计算矩阵行列式的和
  已在每次循环迭代时将存储行列式结果的变量cr_det置零，且对应代码段无累加逻辑，因对Fortran语言及编程基础掌握不熟练，始终未定位到异常原因。

注：代码中意大利语注释可忽略

问题代码 #

module functions
implicit none

contains

recursive function determinant (m, mat2) result (det)
    implicit none
    integer, intent(in)                :: m
    real,  intent(in), dimension(m, m) :: mat2
    integer                            :: x, sgn = -1
    real                               :: submat(m-1, m-1), det
    
    if (m == 1) then
        det = mat2(1, 1)
    else if(m == 0) then
        det = 0
    else if(m==2) then
        det = mat2(1,1)*mat2(2,2) - mat2(1,2)*mat2(2,1)
    else 
        do x =1, m
            submat(1:x-1, 1:m-1) = mat2(1:x-1, 2:m)
            submat(x:m-1, 1:m-1) = mat2(x+1:m, 2:m)
            det = det + (sgn**(x+1))*mat2(x, 1) * determinant(m-1, submat)
        end do
    end if

end function determinant
end module functions


program gaia_leita
use functions
implicit none
real(KIND=16)     :: start, finish
integer           :: n, i, j, h
real, allocatable :: matrix(:,:), sol_vec(:), known_terms(:)
real, allocatable :: cramer_matrix(:,:)
real              :: deter, cr_det = 0


write(*,*) "Questo programma calcola il determinante di una matrice M NxN"
write(*,*) "Inserire la dimensione della matrice:"
read(*,*) n
allocate(matrix(n, n))
write(*,*) "Inserire le entrate della matrice M riga per riga:"
do i = 1, n
    do j = 1, n
       read(*, *) matrix(i, j)
    end do
end do

call cpu_time(start)

deter = determinant(n, matrix)

write(*, *) "Il determinante della matrice M è: |M| = ", deter

write(*,*) "Per risolvere il sistema lineare M*v = b inserire i valori dei termini noti b1...bn:"

allocate(known_terms(n))
allocate(sol_vec(n))
allocate(cramer_matrix(n, n))


do i = 1, n
    read(*, *) known_terms(i)
end do

if(deter == 0) then
    write(*,*) "Il sistema non soddisfa le ipotesi del metodo di Cramer"
else
    cramer_matrix = matrix
    do i = 1, n
        cr_det = 0 
        cramer_matrix = matrix
        cramer_matrix(:, i) = known_terms(:)

        do j = 1, n
            do h=1, n
                write(*,*) cramer_matrix(j, h)
            end do
        end do
        
        cr_det = determinant(n, cramer_matrix)

        write(*, *) cr_det
        write(*,*)
        sol_vec(i) = cr_det/ deter
        write(*,*)
    end do
end if

Write(*,*) "La soluzione del sistema è il vettore x = ( "
do i = 1, n
    write(*,*) sol_vec(i), " "
end do
write(*,*) ")"

deallocate(matrix)
deallocate(known_terms)
deallocate(sol_vec)
deallocate(cramer_matrix)

call cpu_time(finish)
write(*,*) "Time elapsed:", finish-start, "seconds"
read(*,*)
end program gaia_leita

测试运行输出 #

Questo programma calcola il determinante di una matrice M NxN
 Inserire la dimensione della matrice:
3
 Inserire le entrate della matrice M riga per riga:
2
3
6
-5
-3
5
1
0
2
 Il determinante della matrice M è: |M| =    51.0000000
 Per risolvere il sistema lineare M*v = b inserire i valori dei termini noti b1...bn:
3
2
4
   3.00000000
   3.00000000
   6.00000000
   2.00000000
  -3.00000000
   5.00000000
   4.00000000
   0.00000000
   2.00000000
   102.000000


   2.00000000
   3.00000000
   6.00000000
  -5.00000000
   2.00000000
   5.00000000
   1.00000000
   4.00000000
   2.00000000
  -17.0000000


   2.00000000
   3.00000000
   3.00000000
  -5.00000000
  -3.00000000
   2.00000000
   1.00000000
   0.00000000
   4.00000000
   34.0000000


 La soluzione del sistema è il vettore x = (
   2.00000000
 -0.333333343
  0.666666687
 )

测试现象说明：每个替换矩阵的元素打印完成后，下方输出的数值为程序计算的行列式值，仅第一个替换矩阵的行列式结果正确，其余结果均为理论正确值与前序行列式计算结果的和。

根因分析 #

异常由两个Fortran语法细节问题共同导致：

1. 核心错误：返回值未初始化
   当矩阵阶数m≥3进入拉普拉斯展开的递归分支时，函数返回值det没有被初始化为0就直接执行det = det + ...的累加操作。Fortran中普通局部变量不会自动初始化，会直接复用栈空间中残留的上一次函数调用留下的数值，导致每次累加都会叠加上前序行列式计算的残留值，和观察到的错误现象完全吻合。
2. 隐患写法：带初值的局部变量默认save属性
   代码中integer :: x, sgn = -1的写法，在Fortran语法规则下，声明同时赋初值的局部变量默认带有save属性，变量值会在多次函数调用间保留，虽然当前场景下符号计算结果未出错，但属于容易引发其他异常的不规范写法。

修复方法 #

修改determinant函数的递归分支，在循环开始前将det初始化为0，同时优化符号计算逻辑避免save属性隐患：

else 
        det = 0.0  ! 补充累加初始值，修复核心bug
        do x =1, m 
            submat(1:x-1, 1:m-1) = mat2(1:x-1, 2:m)
            submat(x:m-1, 1:m-1) = mat2(x+1:m, 2:m)
            ! 直接计算代数余子式符号，移除带save属性的sgn变量隐患
            det = det + ((-1)**(x+1))*mat2(x, 1) * determinant(m-1, submat)
        end do
    end if

修复后重新运行测试，三个替换矩阵的行列式将返回正确值102、-119、51，对应线性系统解为x=(2, -7/3≈-2.333, 1)，计算结果完全符合数学推导。

内容的提问来源于stack exchange，提问作者Guglielmo