适当的 Lapack/MKL 例程来有效地计算 C = A* inv(B)

计算科学 线性代数 拉帕克 英特尔-mkl
2021-12-12 10:17:23

我知道,从数字的角度来看,计算

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

是完全不同的东西,我们应该使用 TRF 例程分解矩阵,然后使用 TRS/TRI 例程或组合的 gesv 调用求解。

无论如何,是否有一个有效的例程来专门处理

C = A* inv(B)

顺便说一句,我知道一些模板库(例如 Eigen)可能会做得很好,但我对 Lapack/MKL 解决方案更感兴趣,也就是说,也欢迎使用 Eigen 的正确解决方案。

编辑,正如 Christian Clason 指出的那样:C = A* inv(B) 可以通过 gesv(trans(B),trans(A)) 计算,但是,在使用 lapack 的实际实现中,gesv 将覆盖其输入,因此也是如此这意味着我们必须引入临时变量才能在语义上正确计算 C = A * inv(B)。

1个回答

这是您要求的本征示例。它遵循 Christian Clason 评论中的方法,并且使用显式逆进行比较进行相同的计算。

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