Create a LAPACK interface library

/*! \file lapack_calss.h

 *

 * \brief C++ interface to LAPACK calls.

 *

 */

#ifndef INCLUDE_LAPACK_CALLS_H_

#define INCLUDE_LAPACK_CALLS_H_

/*! \brief Invert a complex matrix with double precision elements.

 *

 * Use LAPACKE64 to perform an in-place matrix inversion for a complex

 * matrix with double precision elements.

 *

 * \param mat: Matrix of complex. The matrix to be inverted.

 * \param n: `int` The number of rows and columns of the [n x n] matrix.

 */

void zinvert(std::complex<double> **mat, int n);

#endif

#include <complex>

#include <cstdlib>

#include "lapacke.h"

#ifndef INCLUDE_LAPACK_CALLS_H_

#define INCLUDE_LAPACK_CALLS_H_

#include "../include/lapac_calls.h"

#endif

using namespace std;

void zinvert(std::complex<double> **mat, int n){

  __complex__ double *arr = new __complex__ double[n * n];

  const __complex__ double uim = 1.0di;

  for (int i = 0; i < n; i++) {

    for (int j = 0; j < n; j++) {

      int idx = i * n + j;

      arr[idx] = mat[i][j].real() + uim * mat[i][j].imag();

    }

  }

  int* IPIV = new int[n];

  LAPACKE_zgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);

  LAPACKE_zgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);

  for (int i = 0; i < n; i++) {

    for (int j = 0; j < n; j++) {

      int idx = i * n + j;

      mat[i][j] = complex<double>(__real__ arr[idx], __imag__ arr[idx]);

    }

  }

  delete [] IPIV;

  delete [] arr;

}