Fortran中求解主对角线含0的线性系统的最优算法探讨

我一直用简单的高斯消元法求解线性系统，但当系数矩阵主对角线上有0时，会出现除以0的错误。用Gauss-Jordan变换或Gauss-Seidel迭代法也会遇到同样的问题。作为编程新手，我对线性代数知识也比较生疏，想知道解决这个问题的最佳方案是什么？

我目前使用的高斯消元代码如下：

SUBROUTINE solve_lin_sys(A, c, x, n)
! =====================================================
! Uses gauss elimination and backwards substitution
! to reduce a linear system and solve it.
! Problem (0-div) can arise if there are 0s on main diag.
! =====================================================
  IMPLICIT NONE
  INTEGER:: i, j, k
  REAL*8::fakt, summ
  INTEGER, INTENT(in):: n
  REAL*8, INTENT(inout):: A(n,n)
  REAL*8, INTENT(inout):: c(n)
  REAL*8, INTENT(out):: x(n)

  DO i = 1, n-1                   ! pick the var to eliminate
    DO j = i+1, n                 ! pick the row where to eliminate
      fakt = A(j,i) / A(i,i)      ! elimination factor
      DO k = 1, n                 ! eliminate
        A(j,k) = A(j,k) - A(i,k)*fakt
      END DO
      c(j)=c(j)-c(i)*fakt         ! iterate on known terms
    END DO
  END DO

   ! Actual solving:
    x(n) = c(n) / A(n,n)          ! last variable being solved
    DO i = n-1, 1, -1
      summ = 0.d0
      DO j = i+1, n
        summ = summ + A(i,j)*x(j)
      END DO
      x(i) = (c(i) - summ) / A(i,i)
    END DO

END SUBROUTINE solve_lin_sys

核心解决方案：加入部分选主元（Partial Pivoting） #

在每一步消元前，找到当前列中从第i行到第n行里绝对值最大的元素，把该行和当前的第i行交换（同时交换右端项c的对应行）。这样就能保证消元时的除数（主元）不为0，还能提升计算的数值稳定性。

下面是修改后的带部分选主元的Fortran代码：

SUBROUTINE solve_lin_sys_pivot(A, c, x, n)
! =====================================================
! 带部分选主元的高斯消元法，解决主对角线为0的问题
! =====================================================
  IMPLICIT NONE
  INTEGER:: i, j, k, pivot_row
  REAL*8::fakt, summ, max_val
  INTEGER, INTENT(in):: n
  REAL*8, INTENT(inout):: A(n,n)
  REAL*8, INTENT(inout):: c(n)
  REAL*8, INTENT(out):: x(n)
  REAL*8, DIMENSION(n):: temp_row  ! 临时存储行数据

  DO i = 1, n-1
    ! 1. 寻找当前列的主元（绝对值最大的元素所在行）
    max_val = ABS(A(i,i))
    pivot_row = i
    DO k = i+1, n
      IF (ABS(A(k,i)) > max_val) THEN
        max_val = ABS(A(k,i))
        pivot_row = k
      END IF
    END DO

    ! 2. 如果主元不在当前行，交换行
    IF (pivot_row /= i) THEN
      ! 交换系数矩阵的行
      temp_row = A(i,:)
      A(i,:) = A(pivot_row,:)
      A(pivot_row,:) = temp_row
      ! 交换右端项c的对应元素
      CALL swap(c(i), c(pivot_row))
    END IF

    ! 3. 常规高斯消元步骤
    DO j = i+1, n
      fakt = A(j,i) / A(i,i)
      DO k = i, n  ! 优化：从i开始，前面的元素已经是0，无需计算
        A(j,k) = A(j,k) - A(i,k)*fakt
      END DO
      c(j) = c(j) - c(i)*fakt
    END DO
  END DO

  ! 回代求解
  x(n) = c(n) / A(n,n)
  DO i = n-1, 1, -1
    summ = 0.d0
    DO j = i+1, n
      summ = summ + A(i,j)*x(j)
    END DO
    x(i) = (c(i) - summ) / A(i,i)
  END DO

END SUBROUTINE solve_lin_sys_pivot

! 辅助交换函数
SUBROUTINE swap(a, b)
  IMPLICIT NONE
  REAL*8, INTENT(inout):: a, b
  REAL*8:: temp
  temp = a
  a = b
  b = temp
END SUBROUTINE swap

关键说明： #

• 部分选主元：不仅解决了主对角线为0的问题，还能减少舍入误差，让计算结果更准确。
• 代码优化：消元循环中k从i开始，因为前i-1列的元素已经被消为0，无需重复计算。
• 交换行时要同时处理系数矩阵A和右端项c，保证方程的等价性。

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