15

I am porting an existing code from MATLAB to C++ and have a linear system to solve $xA=b$ (rather than the more typical form $Ax=b$)

The matrix $A$ is dense, and of general form, but is no larger than 1000x1000. So in MATLAB, the solution is found by the mrdivide(b,A) function or the forward-slash notation x = b/A;

How should I solve this in my C++ code using BLAS and LAPACK routines?

I'm familiar with the LAPACK routine DGESV which solves $Ax=b$ for $x$.

So, one thought I had is to do some manipulations using matrix transpose identities:

$(xA)^T=b^T$

$A^T x^T = b^T$

$x^T = (A^T)^{-1} b^T$

Then solve the final form using DGESV operating on the transposed $A^T$. (so cost to transpose $A$ and cost to solve system)

Is there an approach more efficient or otherwise better?

I am working with matrix and vector classes as well as BLAS implementation from the BOOST uBLAS library as well as bindings to the LAPACK library routines. I've been using this setup successfully for other operations and am hoping to find a solution limited to these libraries.

Also, I should note I am only performing this type of operation a few times during code setup, so performance is not a critical concern.

Maybe this MATLAB documentation on mrdivide is helpful for others.

Christian Clason
  • 12,301
  • 3
  • 48
  • 68
NoahR
  • 253
  • 2
  • 4

2 Answers2

10

Trivial answer for square $A$: use dgesvx which solves also for $A^T x = b$ when TRANS = 'T'.

Please note that with BLAS or LAPACK you hardly have to transpose (swapping elements in memory) a matrix: most of the subroutines have a TRANS argument to accommodate for operation on the transpose matrix or on a matrix stored with a different memory layout. (Transposition is equivalent to changing Fortran-contiguous memory layout to a C-contiguos one and viceversa.)

Stefano M
  • 3,839
  • 16
  • 24
  • Thanks for the answer and explanation! I have done very little work with LAPACK and now I know to look for the TRANS option. I'm having trouble getting the TRANS argument to work through boost::numeric::bindings::lapack::gesvx(), but this not part of my question here. If I do have success, I'll come back with a note on how to do it. – NoahR Feb 19 '13 at 18:04
  • I have a working solution using gesvx(), though not without some stumbling along the way. When the TRANS argument is 'T', the LAPACK documentation says gesvx solves $A^T X = B$, but really it solves $A^T X^T = B^T$ because the form of the input arguments $X$ and $B$ are expected to be in their non-transposed form. So, argument $A$ is transposed, while $X$ and $B$ are not. Great, that's more convenient.

    If anyone else stumbles upon this trying to use boost numeric bindings, I'll say I haven't been able to get the transpose interface used in this soln. to work through the bindings.

     – NoahR Feb 20 '13 at 02:33
  • Ok, I found the trick to using the gesvx transpose form through boost::numeric::bindings. Wrap the transposed matrix $A^T$ in the trans() function. This identifies the parameter as transpose type:

    boost::numeric::bindings::lapack::gesvx( FACT, boost::numeric::bindings::trans(Atransposed), af, ipiv, equed, r, c, b, x, rcond, ferr, berr );

     – NoahR Feb 20 '13 at 06:51
1

You can compute the pseudo inverse of $A$, by using say QR decomposition of SVD, which are both included in LAPACK.

\begin{equation} xA=b\\ xQR=b\\ x=bR^{-1}Q^T \end{equation}

This will work for any rectangular $A$.

Gil
  • 382
  • 1
  • 10