15
$\begingroup$

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
10
$\begingroup$

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.)

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – NoahR
    Commented Feb 19, 2013 at 18:04
  • $\begingroup$ 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. $\endgroup$
    – NoahR
    Commented Feb 20, 2013 at 2:33
  • $\begingroup$ 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 ); $\endgroup$
    – NoahR
    Commented Feb 20, 2013 at 6:51
2
$\begingroup$

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$.

Improve this answer
$\endgroup$
1
  • 3
    $\begingroup$ If $A$ is rectangular (not square), then so is $R$, and the expression $R^{-1}$ is undefined. $\endgroup$
    – Geoff Oxberry
    Commented Feb 19, 2013 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.