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

Question

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.

Answer 1

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;
}