如何像MATLAB那样调用UMFPACK

Question

问题

我想解一个一般的线性方程组a*x=b。m乘m矩阵是稀疏的、实的、正方形的、条件较差的、非对称的，但它是奇异的（秩（A）=m-1），因为x仅在加法常数之前是已知的。

我可以通过在三个向量（i、j和v）中指定其非零项来创建矩阵A，使得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;

Numerics

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

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

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

• MATLAB通过UMFPACK V5.4.0对平方、m乘m的原始方程使用LU分解
• MATLAB通过SuiteParseQR1.1.0对m-by-（m+1）修正方程进行QR分解

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

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

我的问题

既然我可以在MATLAB中解决这个问题，我希望在Fortran中解决。不幸的是，MATLAB的mldivide()命令是一个黑匣子，因为我看不出它是如何设置或调用SuiteParse例程的。

假设我在Fortran（90+）中有三个稀疏向量，如下所示，我如何使用SuiteParse解决问题

或者，有人知道有任何F90包装程序可以与UMFPACK接口，使其更容易吗？

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

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在幕后做了什么来设置对SuiteParse的调用？（似乎没有记录……）
• 从Fortran中调用SuiteParse需要做些什么？（如果有足够的时间，我可能会从这个演示中了解到最多，但奇怪的是，它多次调用例程…）

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

Answer 1

这里有两件事。首先，对于在fortran/C中调用UMFPACK的示例，请检查我的一个项目中的源代码，例如

https://github.com/fangq/blit/blob/master/src/blit_ilupcond.f90#L156-L162

这里我使用umf4_f77wrapper.c单元从fortran90调用umfpack

https://github.com/fangq/blit/blob/master/src/umf4_f77wrapper.c

其次，关于反斜杠算子，它本质上是解决Moore Penrose伪反演问题，以获得最小二乘解，见

https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse

基本上，任何矩阵方程A*x=b（包括未确定大小（A，1）<大小（A，2）-或超定-大小（A、1）>大小（A，2）-问题），解决x=Ab与解决相同

CCD_ 16＝>x=inv(A'A)*A'b

然而，你不是通过反演来求解右边的方程，而是应用线性求解器来求解左边的方程(A'A)x=A'b

当原始方程不确定时，也可以应用Sherman–Morrison–Woodbury公式来构造一个等价但更紧凑的解，即x=A'*inv(A'A)*b

https://en.wikipedia.org/wiki/Woodbury_matrix_identity