From 89d2e692ddeed8349ac2814c77e7fb7f1aad1e15 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Mon, 11 Mar 2024 20:31:34 +0100
Subject: [PATCH] Map C1 matrices into contiguous stacks

---
 src/include/Commons.h   | 463 ++++++++++++++++++++--------------------
 src/libnptm/Commons.cpp |  21 +-
 2 files changed, 249 insertions(+), 235 deletions(-)

diff --git a/src/include/Commons.h b/src/include/Commons.h
index e9545444..1994a292 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -36,153 +36,162 @@
  */
 class C1 {
 protected:
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief Maximum order of field expansion.
-	int lm;
-	//! \brief NLMMT. QUESTION: definition?
-	int nlmmt;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief Maximum order of field expansion.
+  int lm;
+  //! \brief NLMMT. QUESTION: definition?
+  int nlmmt;
+  //! \brief Contiguous RMI vector.
+  dcomplex *vec_rmi;
+  //! \brief Contiguous REI vector.
+  dcomplex *vec_rei;
+  //! \brief Contiguous W vector.
+  dcomplex *vec_w;
+  //! \brief Contiguous VINTS vector
+  dcomplex *vec_vints;
+  
 public:
-	//! \brief QUESTION: definition?
-	dcomplex **rmi;
-	//! \brief QUESTION: definition?
-	dcomplex **rei;
-	//! \brief QUESTION: definition?
-	dcomplex **w;
-	//! \brief QUESTION: definition?
-	dcomplex *fsas;
-	//! \brief QUESTION: definition?
-	dcomplex **vints;
-	//! \brief QUESTION: definition?
-	double *sscs;
-	//! \brief QUESTION: definition?
-	double *sexs;
-	//! \brief QUESTION: definition?
-	double *sabs;
-	//! \brief QUESTION: definition?
-	double *sqscs;
-	//! \brief QUESTION: definition?
-	double *sqexs;
-	//! \brief QUESTION: definition?
-	double *sqabs;
-	//! \brief QUESTION: definition?
-	double *gcsv;
-	//! \brief Vector of sphere X coordinates.
-	double *rxx;
-	//! \brief Vector of sphere X coordinates.
-	double *ryy;
-	//! \brief Vector of sphere X coordinates.
-	double *rzz;
-	//! \brief Vector of sphere radii.
-	double *ros;
-	//! \brief Matrix of spherical layer break radii. QUESTION: correct?
-	double **rc;
-	//! \brief Vector of spherical component identifiers.
-	int *iog;
-	//! \brief Vector of numbers of spherical layers.
-	int *nshl;
-	//! \brief QUESTION: definition?
-	dcomplex ***sas;
+  //! \brief QUESTION: definition?
+  dcomplex **rmi;
+  //! \brief QUESTION: definition?
+  dcomplex **rei;
+  //! \brief QUESTION: definition?
+  dcomplex **w;
+  //! \brief QUESTION: definition?
+  dcomplex *fsas;
+  //! \brief QUESTION: definition?
+  dcomplex **vints;
+  //! \brief QUESTION: definition?
+  double *sscs;
+  //! \brief QUESTION: definition?
+  double *sexs;
+  //! \brief QUESTION: definition?
+  double *sabs;
+  //! \brief QUESTION: definition?
+  double *sqscs;
+  //! \brief QUESTION: definition?
+  double *sqexs;
+  //! \brief QUESTION: definition?
+  double *sqabs;
+  //! \brief QUESTION: definition?
+  double *gcsv;
+  //! \brief Vector of sphere X coordinates.
+  double *rxx;
+  //! \brief Vector of sphere X coordinates.
+  double *ryy;
+  //! \brief Vector of sphere X coordinates.
+  double *rzz;
+  //! \brief Vector of sphere radii.
+  double *ros;
+  //! \brief Matrix of spherical layer break radii. QUESTION: correct?
+  double **rc;
+  //! \brief Vector of spherical component identifiers.
+  int *iog;
+  //! \brief Vector of numbers of spherical layers.
+  int *nshl;
+  //! \brief QUESTION: definition?
+  dcomplex ***sas;
 
-	/*! \brief C1 instance constructor.
-	 *
-	 * \param ns: `int` Number of spheres.
-	 * \param l_max: `int` Maximum order of field expansion.
-	 * \param nshl: `int *` Array of number of layers in spheres.
-	 * \param iog: `int *` Vector of spherical units ID numbers.
-	 */
-	C1(int ns, int l_max, int *nshl, int *iog);
+  /*! \brief C1 instance constructor.
+   *
+   * \param ns: `int` Number of spheres.
+   * \param l_max: `int` Maximum order of field expansion.
+   * \param nshl: `int *` Array of number of layers in spheres.
+   * \param iog: `int *` Vector of spherical units ID numbers.
+   */
+  C1(int ns, int l_max, int *nshl, int *iog);
 
-	//! \brief C1 instance destroyer.
-	~C1();
+  //! \brief C1 instance destroyer.
+  ~C1();
 };
 
 /*! \brief Representation of the FORTRAN C2 blocks.
  *
  */
 class C2 {
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief Number of required orders.
-	int nhspo;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief Number of required orders.
+  int nhspo;
 public:
-	//! \brief QUESTION: definition?
-	dcomplex *ris;
-	//! \brief QUESTION: definition?
-	dcomplex *dlri;
-	//! \brief QUESTION: definition?
-	dcomplex *dc0;
-	//! \brief QUESTION: definition?
-	dcomplex *vkt;
-	//! Vector of scaled sizes. QUESTION: correct?
-	double *vsz;
+  //! \brief QUESTION: definition?
+  dcomplex *ris;
+  //! \brief QUESTION: definition?
+  dcomplex *dlri;
+  //! \brief QUESTION: definition?
+  dcomplex *dc0;
+  //! \brief QUESTION: definition?
+  dcomplex *vkt;
+  //! Vector of scaled sizes. QUESTION: correct?
+  double *vsz;
 
-	/*! \brief C2 instance constructor.
-	 *
-	 * \param ns: `int` Number of spheres.
-	 * \param nl: `int`
-	 * \param npnt: `int`
-	 * \param npntts: `int`
-	 */
-	C2(int ns, int nl, int npnt, int npntts);
+  /*! \brief C2 instance constructor.
+   *
+   * \param ns: `int` Number of spheres.
+   * \param nl: `int`
+   * \param npnt: `int`
+   * \param npntts: `int`
+   */
+  C2(int ns, int nl, int npnt, int npntts);
 
-	//! \brief C2 instance destroyer.
-	~C2();
+  //! \brief C2 instance destroyer.
+  ~C2();
 };
 
 /*! \brief Representation of the FORTRAN C3 blocks.
  */
 class C3 {
 public:
-	//! \brief QUESTION: definition?
-	dcomplex tfsas;
-	//! \brief QUESTION: definition?
-	dcomplex **tsas;
-	//! \brief QUESTION: definition?
-	double gcs;
-	//! \brief QUESTION: definition?
-	double scs;
-	//! \brief QUESTION: definition?
-	double ecs;
-	//! \brief QUESTION: definition?
-	double acs;
+  //! \brief QUESTION: definition?
+  dcomplex tfsas;
+  //! \brief QUESTION: definition?
+  dcomplex **tsas;
+  //! \brief QUESTION: definition?
+  double gcs;
+  //! \brief QUESTION: definition?
+  double scs;
+  //! \brief QUESTION: definition?
+  double ecs;
+  //! \brief QUESTION: definition?
+  double acs;
 
-	/*! \brief C3 instance constructor.
-	 */
-	C3();
+  /*! \brief C3 instance constructor.
+   */
+  C3();
 
-	/*! \brief C3 instance destroyer.
-	 */
-	~C3();
+  /*! \brief C3 instance destroyer.
+   */
+  ~C3();
 };
 
 /*! \brief Representation of the FORTRAN C4 blocks.
  */
 struct C4 {
-	//! \brief QUESTION: definition?
-	int litpo;
-	//! \brief QUESTION: definition?
-	int litpos;
-	//! \brief Maximum field expansion order plus one. QUESTION: correct?
-	int lmpo;
-	//! \brief Twice maximum field expansion order plus one. QUESTION: correct?
-	int lmtpo;
-	//! \brief Square of `lmtpo`.
-	int lmtpos;
-	//! \brief QUESTION: definition?
-	int li;
-	//! \brief QUESTION: definition?
-	int nlim;
-	//! \brief QUESTION: definition?
-	int le;
-	//! \brief QUESTION: definition?
-	int nlem;
-	//! \brief Maximum field expansion order. QUESTION: correct?
-	int lm;
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief QUESTION: definition?
-	int nv3j;
+  //! \brief QUESTION: definition?
+  int litpo;
+  //! \brief QUESTION: definition?
+  int litpos;
+  //! \brief Maximum field expansion order plus one. QUESTION: correct?
+  int lmpo;
+  //! \brief Twice maximum field expansion order plus one. QUESTION: correct?
+  int lmtpo;
+  //! \brief Square of `lmtpo`.
+  int lmtpos;
+  //! \brief QUESTION: definition?
+  int li;
+  //! \brief QUESTION: definition?
+  int nlim;
+  //! \brief QUESTION: definition?
+  int le;
+  //! \brief QUESTION: definition?
+  int nlem;
+  //! \brief Maximum field expansion order. QUESTION: correct?
+  int lm;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief QUESTION: definition?
+  int nv3j;
 };
 
 /*! \brief Vectors and matrices that are specific to cluster C1 blocks.
@@ -190,131 +199,131 @@ struct C4 {
  */
 class C1_AddOns {
 protected:
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief QUESTION: definition?
-	int nlemt;
-	//! \brief Maximum expansion order plus one. QUESTION: correct?
-	int lmpo;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief QUESTION: definition?
+  int nlemt;
+  //! \brief Maximum expansion order plus one. QUESTION: correct?
+  int lmpo;
 
-	/*! \brief Allocate the necessary common vectors depending on configuration.
-	 *
-	 * The size of the vectors and matrices defined in various common
-	 * blocks, and particularly in C1, depends on many settings of the
-	 * problem configuration, such as the number of spheres, the number
-	 * of layers the spheres are made of, the field expansion order and
-	 * others. This function collects the calculations needed to infer
-	 * the necessary amount of memory for these configurable elements,
-	 * thus making the class constructor more compact and easier to handle.
-	 *
-	 * \param c4: `C4 *` Pointer to a C4 structure.
-	 */
-	void allocate_vectors(C4 *c4);
+  /*! \brief Allocate the necessary common vectors depending on configuration.
+   *
+   * The size of the vectors and matrices defined in various common
+   * blocks, and particularly in C1, depends on many settings of the
+   * problem configuration, such as the number of spheres, the number
+   * of layers the spheres are made of, the field expansion order and
+   * others. This function collects the calculations needed to infer
+   * the necessary amount of memory for these configurable elements,
+   * thus making the class constructor more compact and easier to handle.
+   *
+   * \param c4: `C4 *` Pointer to a C4 structure.
+   */
+  void allocate_vectors(C4 *c4);
 public:
-	//! \brief QUESTION: definition?
-	dcomplex *vh;
-	//! \brief QUESTION: definition?
-	dcomplex *vj0;
-	//! \brief QUESTION: definition?
-	dcomplex *vj;
-	//! \brief QUESTION: definition?
-	dcomplex *vyhj;
-	//! \brief QUESTION: definition?
-	dcomplex *vyj0;
-	//! \brief QUESTION: definition?
-	dcomplex **am0m;
-	//! \brief QUESTION: definition?
-	dcomplex *vint;
-	//! \brief QUESTION: definition?
-	dcomplex *vintm;
-	//! \brief QUESTION: definition?
-	dcomplex **vints;
-	//! \brief QUESTION: definition?
-	dcomplex *vintt;
-	//! \brief QUESTION: definition?
-	dcomplex **fsac;
-	//! \brief QUESTION: definition?
-	dcomplex **sac;
-	//! \brief QUESTION: definition?
-	dcomplex **fsacm;
-	//! \brief QUESTION: definition?
-	double *scsc;
-	//! \brief QUESTION: definition?
-	dcomplex *scscp;
-	//! \brief QUESTION: definition?
-	double *ecsc;
-	//! \brief QUESTION: definition?
-	double *ecscm;
-	//! \brief QUESTION: definition?
-	double *scscm;
-	//! \brief QUESTION: definition?
-	dcomplex *ecscp;
-	//! \brief QUESTION: definition?
-	dcomplex *scscpm;
-	//! \brief QUESTION: definition?
-	dcomplex *ecscpm;
-	//! \brief QUESTION: definition?
-	double *v3j0;
-	//! \brief QUESTION: definition?
-	double *sscs;
-	//! \brief QUESTION: definition?
-	int **ind3j;
+  //! \brief QUESTION: definition?
+  dcomplex *vh;
+  //! \brief QUESTION: definition?
+  dcomplex *vj0;
+  //! \brief QUESTION: definition?
+  dcomplex *vj;
+  //! \brief QUESTION: definition?
+  dcomplex *vyhj;
+  //! \brief QUESTION: definition?
+  dcomplex *vyj0;
+  //! \brief QUESTION: definition?
+  dcomplex **am0m;
+  //! \brief QUESTION: definition?
+  dcomplex *vint;
+  //! \brief QUESTION: definition?
+  dcomplex *vintm;
+  //! \brief QUESTION: definition?
+  dcomplex **vints;
+  //! \brief QUESTION: definition?
+  dcomplex *vintt;
+  //! \brief QUESTION: definition?
+  dcomplex **fsac;
+  //! \brief QUESTION: definition?
+  dcomplex **sac;
+  //! \brief QUESTION: definition?
+  dcomplex **fsacm;
+  //! \brief QUESTION: definition?
+  double *scsc;
+  //! \brief QUESTION: definition?
+  dcomplex *scscp;
+  //! \brief QUESTION: definition?
+  double *ecsc;
+  //! \brief QUESTION: definition?
+  double *ecscm;
+  //! \brief QUESTION: definition?
+  double *scscm;
+  //! \brief QUESTION: definition?
+  dcomplex *ecscp;
+  //! \brief QUESTION: definition?
+  dcomplex *scscpm;
+  //! \brief QUESTION: definition?
+  dcomplex *ecscpm;
+  //! \brief QUESTION: definition?
+  double *v3j0;
+  //! \brief QUESTION: definition?
+  double *sscs;
+  //! \brief QUESTION: definition?
+  int **ind3j;
 
-	/*! \brief C1_AddOns instance constructor.
-	 *
-	 * \param c4: `C4 *` Pointer to a C4 structure.
-	 */
-	C1_AddOns(C4 *c4);
+  /*! \brief C1_AddOns instance constructor.
+   *
+   * \param c4: `C4 *` Pointer to a C4 structure.
+   */
+  C1_AddOns(C4 *c4);
 
-	//! \brief C1_AddOns instance destroyer.
-	~C1_AddOns();
+  //! \brief C1_AddOns instance destroyer.
+  ~C1_AddOns();
 };
 
 /*! \brief Representation of the FORTRAN C6 blocks.
  */
 class C6 {
 public:
-	//! \brief QUESTION: definition?
-	double *rac3j;
+  //! \brief QUESTION: definition?
+  double *rac3j;
 
-	/*! \brief C6 instance constructor.
-	 *
-	 * \param lmtpo: `int` QUESTION: definition?
-	 */
-	C6(int lmtpo);
+  /*! \brief C6 instance constructor.
+   *
+   * \param lmtpo: `int` QUESTION: definition?
+   */
+  C6(int lmtpo);
 
-	/*! \brief C6 instance destroyer.
-	 */
-	~C6();
+  /*! \brief C6 instance destroyer.
+   */
+  ~C6();
 };
 
 /*! \brief Representation of the FORTRAN C9 blocks.
  */
 class C9 {
 protected:
-	//! \brief Number of rows in the GIS and GLS matrices
-	int gis_size_0;
-	//! \brief Number of rows in the SAM matrix
-	int sam_size_0;
+  //! \brief Number of rows in the GIS and GLS matrices
+  int gis_size_0;
+  //! \brief Number of rows in the SAM matrix
+  int sam_size_0;
 public:
-	//! \brief QUESTION: definition?
-	dcomplex **gis;
-	//! \brief QUESTION: definition?
-	dcomplex **gls;
-	//! \brief QUESTION: definition?
-	dcomplex **sam;
+  //! \brief QUESTION: definition?
+  dcomplex **gis;
+  //! \brief QUESTION: definition?
+  dcomplex **gls;
+  //! \brief QUESTION: definition?
+  dcomplex **sam;
 
-	/*! \brief C9 instance constructor.
-	 *
-	 * \param ndi: `int` QUESTION: definition?
-	 * \param nlem: `int` QUESTION: definition?
-	 * \param ndit: `int` QUESTION: definition?
-	 * \param nlemt: `int` QUESTION: definition?
-	 */
-	C9(int ndi, int nlem, int ndit, int nlemt);
+  /*! \brief C9 instance constructor.
+   *
+   * \param ndi: `int` QUESTION: definition?
+   * \param nlem: `int` QUESTION: definition?
+   * \param ndit: `int` QUESTION: definition?
+   * \param nlemt: `int` QUESTION: definition?
+   */
+  C9(int ndi, int nlem, int ndit, int nlemt);
 
-	/*! \brief C9 instance destroyer.
-	 */
-	~C9();
+  /*! \brief C9 instance destroyer.
+   */
+  ~C9();
 };
 #endif
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index fc9498a8..c88ff579 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -25,25 +25,29 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
   nsph = ns;
   lm = l_max;
 
