lapack_calls.cpp

#include <complex>
#include <complex.h>

#ifndef INCLUDE_LAPACK_CALLS_H_
#include "lapacke.h"
#include "../include/lapack_calls.h"
#endif

void zinvert(std::complex<double> **mat, lapack_int n, int &jer) {
  jer = 0;
  __complex__ double *arr = new __complex__ double[n * n];
  const __complex__ double uim = 1.0*I;
  for (lapack_int i = 0; i < n; i++) {
    for (lapack_int j = 0; j < n; j++) {
      lapack_int idx = i + n * j;
      arr[idx] = mat[j][i].real() + uim * mat[j][i].imag();
    }
  }
  
  lapack_int* IPIV = new lapack_int[n]();
  
  LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, arr, n, IPIV);
  LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, arr, n, IPIV);
  for (lapack_int i = 0; i < n; i++) {
    for (lapack_int j = 0; j < n; j++) {
      lapack_int idx = i + n * j;
      mat[j][i] = std::complex<double>(creal(arr[idx]), cimag(arr[idx]));
    }
  }
  delete [] IPIV;
  delete [] arr;
}