From 4f5ead92a9e6d39ce725559cac3119b3afdf5769 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Wed, 8 Nov 2023 14:58:34 +0100
Subject: [PATCH] Introduce data structures for the cluster case
---
src/include/Commons.h | 175 ++++++++++++++++++++++++++++++++++++-
src/libnptm/Commons.cpp | 185 +++++++++++++++++++++++++++++++++++++---
2 files changed, 345 insertions(+), 15 deletions(-)
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 088dcc74..eef2cdaa 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -7,15 +7,15 @@
* the data blocks in the code, therefore implying the necessity to modify
* the code to adapt it to the input and to recompile before running. C++,
* on the contrary, offers the possibility to represent the necessary data
- * structures as classes that can dinamically instantiate the shared information
+ * structures as classes that can dynamically instantiate the shared information
* in the most convenient format for the current configuration. This approach
* adds an abstraction layer that lifts the need to modify and recompile the
* code depending on the input.
*
*/
-#ifndef SRC_INCLUDE_COMMONS_
-#define SRC_INCLUDE_COMMONS_
+#ifndef INCLUDE_COMMONS_
+#define INCLUDE_COMMONS_
#include <complex>
@@ -89,7 +89,7 @@ public:
* \param ns: `int` Number of spheres.
* \param l_max: `int` Maximum order of field expansion.
*/
- C1(int ns, int l_max);
+ C1(int ns, int l_max, int *nshl, int *iog);
//! \brief C1 instance destroyer.
~C1();
@@ -128,4 +128,171 @@ public:
~C2();
};
+/*! \brief Representation of the FORTRAN C3 blocks.
+ */
+class C3 {
+public:
+ //! \brief QUESTION: definition?
+ std::complex<double> tfsas;
+ //! \brief QUESTION: definition?
+ std::complex<double> **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 destroyer.
+ */
+ ~C3();
+};
+
+/*! \brief Representation of the FORTRAN C4 blocks.
+ */
+struct C4 {
+ //! \brief QUESTION: definition?
+ int litpo;
+ //! \brief QUESTION: definition?
+ int litpos;
+ //! \brief QUESTION: definition?
+ int lmpo;
+ //! \brief QUESTION: definition?
+ int lmtpo;
+ //! \brief QUESTION: definition?
+ 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 Vectors and matrices that are specific to cluster C1 blocks.
+ *
+ */
+class C1_AddOns {
+protected:
+ //! \brief Number of spheres.
+ int nsph;
+ //! \brief QUESTION: definition?
+ int nlemt;
+ //! \brief QUESTION: definition?
+ int lmpo;
+
+ void allocate_vectors(C4 *c4);
+public:
+ //! \brief QUESTION: definition?
+ std::complex<double> *vh;
+ //! \brief QUESTION: definition?
+ std::complex<double> *vj;
+ //! \brief QUESTION: definition?
+ std::complex<double> *vyhj;
+ //! \brief QUESTION: definition?
+ std::complex<double> *vyj0;
+ //! \brief QUESTION: definition?
+ std::complex<double> **am0m;
+ //! \brief QUESTION: definition?
+ std::complex<double> *vint;
+ //! \brief QUESTION: definition?
+ std::complex<double> *vintm;
+ //! \brief QUESTION: definition?
+ std::complex<double> *vintt;
+ //! \brief QUESTION: definition?
+ std::complex<double> **fsac;
+ //! \brief QUESTION: definition?
+ std::complex<double> **sac;
+ //! \brief QUESTION: definition?
+ std::complex<double> **fsacm;
+ //! \brief QUESTION: definition?
+ std::complex<double> *scscp;
+ //! \brief QUESTION: definition?
+ std::complex<double> *ecscp;
+ //! \brief QUESTION: definition?
+ std::complex<double> *scscpm;
+ //! \brief QUESTION: definition?
+ std::complex<double> *ecscpm;
+ //! \brief QUESTION: definition?
+ double *v3j0;
+ //! \brief QUESTION: definition?
+ double *sscs;
+ //! \brief QUESTION: definition?
+ int **ind3j;
+
+ /*! \brief C1_AddOns instance constructor.
+ *
+ * \param ns: `int` Number of spheres.
+ * \param nc: `int`
+ * \param litpo: `int`
+ * \param lmtpo: `int`
+ * \param nv3j: `int`
+ * \param lmpo: `int`
+ * \param li: `int`
+ * \param le: `int`
+ * \param lm: `int`
+ */
+ C1_AddOns(C4 *c4);
+
+ //! \brief C1_AddOns instance destroyer.
+ ~C1_AddOns();
+};
+
+/*! \brief Representation of the FORTRAN C6 blocks.
+ */
+class C6 {
+public:
+ //! \brief QUESTION: definition?
+ double *rac3j;
+
+ /*! \brief C6 instance constructor.
+ *
+ * \param lmtpo: `int` QUESTION: definition?
+ */
+ C6(int lmtpo);
+
+ /*! \brief C6 instance destroyer.
+ */
+ ~C6();
+};
+
+/*! \brief Representation of the FORTRAN C9 blocks.
+ */
+class C9 {
+protected:
+ int gis_size_0;
+ int sam_size_0;
+public:
+ //! \brief QUESTION: definition?
+ std::complex<double> **gis;
+ //! \brief QUESTION: definition?
+ std::complex<double> **gls;
+ //! \brief QUESTION: definition?
+ std::complex<double> **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 destroyer.
+ */
+ ~C9();
+};
#endif
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 809337ea..2140c69a 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1,12 +1,23 @@
/*! \file Commons.cpp
*
+ * DEVELOPMENT NOTE:
+ * The construction of common blocks requires some information
+ * that is stored in configuration objects and is needed to compute
+ * the allocation size of vectors and matrices. Currently, this
+ * information is passed as arguments to the constructors of the
+ * common blocks. A simpler and more logical way to operate is
+ * to design the constructors to take as arguments only pointers
+ * to the configuration objects. These, on their turn, need to
+ * expose methods to access the relevant data in read-only mode.
*/
+#ifndef INCLUDE_COMMONS_H
#include "../include/Commons.h"
+#endif
using namespace std;
-C1::C1(int ns, int l_max) {
+C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
nlmmt = 2 * (l_max * (l_max + 2));
nsph = ns;
lm = l_max;
@@ -21,9 +32,13 @@ C1::C1(int ns, int l_max) {
for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4];
vints = new complex<double>*[nsph];
rc = new double*[nsph];
+ nshl = new int[nsph];
+ iog = new int[nsph];
for (int vi = 0; vi < nsph; vi++) {
- rc[vi] = new double[lm];
+ rc[vi] = new double[_nshl[vi]];
vints[vi] = new complex<double>[16];
+ nshl[vi] = _nshl[vi];
+ iog[vi] = _iog[vi];
}
fsas = new complex<double>[nsph];
sscs = new double[nsph];
@@ -37,8 +52,6 @@ C1::C1(int ns, int l_max) {
ryy = new double[nsph];
rzz = new double[nsph];
ros = new double[nsph];
- iog = new int[nsph];
- nshl = new int[nsph];
sas = new complex<double>**[nsph];
for (int si = 0; si < nsph; si++) {
@@ -51,15 +64,19 @@ C1::C1(int ns, int l_max) {
C1::~C1() {
delete[] rmi;
delete[] rei;
- for (int wi = 1; wi <= nlmmt; wi++) delete[] w[wi];
- for (int vi = 1; vi <= nsph; vi++) {
- delete[] rc[nsph - vi];
- delete[] vints[nsph - vi];
+ for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
+ for (int vi = nsph - 1; vi > - 1; vi--) {
+ delete[] rc[vi];
+ delete[] vints[vi];
}
- for (int si = 1; si <= nsph; si++) {
- delete[] sas[nsph - si][1];
- delete[] sas[nsph - si][0];
+ delete[] rc;
+ delete[] vints;
+ for (int si = nsph - 1; si > -1; si--) {
+ delete[] sas[si][1];
+ delete[] sas[si][0];
+ delete[] sas[si];
}
+ delete[] sas;
delete[] fsas;
delete[] sscs;
delete[] sexs;
@@ -76,6 +93,103 @@ C1::~C1() {
delete[] nshl;
}
+C1_AddOns::C1_AddOns(C4 *c4) {
+ nsph = c4->nsph;
+ lmpo = c4->lmpo;
+ nlemt = 2 * c4->nlem;
+ vh = new complex<double>[(nsph * nsph - 1) * c4->litpo];
+ vj = new complex<double>[nsph * c4->lmtpo];
+ vyhj = new complex<double>[nsph * (nsph - 1) * c4->litpo];
+ vyj0 = new complex<double>[nsph * c4->lmtpo];
+ am0m = new complex<double>*[nlemt];
+ for (int ai = 0; ai < nlemt; ai++) {
+ am0m[ai] = new complex<double>[nlemt];
+ }
+ vint = new complex<double>[16];
+ vintm = new complex<double>[16];
+ vintt = new complex<double>[16];
+ fsac = new complex<double>*[2];
+ sac = new complex<double>*[2];
+ fsacm = new complex<double>*[2];
+ for (int fi = 0; fi < 2; fi++) {
+ fsac[fi] = new complex<double>[2];
+ sac[fi] = new complex<double>[2];
+ fsacm[fi] = new complex<double>[2];
+ }
+ scscp = new complex<double>[2];
+ ecscp = new complex<double>[2];
+ scscpm = new complex<double>[2];
+ ecscpm = new complex<double>[2];
+ allocate_vectors(c4);
+ sscs = new double[nsph];
+}
+
+C1_AddOns::~C1_AddOns() {
+ delete[] sscs;
+ delete[] vyj0;
+ delete[] vyhj;
+ for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
+ delete[] ind3j;
+ delete[] v3j0;
+ delete[] vh;
+ delete[] vj;
+ for (int ai = nlemt - 1; ai > -1; ai--) {
+ delete[] am0m[ai];
+ }
+ delete am0m;
+ delete[] vint;
+ delete[] vintm;
+ delete[] vintt;
+ for (int fi = 1; fi > -1; fi--) {
+ delete[] fsac[fi];
+ delete[] sac[fi];
+ delete[] fsacm[fi];
+ }
+ delete[] fsac;
+ delete[] sac;
+ delete[] fsacm;
+ delete[] scscp;
+ delete[] ecscp;
+ delete[] scscpm;
+ delete[] ecscpm;
+}
+
+void C1_AddOns::allocate_vectors(C4 *c4) {
+ int i = 0;
+ int l1po_max = 0, l2_max = 0;
+ int i_max = 0;
+ // This loop calculates NV3J, the required size of v3j0
+ for (int l1po = 1; l1po <= c4->lmpo; l1po++) {
+ int l1 = l1po - 1;
+ for (int l2 = 1; l2 <= c4->lm; l2++) {
+ if (l1po > l1po_max) l1po_max = l1po;
+ if (l2 > l2_max) l2_max = l2;
+ int iabs = l2 - l1;
+ if (iabs < 0) iabs *= -1;
+ int lmnpo = iabs + 1;
+ int lmxpo = l2 + l1 + 1;
+ int il = 0;
+ int lpo = lmnpo;
+ while (lpo <= lmxpo) {
+ i++;
+ il++;
+ if (i > i_max) i_max = i;
+ lpo += 2;
+ }
+ }
+ } // l1po loop
+ v3j0 = new double[i_max];
+ ind3j = new int*[lmpo];
+ for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm];
+ // This calculates the size of vyhj
+ int nsphmo = c4->nsph - 1;
+ int ivy = nsphmo * nsph * c4->litpos;
+ vyhj = new complex<double>[ivy];
+ // This calculates the size of vyj0
+ ivy = c4->lmtpos * c4->nsph;
+ vyj0 = new complex<double>[ivy];
+}
+
C2::C2(int ns, int nl, int npnt, int npntts) {
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
@@ -94,3 +208,52 @@ C2::~C2() {
delete[] dc0;
delete[] vsz;
}
+
+C3::C3() {
+ tsas = new complex<double>*[2];
+ tsas[0] = new complex<double>[2];
+ tsas[1] = new complex<double>[2];
+ tfsas = complex<double>(0.0, 0.0);
+ gcs = 0.0;
+ scs = 0.0;
+ ecs = 0.0;
+ acs = 0.0;
+}
+
+C3::~C3() {
+ delete[] tsas[1];
+ delete[] tsas[0];
+ delete[] tsas;
+}
+
+C6::C6(int lmtpo) {
+ rac3j = new double[lmtpo];
+}
+
+C6::~C6() {
+ delete[] rac3j;
+}
+
+C9::C9(int ndi, int nlem, int ndit, int nlemt) {
+ gis = new complex<double>*[ndi];
+ gls = new complex<double>*[ndi];
+ for (int gi = 0; gi < ndi; gi++) {
+ gis[gi] = new complex<double>[nlem];
+ gls[gi] = new complex<double>[nlem];
+ }
+ sam = new complex<double>*[ndit];
+ for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt];
+ gis_size_0 = ndi;
+ sam_size_0 = ndit;
+}
+
+C9::~C9() {
+ for (int gi = gis_size_0 - 1; gi > -1; gi--) {
+ delete[] gis[gi];
+ delete[] gls[gi];
+ }
+ delete[] gis;
+ delete[] gls;
+ for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
+ delete[] sam;
+}
--
GitLab