Inverting matrices using LAPACK and time
I am using LAPACK in C ++ to invert a complex matrix. Specifically, two functions that I use:
zgetrf
to decompose LU.
zgetri
for inversion.
Now when I want to optimize my code, I have a timing question. Using the general matrix inversion method with LAPACK (if you have better / faster functions to use, please let me know), is the time of the function (s) independent of the values in the matrix?
For example, would it be faster to invert the unit of a matrix than to invert a densely filled matrix?
Again, I would like to emphasize that I am asking this question regarding the general LAPACK inverse of complex matrices. I am aware of the various tri-band and striped functions that you can use.
I am assuming that all elements of the matrix are complex twins.
Thanks, Kevin
source to share
As Kiran Cooney suggested, LAPACK inverts the identity matrix an order of magnitude faster than a random dense matrix. Using the test below gives me the following results (sample size = 1, but proves the point):
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;
}
source to share