16

问题

我希望求解线性方程组 A*x=b 的一般系统。m×m 矩阵是稀疏的、实数的、方形的、条件较差的、非对称的，但它是奇异的（rank(A)==m-1），因为x只知道一个加性常数。

i我可以通过在三个向量（ 、j和）中指定其非零条目来创建矩阵 A，v使得A(i(k),j(k)) = v(k)：

A = sparse( i, j, v, m, m );

原始方程

我可以解决这个原始方程如下：

x = A \ b;

如果我想要一个唯一的解决方案，我可以在计算非唯一解决方案后施加一个约束（例如，x(4) == 3.14159）：

x = x - x(4) + 3.14159;

修正方程

我可以通过添加一个额外的唯一性约束来创建一个新的全秩矩阵C，如下所示：

% Add the constraint x(4) == 3.14159
extraRow = zeros(1,m);
extraRow(4) = 1.0;
C = [A; extraRow];    % Add to the matrix A
d = [b; 3.14159];     % Add to the RHS vector, b

% Solve C*y = d for y
y = C \ d;

数字

我知道，当我通过x = A \ b或求解这些方程时y = C \ b，MATLAB 将 解释\为mldivide()命令（链接），该命令对矩阵运行一些测试并找出求解例程的最佳算法（有关详细信息，请参阅链接）。

通过设置 MATLAB 的稀疏矩阵求解参数，这些细节在运行时变得冗长spparms('spumoni',2)

当我计算x和/或y时，我注意到以下内容：

  • MATLAB 通过 UMFPACK V5.4.0 将 LU 分解用于平方 m×m 原始方程。
  • MATLAB 通过 SuiteSparseQR 1.1.0 对 m×(m+1) 修改方程使用 QR 分解。

UMFPACK 和 SuiteSparseQR 都是 SuiteSparse 软件包（链接）的一部分。

（有点出乎意料的是，求解修改后的方程会产生比原始方程更多的错误。虽然很重要，但这个错误仍然处于可接受的低阈值。）

我的问题

现在我可以在 MATLAB 中解决这个问题，我希望在 Fortran 中这样做。不幸的是，MATLAB 的mldivide()命令是一个黑匣子，我看不到它是如何设置或调用SuiteSparse例程的。

鉴于我在 Fortran (90+) 中有三个稀疏向量，如下所示，如何使用SuiteSparse解决问题？

或者，是否有人知道与 UMFPACK 接口的任何 F90 包装例程以使这更容易？

如果有人可以帮助解决这个问题，我会非常高兴——无论是原始方程还是修正方程。（如果你帮助一个，我可能会得到另一个。）

subroutine solveSparseMatrixEqnViaSuiteSparse( m, n, nnz, i, j, v, x )
    implicit none

    integer, intent(in)                   :: m     ! sparse matrix rows
    integer, intent(in)                   :: n     ! sparse matrix columns
    integer, intent(in)                   :: nnz   ! number of nonzero entries
    integer, dimension(1:nnz), intent(in) :: i     ! row indices of nonzero entries
    integer, dimension(1:nnz), intent(in) :: j     ! column indices of nonzero entries
    real*8,  dimension(1:nnz), intent(in) :: v     ! values of nonzero entries
    real*8,  dimension(1:n), intent(out)  :: x     ! solution vector

    ! I have no idea what to do next! 

end function solveSparseMatrixEqnViaSuiteSparse

让我感到困惑的是以下几点：

  • MATLAB 在幕后做了什么来设置对 SuiteSparse 的调用？（它似乎没有记录在案......）
  • 从 Fortran 中调用 SuiteSparse 需要做什么？（如果有足够的时间，我可能可以从这个演示中弄清楚大部分，但奇怪的是它多次调用例程......）

注意：尽管我在提出这个问题时考虑到了一个特定的问题（谁不是！？），但我相信这足够普遍，对其他人也有用。

4

0 回答 0