Add auxiliary utilities and write AM matrix to debug output

build/Makefile.am

 LDADD=libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}
 lib_LTLIBRARIES=libnptm/libnptm.la
-libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
+libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
 if BUILDFORTRAN
 PROGS=cluster/edfb_clu cluster/clu cluster/np_cluster sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS
 bin_PROGRAMS=$(PROGS)

build/Makefile.in

@@ -156,8 +156,9 @@ am_libnptm_libnptm_la_OBJECTS = ../src/libnptm/algebraic.lo \
 	../src/libnptm/Configuration.lo ../src/libnptm/file_io.lo \
 	../src/libnptm/lapack_calls.lo ../src/libnptm/logging.lo \
 	../src/libnptm/magma_calls.lo ../src/libnptm/Parsers.lo \
-	../src/libnptm/sph_subs.lo ../src/libnptm/tfrfme.lo \
-	../src/libnptm/TransitionMatrix.lo ../src/libnptm/tra_subs.lo
+	../src/libnptm/sph_subs.lo ../src/libnptm/utils.lo \
+	../src/libnptm/tfrfme.lo ../src/libnptm/TransitionMatrix.lo \
+	../src/libnptm/tra_subs.lo
 libnptm_libnptm_la_OBJECTS = $(am_libnptm_libnptm_la_OBJECTS)
 AM_V_lt = $(am__v_lt_@AM_V@)
 am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@)
@@ -304,6 +305,7 @@ am__depfiles_remade = ../src/cluster/$(DEPDIR)/cluster.Po \
 	../src/libnptm/$(DEPDIR)/sph_subs.Plo \
 	../src/libnptm/$(DEPDIR)/tfrfme.Plo \
 	../src/libnptm/$(DEPDIR)/tra_subs.Plo \
+	../src/libnptm/$(DEPDIR)/utils.Plo \
 	../src/sphere/$(DEPDIR)/np_sphere.Po \
 	../src/sphere/$(DEPDIR)/sphere.Po \
 	../src/testing/$(DEPDIR)/test_TEDF.Po \
@@ -577,7 +579,7 @@ top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
 LDADD = libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}
 lib_LTLIBRARIES = libnptm/libnptm.la
-libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
+libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
 @BUILDFORTRAN_FALSE@PROGS = cluster/np_cluster sphere/np_sphere trapping/np_trapping testing/test_TEDF testing/test_TTMS
 @BUILDFORTRAN_TRUE@PROGS = cluster/edfb_clu cluster/clu cluster/np_cluster sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS
 @BUILDFORTRAN_TRUE@EDFBCLUSOURCES = ../src/cluster/edfb_clu.f
@@ -749,6 +751,8 @@ clean-libLTLIBRARIES:
 	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
 ../src/libnptm/sph_subs.lo: ../src/libnptm/$(am__dirstamp) \
 	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
+../src/libnptm/utils.lo: ../src/libnptm/$(am__dirstamp) \
+	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
 ../src/libnptm/tfrfme.lo: ../src/libnptm/$(am__dirstamp) \
 	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
 ../src/libnptm/TransitionMatrix.lo: ../src/libnptm/$(am__dirstamp) \
@@ -900,6 +904,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/sph_subs.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/tfrfme.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/tra_subs.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/utils.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/sphere/$(DEPDIR)/np_sphere.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/sphere/$(DEPDIR)/sphere.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/testing/$(DEPDIR)/test_TEDF.Po@am__quote@ # am--include-marker
@@ -1268,6 +1273,7 @@ distclean: distclean-am
 	-rm -f ../src/libnptm/$(DEPDIR)/sph_subs.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tfrfme.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tra_subs.Plo
+	-rm -f ../src/libnptm/$(DEPDIR)/utils.Plo
 	-rm -f ../src/sphere/$(DEPDIR)/np_sphere.Po
 	-rm -f ../src/sphere/$(DEPDIR)/sphere.Po
 	-rm -f ../src/testing/$(DEPDIR)/test_TEDF.Po
@@ -1337,6 +1343,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f ../src/libnptm/$(DEPDIR)/sph_subs.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tfrfme.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tra_subs.Plo
+	-rm -f ../src/libnptm/$(DEPDIR)/utils.Plo
 	-rm -f ../src/sphere/$(DEPDIR)/np_sphere.Po
 	-rm -f ../src/sphere/$(DEPDIR)/sphere.Po
 	-rm -f ../src/testing/$(DEPDIR)/test_TEDF.Po

src/cluster/cluster.cpp

@@ -83,6 +83,10 @@
 #include "../include/file_io.h"
 #endif
 
+#ifndef INCLUDE_UTILS_H_
+#include "../include/utils.h"
+#endif
+
 using namespace std;
 
 // I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
@@ -771,6 +775,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   elapsed = interval_end - interval_start;
   message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
   logger->log(message);
+  if (jxi488 == 1) {
+    string am_file_name = "c_am_jxi_1.txt";
+    write_dcomplex_matrix(am_file_name, cid->am, 624, 624);
+  }
   interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
   nvtxRangePush("Invert the matrix");

src/include/utils.h

+/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   A copy of the GNU General Public License is distributed along with
+   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
+ */
+
+/*! \file utils.h
+ *
+ * \brief Definition of auxiliary code utilities.
+ */
+
+#ifndef INCLUDE_UTILS_H_
+#define INCLUDE_UTILS_H_
+
+/*! \brief Write a double complex matrix to a text file.
+ *
+ * \param file_name: `const string&` Name of the file to be written.
+ * \param mat: `dcomplex **` Pointer to the matrix.
+ * \param rows: `np_int` Number of rows in the matrix.
+ * \param columns: `np_int` Number of columns in the matrix.
+ */
+int write_dcomplex_matrix(const std::string& file_name, dcomplex **mat, np_int rows, np_int columns);
+
+#endif

src/libnptm/utils.cpp

+/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   A copy of the GNU General Public License is distributed along with
+   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
+ */
+
+/*! \file utils.cpp
+ *
+ * \brief Implementation of auxiliary code utilities.
+ */
+
+/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   A copy of the GNU General Public License is distributed along with
+   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
+ */
+
+/*! \file utils.h
+ *
+ * \brief Definition of auxiliary code utilities.
+ */
+
+#include <cstdio>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+#ifndef INCLUDE_UTILS_H_
+#include "../include/utils.h"
+#endif
+
+using namespace std;
+
+int write_dcomplex_matrix(const std::string& file_name, dcomplex **mat, np_int rows, np_int columns) {
+  int result = 0;
+  FILE* output = fopen(file_name.c_str(), "w");
+  if (output) {
+    string str_format = "%5d%5d (%13.5lE,%13.5lE)\n";
+    if (sizeof(np_int) > sizeof(int)) {
+      str_format = "%5ld%5ld (%13.5lE,%13.5lE)\n";
+    }
+    for (np_int i = 0; i < rows; i++) {
+      for (np_int j = 0; j < columns; j++) {
+	fprintf(output, str_format.c_str(), i, j, real(mat[i][j]), imag(mat[i][j]));
+      }
+    }
+    fclose(output);
+  } else {
+    // Could not open the output file.
+    result = 1;
+  }
+  return result;
+}