Appropriate Lapack/MKL routines to efficiently compute C = A* inv(B)

I know that, from numerical point of view, computing

Ax = b
B=inv(A), x= B*b


are completely different things, and we should factor the matrix using TRF routine then solve with TRS/TRI routine, or in a combined gesv call.

Anyway is there an efficiently routine to specifically deal with

C = A* inv(B)


BTW, I know some templating libraries(e.g. Eigen) might be able to do it well, but I am more interested in a Lapack/MKL solution, that said, correct solution using Eigen is also welcome.

Edit, as Christian Clason pointed out: C = A* inv(B) can be computed via gesv(trans(B),trans(A)), However, in a real implementation with lapack, gesv will overwrite its input, so does it mean we have to introduce temp variable in order to have C = A * inv(B) correctly computed in semantic sense.

• Note that the columns of $A^{-1}B$ are given by $A^{-1}$ applied to the columns of $B$, so to obtain $X:=A^{-1}B$ you can apply any of the standard methods solve the system $AX=B$ (the LAPACK routines take multiple right hand sides as a matrix, so should work out of the box). For $AB^{-1}$, note that $(AB^{-1})^T = (B^T)^{-1}A^T$. – Christian Clason Sep 17 '16 at 10:49
• @ChristianClason I think you made a good point, so we can still use the same lapack routine with complex conjugate in case of complex value? – lorniper Sep 17 '16 at 13:07
• Yes, real or complex shouldn't make a difference here. – Christian Clason Sep 17 '16 at 13:38
• @ChristianClason maybe you can leave your answer below, on one hand, make sure "real or complex shouldn't make a difference" means trans() for real and conj() for complex, on the other hand, I am still struggling for a semantically reasonable solution because gesv routines overwrite the input A and B. – lorniper Sep 17 '16 at 13:45
• Maybe you could first edit your question to make it a bit more explicit what you actually want to compute; I was guessing that you really needed the matrix $C=AB^{-1}$ (where $B$ is arbitrary and not, as in your first block, the inverse of $A$) for some reason. If all you need it for is to apply it to some vector(s) $x$ later, it is usually far cheaper to first solve $By=x$ and then compute $Ay$. (Similarly, if you want to solve $Cx=b$, you can do this in two stages, where the first is really cheap since you already have the inverse of $B^{-1}$.) – Christian Clason Sep 17 '16 at 14:06

Here is the Eigen example you requested. It follows the approach in Christian Clason's comment and also does the same computation using an explicit inverse for comparison.

#include <iostream>
using std::cout;
using std::endl;
#include <ctime>

#include <Eigen/Core>
#include <Eigen/LU>

void compareInvAndFactor()
{

typedef Eigen::MatrixXd Matrix;
const int n = 1000, m = 600;
Matrix B = Matrix::Random(n, n);
Matrix A = Matrix::Random(m, n);

std::clock_t start = std::clock();
auto Clu = B.transpose().lu().solve(A.transpose()).transpose();
cout << "Elapsed time using LU factorization = " << (std::clock() - start) /
(double)CLOCKS_PER_SEC << " seconds." << endl;

start = std::clock();
auto Cinv = A*B.inverse();
cout << "Elapsed time using inverse = " << (std::clock() - start) /
(double)CLOCKS_PER_SEC << " seconds." << endl;

double diff = (Cinv - Clu).cwiseAbs().maxCoeff();
cout << "max difference in C matrices =" << diff << endl;
}

• The internal implementation for A*B.inverse() is really inverse B first and then do the gemm? – lorniper Sep 17 '16 at 15:14
• Yes, in C++ the precedence of the dot(.) operator is much higher than *. – Bill Greene Sep 17 '16 at 15:31