From 1487499dde2734d61ec64e2738e3f45c9abf4599 Mon Sep 17 00:00:00 2001
From: "Mulas, Giacomo" <gmulas@oa-cagliari.inaf.it>
Date: Wed, 22 Jan 2025 16:36:13 +0100
Subject: [PATCH] clean up DEBUG_REFINE stuff
---
src/cluster/cluster.cpp | 2 +-
src/include/algebraic.h | 4 ++-
src/include/magma_calls.h | 6 +++--
src/inclusion/inclusion.cpp | 2 +-
src/libnptm/algebraic.cpp | 8 ++++--
src/libnptm/magma_calls.cpp | 49 ++++++++++++++++++++++++++++---------
6 files changed, 53 insertions(+), 18 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 8950265a..3a5dfd19 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -841,7 +841,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
#endif
// we put the accuracygoal in, get the actual accuracy back out
double actualaccuracy = cid->accuracygoal;
- invert_matrix(cid->am, ndit, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, mxndm, cid->proc_device);
+ invert_matrix(cid->am, ndit, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device);
// in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not
#ifdef USE_REFINEMENT
if (cid->refinemode==2) {
diff --git a/src/include/algebraic.h b/src/include/algebraic.h
index 344e78c6..508cf0a3 100644
--- a/src/include/algebraic.h
+++ b/src/include/algebraic.h
@@ -29,6 +29,8 @@
#ifndef INCLUDE_ALGEBRAIC_H_
#define INCLUDE_ALGEBRAIC_H_
+using namespace std;
+
/*! \brief Perform in-place matrix inversion.
*
* \param mat: `complex double **` The matrix to be inverted (must be a square matrix).
@@ -38,6 +40,6 @@
* optional, defaults to 0).
* \param target_device: `int` ID of target GPU, if available (defaults to 0).
*/
-void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, np_int max_size=0, int target_device=0);
+void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, const string& output_path, int jxi488, np_int max_size=0, int target_device=0);
#endif
diff --git a/src/include/magma_calls.h b/src/include/magma_calls.h
index a38c6148..1e1ccff6 100644
--- a/src/include/magma_calls.h
+++ b/src/include/magma_calls.h
@@ -20,6 +20,8 @@
*
*/
+#include <string>
+
#ifndef INCLUDE_MAGMA_CALLS_H_
#define INCLUDE_MAGMA_CALLS_H_
@@ -57,7 +59,7 @@ void magma_zinvert1(dcomplex * &inva, np_int n, int &jer, int device_id);
* \param accuracygoal: `double` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy
* \param device_id: `int` ID of the device for matrix inversion offloading.
*/
-void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id);
+void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488);
/*! \brief apply iterative refinement of the solution of a matrix inversion
*
@@ -71,6 +73,6 @@ void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxrefite
* \param accuracygoal: `double` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy
* \param device_id: `int` ID of the device for matrix inversion offloading.
*/
-void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id);
+void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488);
#endif
diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 8b01327c..fb9fab49 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -1480,7 +1480,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
#endif
// we the accuracygoal in, get the actual accuracy back out
double actualaccuracy = cid->accuracygoal;
- invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, mxndm, cid->proc_device);
+ invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device);
// in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not
#ifdef USE_REFINEMENT
if (cid->refinemode==2) {
diff --git a/src/libnptm/algebraic.cpp b/src/libnptm/algebraic.cpp
index 1df84050..e0363e21 100644
--- a/src/libnptm/algebraic.cpp
+++ b/src/libnptm/algebraic.cpp
@@ -18,6 +18,10 @@
*
* \brief Implementation of algebraic functions with different call-backs.
*/
+
+#include <string>
+using namespace std;
+
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif
@@ -58,14 +62,14 @@ extern void lucin(dcomplex **mat, np_int max_size, np_int size, int &ier);
using namespace std;
-void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, np_int max_size, int target_device) {
+void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, const string& output_path, int jxi488, np_int max_size, int target_device) {
ier = 0;
#ifdef USE_MAGMA
#ifdef USE_REFINEMENT
// try using the iterative refinement to obtain a more accurate solution
// we pass to magma_zinvert_and_refine() the accuracygoal in, get the actual
// accuracy back out
- magma_zinvert_and_refine(mat, size, ier, maxrefiters, accuracygoal, refinemode, target_device);
+ magma_zinvert_and_refine(mat, size, ier, maxrefiters, accuracygoal, refinemode, target_device, output_path, jxi488);
#else
magma_zinvert(mat, size, ier, target_device);
#endif
diff --git a/src/libnptm/magma_calls.cpp b/src/libnptm/magma_calls.cpp
index e1f916e4..22e32c89 100644
--- a/src/libnptm/magma_calls.cpp
+++ b/src/libnptm/magma_calls.cpp
@@ -18,6 +18,9 @@
*
* \brief Implementation of the interface with MAGMA libraries.
*/
+
+using namespace std;
+
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif
@@ -30,8 +33,8 @@
// hand define some preprocessor flags
//#define USE_ZGESV_GPU 1
-#define USE_ZGESV_RBT 1
-#define DEBUG_REFINE 1
+//#define USE_ZGESV_RBT 1
+//#define DEBUG_REFINE 1
#ifdef DEBUG_REFINE
#include <hdf5.h>
@@ -141,7 +144,7 @@ void magma_zinvert1(dcomplex * &inva, np_int n, int &jer, int device_id) {
jer = (int)err;
}
-void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxiters, double &accuracygoal, int refinemode, int device_id) {
+void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488) {
magma_int_t err = MAGMA_SUCCESS;
magma_queue_t queue = NULL;
magma_device_t dev = (magma_device_t)device_id;
@@ -174,10 +177,14 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
magmaDoubleComplex *abuf = (magmaDoubleComplex *) ambuf[0];
memcpy(ambuf[0], a_orig, n*n*sizeof(dcomplex));
VirtualAsciiFile *outam = new VirtualAsciiFile();
- string outam_name = "c_AM0.txt";
+ string outam_name = output_path + "/c_AM4_JXI" + to_string(jxi488) + ".txt";
sprintf(virtual_line, " AM matrix before inversion\n");
outam->append_line(virtual_line);
+#ifdef USE_ILP64
+ sprintf(virtual_line, " %ld\n", n);
+#else
sprintf(virtual_line, " %d\n", n);
+#endif
outam->append_line(virtual_line);
sprintf(virtual_line, " I1+1 I2+1 Real Imag\n");
outam->append_line(virtual_line);
@@ -186,10 +193,14 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
delete outam;
memcpy(ambuf[0], a, n*n*sizeof(dcomplex));
outam = new VirtualAsciiFile();
- outam_name = "c_AM1_0.txt";
+ outam_name = output_path + "/c_AM5_JXI" + to_string(jxi488) + "_iter0.txt";
sprintf(virtual_line, " AM inverse before refinement\n");
outam->append_line(virtual_line);
+#ifdef USE_ILP64
+ sprintf(virtual_line, " %ld\n", n);
+#else
sprintf(virtual_line, " %d\n", n);
+#endif
outam->append_line(virtual_line);
sprintf(virtual_line, " I1+1 I2+1 Real Imag\n");
outam->append_line(virtual_line);
@@ -235,10 +246,14 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
#ifdef DEBUG_REFINE
magma_zgetmatrix(m, m, d_a_residual, m, abuf, m, queue); // copy residual to buffer
outam = new VirtualAsciiFile();
- outam_name = "c_AM2_0.txt";
+ outam_name = output_path + "/c_AM6_JXI" + to_string(jxi488) + "_iter0.txt";
sprintf(virtual_line, " AM residual before refinement\n");
outam->append_line(virtual_line);
+#ifdef USE_ILP64
+ sprintf(virtual_line, " %ld\n", n);
+#else
sprintf(virtual_line, " %d\n", n);
+#endif
outam->append_line(virtual_line);
sprintf(virtual_line, " I1+1 I2+1 Real Imag\n");
outam->append_line(virtual_line);
@@ -266,10 +281,14 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
#ifdef DEBUG_REFINE
magma_zgetmatrix(m, m, d_a_refine, m, abuf, m, queue); // copy refinement to buffer
outam = new VirtualAsciiFile();
- outam_name = "c_AM3_" + to_string(iter+1) + ".txt";
+ outam_name = output_path + "/c_AM7_JXI" + to_string(jxi488) + "_iter" + to_string(iter+1) + ".txt";
sprintf(virtual_line, " AM correction refinement %d\n", iter+1);
outam->append_line(virtual_line);
+#ifdef USE_ILP64
+ sprintf(virtual_line, " %ld\n", n);
+#else
sprintf(virtual_line, " %d\n", n);
+#endif
outam->append_line(virtual_line);
sprintf(virtual_line, " I1+1 I2+1 Real Imag\n");
outam->append_line(virtual_line);
@@ -282,10 +301,14 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
#ifdef DEBUG_REFINE
magma_zgetmatrix(m, m, d_a, m, abuf, m, queue); // copy new inverse to buffer
outam = new VirtualAsciiFile();
- outam_name = "c_AM1_" + to_string(iter+1) + ".txt";
+ outam_name = output_path + "/c_AM5_JXI" + to_string(jxi488) + "_iter" + to_string(iter+1) + ".txt";
sprintf(virtual_line, " AM inverse after refinement %d\n", iter+1);
outam->append_line(virtual_line);
+#ifdef USE_ILP64
+ sprintf(virtual_line, " %ld\n", n);
+#else
sprintf(virtual_line, " %d\n", n);
+#endif
outam->append_line(virtual_line);
sprintf(virtual_line, " I1+1 I2+1 Real Imag\n");
outam->append_line(virtual_line);
@@ -300,10 +323,14 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
#ifdef DEBUG_REFINE
magma_zgetmatrix(m, m, d_a_residual, m, abuf, m, queue); // copy residual to buffer
outam = new VirtualAsciiFile();
- outam_name = "c_AM2_" + to_string(iter+1) + ".txt";
+ outam_name = output_path + "/c_AM6_JXI" + to_string(jxi488) + "_iter" + to_string(iter+1) + ".txt";
sprintf(virtual_line, " AM residual after refinement %d\n", iter+1);
outam->append_line(virtual_line);
+#ifdef USE_ILP64
+ sprintf(virtual_line, " %ld\n", n);
+#else
sprintf(virtual_line, " %d\n", n);
+#endif
outam->append_line(virtual_line);
sprintf(virtual_line, " I1+1 I2+1 Real Imag\n");
outam->append_line(virtual_line);
@@ -345,13 +372,13 @@ void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxi
jer = (int)err;
}
-void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxiters, double &accuracygoal, int refinemode, int device_id) {
+void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488) {
if (maxiters>0) {
dcomplex *inva = new dcomplex[n*n]();
dcomplex *aorig = mat[0];
memcpy(inva, aorig, n*n*sizeof(dcomplex));
magma_zinvert1(inva, n, jer, device_id);
- magma_refine(aorig, inva, n, jer, maxiters, accuracygoal, refinemode, device_id);
+ magma_refine(aorig, inva, n, jer, maxiters, accuracygoal, refinemode, device_id, output_path, jxi488);
memcpy(aorig, inva, n*n*sizeof(dcomplex));
delete[] inva;
}
--
GitLab