如何像MATLAB一样调用UMFPACK？

Question

问题

我希望解决一般的线性方程组A*x = b.m-by-m矩阵是稀疏的,实数,正方形,条件稍差,非对称,但它是奇异的(秩(A)== m-1),因为x只知道一个加性常数.

我可以通过在三个矢量(指定其非零项创建矩阵A i,j和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;

NUMERICS

据我所知,当我经由任一解决这些方程x = A \ b或y = C \ b,MATLAB解释\作为mldivide()命令(链路),它运行在基质一些测试,计算出最佳的算法来求解例程(参见链接的详细信息).

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

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

• MATLAB使用通过UMFPACK V5.4.0进行LU分解的方形,m-by-m原始方程.
• MATLAB使用SuiteSparseQR 1.1.0进行QR分解,得到m-by-(m + 1)修正方程.

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需要做些什么？(如果有足够的时间,我可能会从这个演示中找到最多的东西,但奇怪的是它多次调用例程....)

注意:虽然我在考虑这个问题时考虑到了一个特定的问题(谁不是!？),但我相信这也足以对其他人有用.

Answer 1

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

\n

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

\n

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

\n

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

\n

其次，关于反斜杠算子，本质上是求解Moore-Penrose伪反演问题，得到最小二乘解，参见

\n

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

\n

基本上，任何矩阵方程 A*x=b （包括欠定 - size(A,1)<size(A,2) - 或超定 - size(A,1)>size(A,2) - 问题），解决x=A\\b与解决相同

\n

(A\'A) x = A\'b => x=inv(A\'A)*A\'b

\n

但是，您不能通过反演来求解右方程，而是应用线性求解器求解左方程(A\'A)x=A\'b

\n

当原始方程欠定时，您还可以应用 Sherman\xe2\x80\x93Morrison\xe2\x80\x93Woodbury 公式构造等效但更紧凑的解，即x=A\'*inv(A\'A)*b

\n

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

\n