Introduce data structures for the cluster case

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

src/libnptm/Commons.cpp

 /*! \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;
+}