17

The problem

I wish to solve a general system of linear equations A*x=b. The m-by-m matrix is sparse, real, square, somewhat poorly-conditioned, non-symmetric, but it is singular (rank(A)==m-1) because x is known only up to an additive constant.

I can create the matrix A by specifying its nonzero entries in three vectors (i, j, and v) such that A(i(k),j(k)) = v(k):

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

Original Equation

I can solve this original equation as follows:

x = A \ b;

If I want a unique solution, I can impose a constraint (say, x(4) == 3.14159) after calculating the non-unique solution:

x = x - x(4) + 3.14159;

Modified Equation

I can make a new, full-rank matrix C by adding an additional uniqueness constraint, as follows: 

% 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

I understand that when I solve these equations via either x = A \ b or y = C \ b, MATLAB interprets the \ as the mldivide() command (link), which runs some tests on the matrix and figures out an optimal algorithm to solve the routine (see link for details).

These details are made verbose at runtime by setting MATLAB's sparse matrix solving parameters via spparms('spumoni',2)

When I compute x and/or y, I notice the following:

  • MATLAB uses LU decomposition via UMFPACK V5.4.0 for the square, m-by-m original equation.
  • MATLAB uses QR decomposition via SuiteSparseQR 1.1.0 for the m-by-(m+1) modified equation.

Both UMFPACK and SuiteSparseQR are part of the SuiteSparse software package (link).

(Somewhat unexpectedly, solving the modified equation gives more error than the original equation. While significant, this error is still at an acceptably-low threshold.)

My Problem

Now that I can solve this problem in MATLAB, I wish to do so in Fortran. Unfortunately, MATLAB's mldivide() command is a black-box in that I cannot see how it sets up or calls the SuiteSparse routines.

Given that I have the three sparse vectors in Fortran (90+) as shown below, how can I solve the problem using SuiteSparse?

Alternatively, does anyone know of any F90 wrapper routines that interface to UMFPACK to make this easier?

I'm more than happy if anyone can help with this--either with the Original Equation or with the Modified Equation. (If you help with one, I could probably get the other.) 

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

What confuse me are the following:

  • What is MATLAB doing behind the scenes to set up the call to SuiteSparse? (It does not appear to be documented....)
  • What needs to be done to call SuiteSparse from within Fortran? (I could probably figure most out from this demo, given enough time, but it is odd that it calls routines several times....)

Note: Although I'm asking this question with a particular problem in mind (who isn't!?), I believe this is general enough to be useful to others as well.

jvriesem
  • 1,859
  • 3
  • 18
  • 40
  • 4
    Maybe you can take a look at what [Octave is doing](http://hg.savannah.gnu.org/hgweb/octave/file/4abcc0969ebd/liboctave/array/dSparse.cc)... – Alexander Vogt Mar 30 '16 at 20:49
  • 1
    Just as a side note: [Using i and j as variables in Matlab](http://stackoverflow.com/questions/14790740/using-i-and-j-as-variables-in-matlab) – Matthias W. Mar 30 '16 at 22:05
  • 3
    @MatthiasW. In recent versions of Matlab it doesn't matter anymore, there are no consequences in terms of speed. The better advice would be *not to use i and j as imaginary units* (but `1i` and `1j`) - as that makes a difference. – Robert Seifert Mar 31 '16 at 06:53

1 Answers1

0

Two things here. First, for examples calling UMFPACK in fortran/C, please checkout the source code in one of my projects, for example,

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

here I am call umfpack from fortran90 using a umf4_f77wrapper.c unit

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

Secondly, regarding the back-slash operator, it is essentially solving the Moore-Penrose pseudo-inversion problem to get a least-square solution, see

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

basically, any matrix equation A*x=b (including under-determined - size(A,1)<size(A,2) - or overdetermined - size(A,1)>size(A,2) - problems), solving x=A\b is the same as solving

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

however, you don't solve the right-equation by taking an inversion, but apply a linear solver solving the left equation (A'A)x=A'b

you can also apply the Sherman–Morrison–Woodbury formula to construct a equivalent, but more compact, solution when the original equation is under-determined, which is x=A'*inv(A'A)*b

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

FangQ
  • 1,444
  • 10
  • 18
  • if A is rank-deficient or ill-posed, you must add regularization term (Tichonov for example) or use truncated SVD or other algorithms to make the left-hand-side (A'A or AA') invertable. – FangQ Jul 02 '23 at 21:10