使用 LAPACK 和时序求逆矩阵
Inverting Matrices using LAPACK and timing
我在 C++ 上使用 LAPACK 来求逆复杂矩阵。具体来说,我使用的两个函数是:
zgetrf
用于 LU 分解。
zgetri
为反转。
现在作为我优化代码的目标,我有一个关于时间的问题。使用 LAPACK 的通用矩阵求逆方法(如果您有 better/quicker 函数要使用,请告诉我),函数的时间是否独立于矩阵中的值?
例如,反转单位矩阵会比反转密集矩阵更快吗?
再次强调,我问的这个问题是关于复数矩阵的一般 LAPACK 求逆。我知道可以使用的各种三角函数和带状函数。
我假设矩阵中的所有元素都是复数双精度数。
谢谢,
凯文
正如 Kieran Cooney 假设的那样,LAPACK 反转单位矩阵的速度比随机密集矩阵快几个数量级。使用下面的测试给出了以下结果(样本量 = 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;
}
我在 C++ 上使用 LAPACK 来求逆复杂矩阵。具体来说,我使用的两个函数是:
zgetrf
用于 LU 分解。
zgetri
为反转。
现在作为我优化代码的目标,我有一个关于时间的问题。使用 LAPACK 的通用矩阵求逆方法(如果您有 better/quicker 函数要使用,请告诉我),函数的时间是否独立于矩阵中的值?
例如,反转单位矩阵会比反转密集矩阵更快吗?
再次强调,我问的这个问题是关于复数矩阵的一般 LAPACK 求逆。我知道可以使用的各种三角函数和带状函数。
我假设矩阵中的所有元素都是复数双精度数。
谢谢, 凯文
正如 Kieran Cooney 假设的那样,LAPACK 反转单位矩阵的速度比随机密集矩阵快几个数量级。使用下面的测试给出了以下结果(样本量 = 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;
}