3

我在 C++ 上使用 LAPACK 来反转一个复杂的矩阵。具体来说，我使用的两个功能是：

zgetrf用于 LU 分解。

zgetri为反转。

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

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

再次，我想强调一下，我问的是关于复矩阵的一般 LAPACK 求逆的问题。我知道可以使用的各种三对角和带状功能。

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

谢谢，凯文

4

1 回答 1

2

正如 Kieran Cooney 所假设的，LAPACK 反转单位矩阵的速度比随机密集矩阵快几个数量级。使用下面的测试给了我以下结果（样本大小 = 1，但证明了这一点）：

调整大小的
信息：0
总时间（随机）= 2389 毫秒。
信息：0
总时间（身份）= 14 毫秒。 

#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;
}
于 2015-08-10T13:52:48.887 回答