利用LAPACK和时序反演矩阵

Inverting Matrices using LAPACK and timing

本文关键字：LAPACK 利用      更新时间：2023-10-16

我在c++上使用LAPACK来反转一个复矩阵。具体来说，我使用的两个函数是:

zgetrf为LU分解。

zgetri为反演。

现在我的目标是优化我的代码，我有一个关于时间的问题。使用LAPACK的一般矩阵反演方法(如果您有更好/更快的函数可以使用，请告诉我)，函数的时间是否与矩阵中的值无关?

例如，对单位矩阵求逆是否比对密集矩阵求逆更快?

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

我假设矩阵中的所有元素都是复双精度。

谢谢,凯文

正如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 << "Resizedn";
    // 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;
}
相关文章：