diff --git a/src/Commons.cpp b/src/Commons.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d5141e077a825428ee0bca173705f4a96caa1f0b --- /dev/null +++ b/src/Commons.cpp @@ -0,0 +1,96 @@ +/*! \file Commons.cpp + * + */ + +#include "include/Commons.h" + +using namespace std; + +C1::C1(int ns, int l_max) { + nlmmt = 2 * (l_max * (l_max + 2)); + nsph = ns; + lm = l_max; + + rmi = new complex<double>*[lm]; + rei = new complex<double>*[lm]; + for (int ri = 0; ri < lm; ri++) { + rmi[ri] = new complex<double>[nsph]; + rei[ri] = new complex<double>[nsph]; + } + w = new complex<double>*[nlmmt]; + for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]; + vints = new complex<double>*[nsph]; + rc = new double*[nsph]; + for (int vi = 0; vi < nsph; vi++) { + rc[vi] = new double[lm]; + vints[vi] = new complex<double>[16]; + } + fsas = new complex<double>[nsph]; + sscs = new double[nsph]; + sexs = new double[nsph]; + sabs = new double[nsph]; + sqscs = new double[nsph]; + sqexs = new double[nsph]; + sqabs = new double[nsph]; + gcsv = new double[nsph];; + rxx = new double[nsph]; + 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++) { + sas[si] = new complex<double>*[2]; + sas[si][0] = new complex<double>[2]; + sas[si][1] = new complex<double>[2]; + } +} + +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 si = 1; si <= nsph; si++) { + delete[] sas[nsph - si][1]; + delete[] sas[nsph - si][0]; + } + delete[] fsas; + delete[] sscs; + delete[] sexs; + delete[] sabs; + delete[] sqscs; + delete[] sqexs; + delete[] sqabs; + delete[] gcsv; + delete[] rxx; + delete[] ryy; + delete[] rzz; + delete[] ros; + delete[] iog; + delete[] nshl; +} + +C2::C2(int ns, int nl, int npnt, int npntts) { + nsph = ns; + int max_n = (npnt > npntts) ? npnt : npntts; + nhspo = 2 * max_n - 1; + ris = new complex<double>[nhspo]; + dlri = new complex<double>[nhspo]; + vkt = new complex<double>[nsph]; + dc0 = new complex<double>[nl]; + vsz = new double[nsph]; +} + +C2::~C2() { + delete[] ris; + delete[] dlri; + delete[] vkt; + delete[] dc0; + delete[] vsz; +} diff --git a/src/Makefile b/src/Makefile index 9e7b374bd0cba51d70cc207467d4b1d145ce686e..8e30f71de441b29d745ba944206be9cb71290550 100644 --- a/src/Makefile +++ b/src/Makefile @@ -8,8 +8,8 @@ all: $(SUBDIRS) $(SUBDIRS): $(MAKE) -C $@ -np_sphere: $(BUILLDDIR)/sphere/np_sphere.o $(BUILLDDIR)/sphere/Configuration.o $(BUILLDDIR)/sphere/Parsers.o $(BUILLDDIR)/sphere/sphere.o - $(CC) -o $(BUILDDIR)/sphere/np_sphere $(BUILDDIR)/sphere/np_sphere.o $(BUILDDIR)/sphere/Configuration.o $(BUILDDIR)/sphere/Parsers.o $(BUILDDIR)/sphere/sphere.o +np_sphere: $(BUILDDIR)/sphere/np_sphere.o $(BUILDDIR)/sphere/Commons.o $(BUILDDIR)/sphere/Configuration.o $(BUILDDIR)/sphere/Parsers.o $(BUILDDIR)/sphere/sphere.o + $(CC) -o $(BUILDDIR)/sphere/np_sphere $(BUILDDIR)/sphere/np_sphere.o $(BUILDDIR)/sphere/Commons.o $(BUILDDIR)/sphere/Configuration.o $(BUILDDIR)/sphere/Parsers.o $(BUILDDIR)/sphere/sphere.o clean: rm -f $(BUILDDIR)/cluster/*.o @@ -23,14 +23,17 @@ wipe: .PHONY: all $(SUBDIRS) -$(BUILLDDIR)/sphere/np_sphere.o: +$(BUILDDIR)/sphere/np_sphere.o: $(CC) -c np_sphere.cpp -o $(BUILDDIR)/sphere/np_sphere.o -$(BUILLDDIR)/sphere/Configuration.o: +$(BUILDDIR)/sphere/Commons.o: + $(CC) -c Commons.cpp -o $(BUILDDIR)/sphere/Commons.o + +$(BUILDDIR)/sphere/Configuration.o: $(CC) -c Configuration.cpp -o $(BUILDDIR)/sphere/Configuration.o -$(BUILLDDIR)/sphere/Parsers.o: +$(BUILDDIR)/sphere/Parsers.o: $(CC) -c Parsers.cpp -o $(BUILDDIR)/sphere/Parsers.o -$(BUILLDDIR)/sphere/sphere.o: +$(BUILDDIR)/sphere/sphere.o: $(CC) -c sphere/sphere.cpp -o $(BUILDDIR)/sphere/sphere.o diff --git a/src/include/Commons.h b/src/include/Commons.h new file mode 100644 index 0000000000000000000000000000000000000000..3d176820ea8b152f5a65ee6fd3a8ba82b52f8ab8 --- /dev/null +++ b/src/include/Commons.h @@ -0,0 +1,103 @@ +/*! \file Commons.h + * + * \brief C++ porting of common data structures. + * + * Many functions of the original FORTRAN code share complex data blocks in + * form of COMMON blocks. This poses the limit of freezing the structure of + * 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 + * 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_ + +#include <complex> + +/*! \brief Representation of the FORTRAN C1 common blocks. + * + * C1 common blocks are used to store vector field expansions and geometric + * properties, such as sphere sizes, positions and cross-sections. These are + * used by functions such as `aps`, `diel`, `pwma`, `rabas`, `sscr0`, `sscr2`, + * `wmamp`, `wmasp` and `dme`. QUESTION: correct? + * + * Due to the necessity to share the class contents with many functions that + * need to read and update various fields, all shared members of C1 are declared + * public (i.e., the class is just a structure with more advanced constructor + * and destroyer). Further development may go in the direction of creating + * a better encapsulation, either by unpacking the contents (recommended) or + * by creating public methods to access protected fields (standard way, but not + * recommended for HPC applications). + */ +class C1 { +protected: + //! \brief Number of spheres. + int nsph; + //! \brief Maximum order of field expansion. + int lm; + //! \brief NLMMT. QUESTION: definition? + int nlmmt; +public: + std::complex<double> **rmi; + std::complex<double> **rei; + std::complex<double> **w; + std::complex<double> *fsas; + std::complex<double> **vints; + double *sscs; + double *sexs; + double *sabs; + double *sqscs; + double *sqexs; + double *sqabs; + double *gcsv; + double *rxx; + double *ryy; + double *rzz; + double *ros; + double **rc; + int *iog; + int *nshl; + std::complex<double> ***sas; + + /*! \brief C1 instance constructor. + * + * \param ns: `int` Number of spheres. + * \param l_max: `int` Maximum order of field expansion. + */ + C1(int ns, int l_max); + + //! \brief C1 instance destroyer. + ~C1(); +}; + +/*! \brief Representation of the FORTRAN C2 blocks. + * + */ +class C2 { + int nsph, nhspo; +public: + std::complex<double> *ris; + std::complex<double> *dlri; + std::complex<double> *dc0; + std::complex<double> *vkt; + 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 destroyer. + ~C2(); +}; + +#endif diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h index fffd1085e37d8ed6b3a0606efd3b6be29e58a1ac..d379701bdb8bcc4faf844a33b64bbcd5bab45c03 100644 --- a/src/include/sph_subs.h +++ b/src/include/sph_subs.h @@ -15,39 +15,7 @@ #define SRC_INCLUDE_SPH_SUBS_H_ #include <complex> - -//! \brief A structure representing the common block C1 -struct C1 { - std::complex<double> rmi[40][2]; - std::complex<double> rei[40][2]; - std::complex<double> w[3360][4]; - std::complex<double> fsas[2]; - std::complex<double> sas[2][2][2]; - std::complex<double> vints[2][16]; - double sscs[2]; - double sexs[2]; - double sabs[2]; - double sqscs[2]; - double sqexs[2]; - double sqabs[2]; - double gcsv[2]; - double rxx[1]; - double ryy[1]; - double rzz[1]; - double ros[2]; - double rc[2][8]; - int iog[2]; - int nshl[2]; -}; - -//! \brief A structure representing the common block C1 -struct C2 { - std::complex<double> ris[999]; - std::complex<double> dlri[999]; - std::complex<double> dc0[5]; - std::complex<double> vkt[2]; - double vsz[2]; -}; +#include "Commons.h" /*! \brief Conjugate of a double precision complex number * @@ -776,7 +744,7 @@ void rkt( yy = y1 + dy1 * step; c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step; y1= yy + (c11 + c12 + c13) * step / 6.0; - dy1 += (0.5 * c11 + c12 + c13 + 0.5 *c14) /3.0; + dy1 += (0.5 * c11 + c12 + c13 + 0.5 * c14) /3.0; cy1 -= cdy1 * c2->dlri[jpnt - 1]; cdy1 += 2.0 * c2->dlri[jpnt - 1]; c21 = (cy1 * y2 + cdy1 * dy2) * step; diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp index cf037183e87dc5babdc5842ee8c77ed88533fe45..5d92b19255106177f5dc30201704280ecc4a3e90 100644 --- a/src/sphere/sphere.cpp +++ b/src/sphere/sphere.cpp @@ -49,7 +49,7 @@ void sphere() { complex<double> *vint = new complex<double>[16]; int jw; int nsph = gconf->number_of_spheres; - C1 *c1 = new C1; + C1 *c1 = new C1(nsph, gconf->l_max); for (int i = 0; i < nsph; i++) { c1->iog[i] = sconf->iog_vec[i]; c1->nshl[i] = sconf->nshl_vec[i]; @@ -61,7 +61,7 @@ void sphere() { c1->rei[i][0] = complex<double>(0.0, 0.0); c1->rei[i][1] = complex<double>(0.0, 0.0); } - C2 *c2 = new C2; + C2 *c2 = new C2(nsph, 5, gconf->npnt, gconf->npntts); argi = new double[1]; args = new double[1]; gaps = new double[2];