c++ - 使用 LAPACK 和时序求逆矩阵

Question

我在 C++ 上使用 LAPACK 来求逆复杂矩阵。具体来说，我使用的两个函数是:

zgetrf 用于 LU 分解。

zgetri 用于反转。

现在作为我优化代码的目标，我有一个关于时间的问题。使用带有 LAPACK 的通用矩阵求逆方法(如果您有更好/更快的函数，请告诉我)，函数的时间是否独立于矩阵中的值？

例如，反转单位矩阵会比反转密集矩阵更快吗？

再次强调，我问的这个问题是关于复数矩阵的一般 LAPACK 求逆。我知道可以使用的各种三角函数和带状函数。

我假设矩阵中的所有元素都是复数 double 值。

谢谢，凯文

Answer 1

正如 Kieran Cooney 所假设的那样，LAPACK 反转单位矩阵的速度比随 secret 集矩阵快几个数量级。使用下面的测试给出了以下结果(样本量 = 1，但证明了这一点):

Resized
Info: 0
Total Time (random) = 2389 milliseconds.
Info: 0
Total Time (identity) = 14 milliseconds.

#include "lapacke.h"

#include <iostream>
#include <vector>
#include <Eigen/Core>
#include <chrono>

lapack_int getSize(lapack_int n, lapack_complex_double* a,
    const lapack_int* ipiv, lapack_complex_double* work)
{
    lapack_complex_double resizetome;
    lapack_int hello = -1;
    lapack_int info = -1;

    LAPACK_zgetri(&n, a, &n, ipiv, &resizetome, &hello, &info);

    return lapack_int(resizetome.real());

}
void invert(lapack_int n, lapack_complex_double* a,
    lapack_int* ipiv, lapack_complex_double* work, lapack_int lwork, lapack_int *info)
{
    // LU factor
    LAPACK_zgetrf(&n, &n, a, &n, ipiv, info);

    // Invert
    LAPACK_zgetri(&n, a, &n, ipiv, work, &lwork, info);
}

int main(int argc, char* argv[]) {

    int sz = 1000;

    int ln = sz;
    int llda = sz;
    int lipiv = 1;
    int llwork = -1;
    int linfo = 0;

    srand(time(NULL));

    typedef Eigen::MatrixXcd lapackMat;
    lapackMat ident = lapackMat::Identity(sz, sz).eval();
    lapackMat randm = lapackMat::Random(sz, sz);
    lapackMat work = lapackMat::Zero(1, 1);
    Eigen::VectorXi ipvt(sz);
    randm;

    work.resize(1,
        getSize(ln, randm.data(), ipvt.data(), work.data())
        );

    std::cout << "Resized\n";

    // Timing for random matrix
    {
        auto startTime = std::chrono::high_resolution_clock::now();

        invert(ln, randm.data(), ipvt.data(), work.data(), llwork, &linfo);

        auto endTime = std::chrono::high_resolution_clock::now();

        std::cout << "Info: " << linfo << "\n";

        std::cout << "Total Time (random) = " <<
            std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count()
            << " milliseconds.\n";
    }

    // Timing for identity matrix
    {
        auto startTime = std::chrono::high_resolution_clock::now();

        invert(ln, ident.data(), ipvt.data(), work.data(), llwork, &linfo);

        auto endTime = std::chrono::high_resolution_clock::now();

        std::cout << "Info: " << linfo << "\n";

        std::cout << "Total Time (identity) = " <<
            std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count()
            << " milliseconds.\n";

    }

    return 0;
}