From 69932a5815e2e14c279f1a9ea0dce640a8ea6af1 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Mon, 4 Mar 2024 15:49:31 +0100
Subject: [PATCH] Create a LAPACK interface library
---
src/include/lapack_calls.h | 20 ++++++++++++++++++++
src/libnptm/lapack_calls.cpp | 34 ++++++++++++++++++++++++++++++++++
2 files changed, 54 insertions(+)
create mode 100644 src/include/lapack_calls.h
create mode 100644 src/libnptm/lapack_calls.cpp
diff --git a/src/include/lapack_calls.h b/src/include/lapack_calls.h
new file mode 100644
index 00000000..0270a05c
--- /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 00000000..ded76f3c
--- /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;
+}
--
GitLab