title: How to call UMFPACK as MATLAB does?
body: <h1>The problem</h1>
  
  <p>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 <code>x</code> is known only up to an additive constant. </p>
  
  <p>I can create the matrix A by specifying its nonzero entries in three vectors (<code>i</code>, <code>j</code>, and <code>v</code>) such that <code>A(i(k),j(k)) = v(k)</code>:</p>
  
  <pre><code>A = sparse( i, j, v, m, m );
  </code></pre>
  
  <h3>Original Equation</h3>
  
  <p>I can solve this original equation as follows:</p>
  
  <pre><code>x = A \ b;
  </code></pre>
  
  <p>If I want a unique solution, I can impose a constraint (say, x(4) == 3.14159) after calculating the non-unique solution:</p>
  
  <pre><code>x = x - x(4) + 3.14159;
  </code></pre>
  
  <h3>Modified Equation</h3>
  
  <p>I can make a new, full-rank matrix <code>C</code> by adding an additional uniqueness constraint, as follows: </p>
  
  <pre><code>% 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;
  </code></pre>
  
  <h1>Numerics</h1>
  
  <p>I understand that when I solve these equations via either <code>x = A \ b</code> or <code>y = C \ b</code>, MATLAB interprets the <code>\</code> as the <code>mldivide()</code> command (<a href="http://www.mathworks.com/help/matlab/ref/mldivide.html" rel="noreferrer">link</a>), which runs some tests on the matrix and figures out an optimal algorithm to solve the routine (see link for details). </p>
  
  <p>These details are made verbose at runtime by setting MATLAB's sparse matrix solving parameters via <code>spparms('spumoni',2)</code></p>
  
  <p>When I compute <code>x</code> and/or <code>y</code>, I notice the following: </p>
  
  <ul>
  <li>MATLAB uses LU decomposition via UMFPACK V5.4.0 for the square, m-by-m original equation. </li>
  <li>MATLAB uses QR decomposition via SuiteSparseQR 1.1.0 for the m-by-(m+1) modified equation. </li>
  </ul>
  
  <p>Both UMFPACK and SuiteSparseQR are part of the SuiteSparse software package (<a href="http://faculty.cse.tamu.edu/davis/suitesparse.html" rel="noreferrer">link</a>). </p>
  
  <p>(Somewhat unexpectedly, solving the modified equation gives more error than the original equation. While significant, this error is still at an acceptably-low threshold.)</p>
  
  <h1>My Problem</h1>
  
  <p>Now that I can solve this problem in MATLAB, I wish to do so in Fortran. Unfortunately, MATLAB's <code>mldivide()</code> command is a black-box in that I cannot see how it sets up or calls the <em>SuiteSparse</em> routines. </p>
  
  <p><strong>Given that I have the three sparse vectors in Fortran (90+) as shown below, how can I solve the problem using <em>SuiteSparse</em>?</strong> </p>
  
  <p>Alternatively, does anyone know of any F90 wrapper routines that interface to UMFPACK to make this easier? </p>
  
  <p>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.) </p>
  
  <pre><code>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
  </code></pre>
  
  <p>What confuse me are the following:</p>
  
  <ul>
  <li>What is MATLAB doing behind the scenes to set up the call to SuiteSparse? (It does not appear to be documented....)</li>
  <li>What needs to be done to call SuiteSparse from within Fortran? (I could probably figure most out from <a href="https://github.com/PetterS/SuiteSparse/blob/master/UMFPACK/Demo/umf4zhb.f" rel="noreferrer">this demo</a>, given enough time, but it is odd that it calls routines several times....)</li>
  </ul>
  
  <p>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. </p>