diff --git a/src/include/lapack_calls.h b/src/include/lapack_calls.h new file mode 100644 index 0000000000000000000000000000000000000000..0270a05c148944960c8091b2c45f278c7981db8e --- /dev/null +++ b/src/include/lapack_calls.h @@ -0,0 +1,20 @@ +/*! \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 diff --git a/src/libnptm/lapack_calls.cpp b/src/libnptm/lapack_calls.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ded76f3c5a35f7b122f9aceb8d154c779927b4a6 --- /dev/null +++ b/src/libnptm/lapack_calls.cpp @@ -0,0 +1,34 @@ +#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; +}