+  vec_rmi = new dcomplex[nsph * lm]();
+  vec_rei = new dcomplex[nsph * lm]();
   rmi = new dcomplex*[lm];
   rei = new dcomplex*[lm];
   for (int ri = 0; ri < lm; ri++) {
-    rmi[ri] = new dcomplex[nsph]();
-    rei[ri] = new dcomplex[nsph]();
+    rmi[ri] = &(vec_rmi[nsph * ri]);
+    rei[ri] = &(vec_rei[nsph * ri]);
   }
+  vec_w = new dcomplex[4 * nlmmt]();
   w = new dcomplex*[nlmmt];
-  for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[4]();
+  for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
   int configurations = 0;
   for (int ci = 1; ci <= nsph; ci++) {
     if (_iog[ci - 1] >= ci) configurations++;
   }
+  vec_vints = new dcomplex[nsph * 16]();
   vints = new dcomplex*[nsph];
   rc = new double*[configurations];
   nshl = new int[configurations]();
   iog = new int[nsph]();
   int conf_index = 0;
   for (int vi = 0; vi < nsph; vi++) {
-    vints[vi] = new dcomplex[16]();
+    vints[vi] = &(vec_vints[vi * 16]);
     iog[vi] = _iog[vi];
     if (iog[vi] >= vi + 1) {
       nshl[conf_index] = _nshl[conf_index];
@@ -73,9 +77,12 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
 }
 
 C1::~C1() {
+  delete[] vec_rmi;
+  delete[] vec_rei;
   delete[] rmi;
   delete[] rei;
-  for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
+  delete[] vec_w;
+  delete[] w;
   int conf_index = 0;
   for (int ci = 1; ci <= nsph; ci++) {
     if (iog[ci] >= ci) {
@@ -84,9 +91,7 @@ C1::~C1() {
     }
   }
   delete[] rc;
-  for (int vi = nsph - 1; vi > - 1; vi--) {
-    delete[] vints[vi];
-  }
+  delete[] vec_vints;
   delete[] vints;
   for (int si = nsph - 1; si > -1; si--) {
     delete[] sas[si][1];
-- 
GitLab