CODEWITHSUNDEEP
c++matrixlapackinversion
Kevin Synnott
Kevin Synnott

Reputation: 69

Inverting Matrices using LAPACK and timing

I am using LAPACK on C++ to invert a complex matrix. Specifically, the two functions I am using are:

zgetrf for the LU decomposition.

zgetri for the inversion.

Now as my goal to optimize my code, I have a question about the timing. Using a general matrix inversion method with LAPACK (if you have better/quicker functions to use please let me know), is the timing of the function(s) independent of the values in the matrix?

For example, would it be quicker to invert the identity matrix than to invert a densely populated matrix?

Again, I would like to stress that I am asking this question with regards to general LAPACK inversion of complex matrices. I know about the various tri-diagional and banded function one can use.

I am assuming that all the elements in the matrix are complex doubles.

Thanks, Kevin

Upvotes: 3

Views: 972

Answers (1)

Avi Ginsburg
Avi Ginsburg

Reputation: 10596

As Kieran Cooney assumed, LAPACK inverts the identity matrix orders 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;
}

Upvotes: 2

Related Questions