diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 7d2533eac1e1f313fcd14b29e8805c0aeb229534..8c499689ef27676671b5d05f2a9ddd3d0d117b05 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -95,9 +95,11 @@
#include "../include/outputs.h"
#endif
-using namespace std;
+#ifndef INCLUDE_ITERATION_DATA_H_
+#include "../include/IterationData.h"
+#endif
-// 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...
+using namespace std;
/*! \brief Main calculation loop.
*
@@ -1531,3 +1533,536 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
return jer;
}
+
+// >>> IMPLEMENTATION OF ClusterIterationData CLASS <<<
+ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count) {
+ c1 = new ParticleDescriptorCluster(gconf, sconf);
+ const int ndi = c1->nsph * c1->nlim;
+ const np_int ndit = 2 * ndi;
+ gaps = new double[c1->nsph]();
+ tqev = new double[3]();
+ tqsv = new double[3]();
+ tqse = new double*[2];
+ tqspe = new dcomplex*[2];
+ tqss = new double*[2];
+ tqsps = new dcomplex*[2];
+ tqce = new double*[2];
+ tqcpe = new dcomplex*[2];
+ tqcs = new double*[2];
+ tqcps = new dcomplex*[2];
+ for (int ti = 0; ti < 2; ti++) {
+ tqse[ti] = new double[c1->nsph]();
+ tqspe[ti] = new dcomplex[c1->nsph]();
+ tqss[ti] = new double[c1->nsph]();
+ tqsps[ti] = new dcomplex[c1->nsph]();
+ tqce[ti] = new double[3]();
+ tqcpe[ti] = new dcomplex[3]();
+ tqcs[ti] = new double[3]();
+ tqcps[ti] = new dcomplex[3]();
+ }
+ gapv = new double[3]();
+ gapp = new dcomplex*[3];
+ gappm = new dcomplex*[3];
+ gap = new double*[3];
+ gapm = new double*[3];
+ for (int gi = 0; gi < 3; gi++) {
+ gapp[gi] = new dcomplex[2]();
+ gappm[gi] = new dcomplex[2]();
+ gap[gi] = new double[2]();
+ gapm[gi] = new double[2]();
+ }
+ u = new double[3]();
+ us = new double[3]();
+ un = new double[3]();
+ uns = new double[3]();
+ up = new double[3]();
+ ups = new double[3]();
+ unmp = new double[3]();
+ unsmp = new double[3]();
+ upmp = new double[3]();
+ upsmp = new double[3]();
+ argi = new double[1]();
+ args = new double[1]();
+ duk = new double[3]();
+ cextlr = new double*[4];
+ cext = new double*[4];
+ cmullr = new double*[4];;
+ cmul = new double*[4];
+ for (int ci = 0; ci < 4; ci++) {
+ cextlr[ci] = new double[4]();
+ cext[ci] = new double[4]();
+ cmullr[ci] = new double[4]();
+ cmul[ci] = new double[4]();
+ }
+ zpv = new double***[c1->lm];
+ for (int zi = 0; zi < c1->lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ for (int zk = 0; zk < 2; zk++) {
+ zpv[zi][zj][zk] = new double[2]();
+ }
+ }
+ }
+ am_vector = new dcomplex[ndit * ndit]();
+ am = new dcomplex*[ndit];
+ for (int ai = 0; ai < ndit; ai++) {
+ am[ai] = (am_vector + ai * ndit);
+ }
+
+ arg = 0.0 + 0.0 * I;
+ // These are suspect initializations
+ scan = 0.0;
+ cfmp = 0.0;
+ sfmp = 0.0;
+ cfsp = 0.0;
+ sfsp = 0.0;
+ // End of suspect initializations
+ wn = sconf->wp / 3.0e8;
+ xip = sconf->xip;
+ sqsfi = 1.0;
+ vk = 0.0;
+ number_of_scales = sconf->number_of_scales;
+ xiblock = (int) ceil(((double) (sconf->number_of_scales-1))/((double) mpidata->nprocs));
+ lastxi = ((mpidata->rank+1) * xiblock)+1;
+ firstxi = lastxi-xiblock+1;
+ if (lastxi > sconf->number_of_scales) lastxi = sconf->number_of_scales;
+
+#ifdef USE_MAGMA
+ proc_device = mpidata->rank % device_count;
+#else
+ proc_device = 0;
+#endif
+
+ // In the first iteration, if refinement is enabled, determine the number of refinement iterations required to arrive at the target accuracy (if achievable in a reasonable number of iterations)
+ refinemode = 2;
+ // maxrefiters and accuracygoal should be configurable and preferably set somewhere else
+ maxrefiters = 20;
+ accuracygoal = 1e-6;
+}
+
+ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
+ c1 = new ParticleDescriptorCluster(reinterpret_cast<ParticleDescriptorCluster &>(*(rhs.c1)));
+ const int ndi = c1->nsph * c1->nlim;
+ const np_int ndit = 2 * ndi;
+ gaps = new double[c1->nsph]();
+ for (int gi = 0; gi < c1->nsph; gi++) gaps[gi] = rhs.gaps[gi];
+ tqev = new double[3]();
+ tqsv = new double[3]();
+ for (int ti = 0; ti < 3; ti++) {
+ tqev[ti] = rhs.tqev[ti];
+ tqsv[ti] = rhs.tqsv[ti];
+ }
+ tqse = new double*[2];
+ tqspe = new dcomplex*[2];
+ tqss = new double*[2];
+ tqsps = new dcomplex*[2];
+ tqce = new double*[2];
+ tqcpe = new dcomplex*[2];
+ tqcs = new double*[2];
+ tqcps = new dcomplex*[2];
+ for (int ti = 0; ti < 2; ti++) {
+ tqse[ti] = new double[c1->nsph]();
+ tqspe[ti] = new dcomplex[c1->nsph]();
+ tqss[ti] = new double[c1->nsph]();
+ tqsps[ti] = new dcomplex[c1->nsph]();
+ for (int tj = 0; tj < c1->nsph; tj++) {
+ tqse[ti][tj] = rhs.tqse[ti][tj];
+ tqspe[ti][tj] = rhs.tqspe[ti][tj];
+ tqss[ti][tj] = rhs.tqss[ti][tj];
+ tqsps[ti][tj] = rhs.tqsps[ti][tj];
+ }
+ tqce[ti] = new double[3]();
+ tqcpe[ti] = new dcomplex[3]();
+ tqcs[ti] = new double[3]();
+ tqcps[ti] = new dcomplex[3]();
+ for (int tj = 0; tj < 3; tj++) {
+ tqce[ti][tj] = rhs.tqce[ti][tj];
+ tqcpe[ti][tj] = rhs.tqcpe[ti][tj];
+ tqcs[ti][tj] = rhs.tqcs[ti][tj];
+ tqcps[ti][tj] = rhs.tqcps[ti][tj];
+ }
+ }
+ gapv = new double[3]();
+ gapp = new dcomplex*[3];
+ gappm = new dcomplex*[3];
+ gap = new double*[3];
+ gapm = new double*[3];
+ for (int gi = 0; gi < 3; gi++) {
+ gapv[gi] = rhs.gapv[gi];
+ gapp[gi] = new dcomplex[2]();
+ gappm[gi] = new dcomplex[2]();
+ gap[gi] = new double[2]();
+ gapm[gi] = new double[2]();
+ for (int gj = 0; gj < 2; gj++) {
+ gapp[gi][gj] = rhs.gapp[gi][gj];
+ gappm[gi][gj] = rhs.gappm[gi][gj];
+ gap[gi][gj] = rhs.gap[gi][gj];
+ gapm[gi][gj] = rhs.gapm[gi][gj];
+ }
+ }
+ u = new double[3]();
+ us = new double[3]();
+ un = new double[3]();
+ uns = new double[3]();
+ up = new double[3]();
+ ups = new double[3]();
+ unmp = new double[3]();
+ unsmp = new double[3]();
+ upmp = new double[3]();
+ upsmp = new double[3]();
+ duk = new double[3]();
+ for (int ui = 0; ui < 3; ui++) {
+ u[ui] = rhs.u[ui];
+ us[ui] = rhs.us[ui];
+ un[ui] = rhs.un[ui];
+ uns[ui] = rhs.uns[ui];
+ up[ui] = rhs.up[ui];
+ ups[ui] = rhs.ups[ui];
+ unmp[ui] = rhs.unmp[ui];
+ unsmp[ui] = rhs.unsmp[ui];
+ upmp[ui] = rhs.upmp[ui];
+ upsmp[ui] = rhs.upsmp[ui];
+ duk[ui] = rhs.duk[ui];
+ }
+ argi = new double[1]();
+ args = new double[1]();
+ argi[0] = rhs.argi[0];
+ args[0] = rhs.args[0];
+ cextlr = new double*[4];
+ cext = new double*[4];
+ cmullr = new double*[4];;
+ cmul = new double*[4];
+ for (int ci = 0; ci < 4; ci++) {
+ cextlr[ci] = new double[4]();
+ cext[ci] = new double[4]();
+ cmullr[ci] = new double[4]();
+ cmul[ci] = new double[4]();
+ for (int cj = 0; cj < 4; cj++) {
+ cextlr[ci][cj] = rhs.cextlr[ci][cj];
+ cext[ci][cj] = rhs.cext[ci][cj];
+ cmullr[ci][cj] = rhs.cmullr[ci][cj];
+ cmul[ci][cj] = rhs.cmul[ci][cj];
+ }
+ }
+ zpv = new double***[c1->lm];
+ for (int zi = 0; zi < c1->lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ for (int zk = 0; zk < 2; zk++) {
+ zpv[zi][zj][zk] = new double[2]();
+ zpv[zi][zj][zk][0] = rhs.zpv[zi][zj][zk][0];
+ zpv[zi][zj][zk][1] = rhs.zpv[zi][zj][zk][1];
+ }
+ }
+ }
+ am_vector = new dcomplex[ndit * ndit]();
+ for (np_int ai = 0; ai < ndit * ndit; ai++) am_vector[ai] = rhs.am_vector[ai];
+ am = new dcomplex*[ndit];
+ for (np_int ai = 0; ai < ndit; ai++) {
+ am[ai] = (am_vector + ai * ndit);
+ }
+
+ arg = rhs.arg;
+ // These are suspect initializations
+ scan = rhs.scan;
+ cfmp = rhs.cfmp;
+ sfmp = rhs.sfmp;
+ cfsp = rhs.cfsp;
+ sfsp = rhs.sfsp;
+ // End of suspect initializations
+ wn = rhs.wn;
+ xip = rhs.xip;
+ sqsfi = rhs.sqsfi;
+ vk = rhs.vk;
+ firstxi = rhs.firstxi;
+ lastxi = rhs.lastxi;
+ xiblock = rhs.xiblock;
+ number_of_scales = rhs.number_of_scales;
+
+ proc_device = rhs.proc_device;
+ refinemode = rhs.refinemode;
+ maxrefiters = rhs.maxrefiters;
+ accuracygoal = rhs.accuracygoal;
+}
+
+#ifdef MPI_VERSION
+ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int device_count) {
+ c1 = new ParticleDescriptorCluster(mpidata);
+ const int ndi = c1->nsph * c1->nlim;
+ const np_int ndit = 2 * ndi;
+ gaps = new double[c1->nsph]();
+ MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ tqev = new double[3]();
+ tqsv = new double[3]();
+ MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ tqse = new double*[2];
+ tqspe = new dcomplex*[2];
+ tqss = new double*[2];
+ tqsps = new dcomplex*[2];
+ tqce = new double*[2];
+ tqcpe = new dcomplex*[2];
+ tqcs = new double*[2];
+ tqcps = new dcomplex*[2];
+ for (int ti = 0; ti < 2; ti++) {
+ tqse[ti] = new double[c1->nsph]();
+ tqspe[ti] = new dcomplex[c1->nsph]();
+ tqss[ti] = new double[c1->nsph]();
+ tqsps[ti] = new dcomplex[c1->nsph]();
+ MPI_Bcast(tqse[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqspe[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqss[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqsps[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ tqce[ti] = new double[3]();
+ tqcpe[ti] = new dcomplex[3]();
+ tqcs[ti] = new double[3]();
+ tqcps[ti] = new dcomplex[3]();
+ MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ }
+ gapv = new double[3]();
+ gapp = new dcomplex*[3];
+ gappm = new dcomplex*[3];
+ gap = new double*[3];
+ gapm = new double*[3];
+ MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ for (int gi = 0; gi < 3; gi++) {
+ gapp[gi] = new dcomplex[2]();
+ gappm[gi] = new dcomplex[2]();
+ gap[gi] = new double[2]();
+ gapm[gi] = new double[2]();
+ MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ }
+ u = new double[3]();
+ us = new double[3]();
+ un = new double[3]();
+ uns = new double[3]();
+ up = new double[3]();
+ ups = new double[3]();
+ unmp = new double[3]();
+ unsmp = new double[3]();
+ upmp = new double[3]();
+ upsmp = new double[3]();
+ duk = new double[3]();
+ MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ argi = new double[1]();
+ args = new double[1]();
+ MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ cextlr = new double*[4];
+ cext = new double*[4];
+ cmullr = new double*[4];;
+ cmul = new double*[4];
+ for (int ci = 0; ci < 4; ci++) {
+ cextlr[ci] = new double[4]();
+ cext[ci] = new double[4]();
+ cmullr[ci] = new double[4]();
+ cmul[ci] = new double[4]();
+ MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ }
+ zpv = new double***[c1->lm];
+ for (int zi = 0; zi < c1->lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ for (int zk = 0; zk < 2; zk++) {
+ zpv[zi][zj][zk] = new double[2]();
+ MPI_Bcast(zpv[zi][zj][zk], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ }
+ }
+ }
+ am_vector = new dcomplex[ndit * ndit]();
+ am = new dcomplex*[ndit];
+ for (np_int ai = 0; ai < ndit; ai++) {
+ am[ai] = (am_vector + ai * ndit);
+ MPI_Bcast(am[ai], ndit, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ }
+ MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ lastxi = ((mpidata->rank+1) * xiblock)+1;
+ firstxi = lastxi-xiblock+1;
+ if (lastxi > number_of_scales) lastxi = number_of_scales;
+
+#ifdef USE_MAGMA
+ proc_device = mpidata->rank % device_count;
+#else
+ proc_device = 0;
+#endif
+ MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+}
+
+void ClusterIterationData::mpibcast(const mixMPI *mpidata) {
+ c1->mpibcast(mpidata);
+ const int ndi = c1->nsph * c1->nlim;
+ const np_int ndit = 2 * ndi;
+ MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ for (int ti = 0; ti < 2; ti++) {
+ MPI_Bcast(tqse[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqspe[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqss[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqsps[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ }
+ MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ for (int gi = 0; gi < 3; gi++) {
+ MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ }
+ MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ for (int ci = 0; ci < 4; ci++) {
+ MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ }
+ for (int zi = 0; zi < c1->lm; zi++) {
+ for (int zj = 0; zj < 3; zj++) {
+ for (int zk = 0; zk < 2; zk++) {
+ MPI_Bcast(zpv[zi][zj][zk], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ }
+ }
+ }
+ // since MPI expects an int argument for the number of elements to transfer in one go, transfer am one row at a time
+ for (int ai = 0; ai < ndit; ai++) {
+ MPI_Bcast(am[ai], ndit, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ }
+ MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+}
+#endif
+
+ClusterIterationData::~ClusterIterationData() {
+ const int nsph = c1->nsph;
+ delete[] am_vector;
+ delete[] am;
+ for (int zi = c1->lm - 1; zi > -1; zi--) {
+ for (int zj = 2; zj > -1; zj--) {
+ delete[] zpv[zi][zj][1];
+ delete[] zpv[zi][zj][0];
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
+ delete c1;
+ delete[] gaps;
+ for (int ti = 1; ti > -1; ti--) {
+ delete[] tqse[ti];
+ delete[] tqss[ti];
+ delete[] tqspe[ti];
+ delete[] tqsps[ti];
+ delete[] tqce[ti];
+ delete[] tqcpe[ti];
+ delete[] tqcs[ti];
+ delete[] tqcps[ti];
+ }
+ delete[] tqse;
+ delete[] tqss;
+ delete[] tqspe;
+ delete[] tqsps;
+ delete[] tqce;
+ delete[] tqcpe;
+ delete[] tqcs;
+ delete[] tqcps;
+ delete[] tqev;
+ delete[] tqsv;
+ for (int gi = 2; gi > -1; gi--) {
+ delete[] gapp[gi];
+ delete[] gappm[gi];
+ delete[] gap[gi];
+ delete[] gapm[gi];
+ }
+ delete[] gapp;
+ delete[] gappm;
+ delete[] gap;
+ delete[] gapm;
+ delete[] gapv;
+ delete[] u;
+ delete[] us;
+ delete[] un;
+ delete[] uns;
+ delete[] up;
+ delete[] ups;
+ delete[] unmp;
+ delete[] unsmp;
+ delete[] upmp;
+ delete[] upsmp;
+ delete[] argi;
+ delete[] args;
+ delete[] duk;
+ for (int ci = 3; ci > -1; ci--) {
+ delete[] cextlr[ci];
+ delete[] cext[ci];
+ delete[] cmullr[ci];
+ delete[] cmul[ci];
+ }
+ delete[] cextlr;
+ delete[] cext;
+ delete[] cmullr;
+ delete[] cmul;
+}
+// >>> END OF ClusterIterationData CLASS IMPLEMENTATION <<<
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 2012d7f86ce35a1d3c770e38eb4ac16087b9594a..a2e8fee9642773c404301ce3056a20b9881dab47 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -65,162 +65,6 @@ public:
~mixMPI();
};
-/*! \brief A data structure representing the information used for a single scale
- * of the CLUSTER case.
- */
-class ClusterIterationData {
-public:
- //! \brief Pointer to a ParticleDescriptor structure.
- ParticleDescriptor *c1;
- //! \brief Vector of geometric asymmetry factors.
- double *gaps;
- //! \brief Components of extinction contribution to radiation torque on a single sphere along k.
- double **tqse;
- //! \brief Components of polarized extinction contribution to radiation torque on a single sphere along k.
- dcomplex **tqspe;
- //! \brief Components of scattering contribution to radiation torque on a single sphere along k.
- double **tqss;
- //! \brief Components of polarized scattering contribution to radiation torque on a single sphere along k.
- dcomplex **tqsps;
- //! \brief L-dependent coefficients of the geometric asymmetry parameter.
- double ****zpv;
- //! \brief Mean geometric asymmetry parameters.
- double **gapm;
- //! \brief Mean geometric asymmetry parameters referred to polarization plane.
- dcomplex **gappm;
- //! \brief Imaginary part of the harmonic functions argument.
- double *argi;
- //! \brief Argument of the harmonic functions referred to the scattering plane.
- double *args;
- //! \brief Geometric asymmetry parameters.
- double **gap;
- //! \brief Geometric asymmetry parameters referred to polarization plane.
- dcomplex **gapp;
- //! \brief Components of extinction contribution to radiation torque on the cluster along k.
- double **tqce;
- //! \brief Components of extinction contribution to radiation torque on the cluster along k referred to polarization plane.
- dcomplex **tqcpe;
- //! \brief Components of scattering contribution to radiation torque on the cluster along k.
- double **tqcs;
- //! \brief Components of scattering contribution to radiation torque on the cluster along k referred to polarization plane.
- dcomplex **tqcps;
- //! \brief Variation of unitary radiation vector. QUESTION: correct?
- double *duk;
- //! \brief Cluster extinction cross-section components referred to scattering plane.
- double **cextlr;
- //! \brief Cluster extinction cross-section components referred to meridional plane.
- double **cext;
- //! \brief Cluster Mueller Transformation Matrix components referred to scattering plane.
- double **cmullr;
- //! \brief Cluster Mueller Transformation Matrix components referred to meridional plane.
- double **cmul;
- //! \brief Geometric asymmetry parameter components.
- double *gapv;
- //! \brief Radiation extinction torque components.
- double *tqev;
- //! \brief Radiation scattering torque components.
- double *tqsv;
- //! \brief Incident unitary vector components.
- double *u;
- //! \brief Scattered unitary vector components.
- double *us;
- //! \brief Normal unitary vector components.
- double *un;
- //! \brief Normal scattered unitary vector components.
- double *uns;
- //! \brief Incident unitary vector components on polarization plane.
- double *up;
- //! \brief Scattered unitary vector components on polarization plane.
- double *ups;
- //! \brief Mean unitary vector components normal to polarization plane.
- double *unmp;
- //! \brief Mean scattered unitary vector components normal to polarization plane.
- double *unsmp;
- //! \brief Mean incident unitary vector components on polarization plane.
- double *upmp;
- //! \brief Mean scattered unitary vector components on polarization plane.
- double *upsmp;
- //! \brief Scattering angle.
- double scan;
- //! \brief Control parameter on incidence direction referred to meridional plane.
- double cfmp;
- //! \brief Control parameter on scattering direction referred to meridional plane.
- double sfmp;
- //! \brief Control parameter on incidence direction referred to scattering plane.
- double cfsp;
- //! \brief Control parameter on scattering direction referred to scattering plane.
- double sfsp;
- //! \brief SQSFI = XI^-2
- double sqsfi;
- //! \brief Vectorized scattering coefficient matrix.
- dcomplex *am_vector;
- //! \brief Scattering coefficient matrix.
- dcomplex **am;
- //! \brief Argument of harmonic functions. QUESTION: correct?
- dcomplex arg;
- //! \brief Vacuum magnitude of wave vector.
- double vk;
- //! \brief Wave number.
- double wn;
- //! \brief Normalization scale. QUESTION: correct?
- double xip;
- //! \brief Number of scales (wavelengths) to be computed.
- int number_of_scales;
- //! \brief Size of the block of scales handled by the current process.
- int xiblock;
- //! \brief Index of the first scale handled by the current process.
- int firstxi;
- //! \brief Index of the last scale handled by the current process.
- int lastxi;
- //! \brief ID of the GPU used by one MPI process.
- int proc_device;
- //! \brief Refinement mode selction flag.
- int refinemode;
- //! \brief Maximum number of refinement iterations.
- int maxrefiters;
- //! \brief Required accuracy level.
- double accuracygoal;
-
- /*! \brief `ClusterIterationData` default instance constructor.
- *
- * \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
- * \param sconf: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` object.
- * \param mpidata: `mixMPI *` Pointer to a `mixMPI` object.
- * \param device_count: `const int` Number of offload devices available on the system.
- */
- ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count);
-
- /*! \brief `ClusterIterationData` copy constructor.
- *
- * \param rhs: `const ClusterIterationData &` Reference to the `ClusterIterationData` object to be copied.
- */
- ClusterIterationData(const ClusterIterationData& rhs);
-
-#ifdef MPI_VERSION
- /*! \brief `ClusterIterationData` MPI constructor.
- *
- * \param mpidata: `const mixMPI *` Pointer to a `mixMPI` instance.
- * \param device_count: `const int` Number of offload devices available on the system.
- */
- ClusterIterationData(const mixMPI *mpidata, const int device_count);
-
- /*! \brief Broadcast over MPI the ClusterIterationData instance from MPI process 0 to all others.
- *
- * When using MPI, the initial ClusterIterationData instance created by MPI process 0
- * needs to be replicated on all other processes. This function sends it using
- * MPI broadcast calls. The MPI broadcast calls in this function must match those
- * in the constructor using the mixMPI pointer.
- *
- * \param mpidata: `mixMPI *` Pointer to the mpi structure used to do the MPI broadcast.
- */
- void mpibcast(const mixMPI *mpidata);
-#endif // MPI_VERSION
-
- /*! \brief `ClusterIterationData` instance destroyer.
- */
- ~ClusterIterationData();
-};
-
/*! \brief Basic data structure describing the particle model and its interaction with fields.
*
* This class forms a base of the data structure collections that are used by the
diff --git a/src/include/IterationData.h b/src/include/IterationData.h
new file mode 100644
index 0000000000000000000000000000000000000000..46ebcee36ab9666511cd005795e035c32f22036d
--- /dev/null
+++ b/src/include/IterationData.h
@@ -0,0 +1,184 @@
+/* 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 IterationData.h
+ *
+ * \brief Multi-process communication data structures.
+ *
+ */
+
+#ifndef INCLUDE_ITERATION_DATA_H_
+#define INCLUDE_ITERATION_DATA_H_
+
+// >>> DEFINITION OF ClusterIterationData CLASS <<<
+/*! \brief A data structure representing the information used for a single scale
+ * of the CLUSTER case.
+ */
+class ClusterIterationData {
+public:
+ //! \brief Pointer to a ParticleDescriptor structure.
+ ParticleDescriptor *c1;
+ //! \brief Vector of geometric asymmetry factors.
+ double *gaps;
+ //! \brief Components of extinction contribution to radiation torque on a single sphere along k.
+ double **tqse;
+ //! \brief Components of polarized extinction contribution to radiation torque on a single sphere along k.
+ dcomplex **tqspe;
+ //! \brief Components of scattering contribution to radiation torque on a single sphere along k.
+ double **tqss;
+ //! \brief Components of polarized scattering contribution to radiation torque on a single sphere along k.
+ dcomplex **tqsps;
+ //! \brief L-dependent coefficients of the geometric asymmetry parameter.
+ double ****zpv;
+ //! \brief Mean geometric asymmetry parameters.
+ double **gapm;
+ //! \brief Mean geometric asymmetry parameters referred to polarization plane.
+ dcomplex **gappm;
+ //! \brief Imaginary part of the harmonic functions argument.
+ double *argi;
+ //! \brief Argument of the harmonic functions referred to the scattering plane.
+ double *args;
+ //! \brief Geometric asymmetry parameters.
+ double **gap;
+ //! \brief Geometric asymmetry parameters referred to polarization plane.
+ dcomplex **gapp;
+ //! \brief Components of extinction contribution to radiation torque on the cluster along k.
+ double **tqce;
+ //! \brief Components of extinction contribution to radiation torque on the cluster along k referred to polarization plane.
+ dcomplex **tqcpe;
+ //! \brief Components of scattering contribution to radiation torque on the cluster along k.
+ double **tqcs;
+ //! \brief Components of scattering contribution to radiation torque on the cluster along k referred to polarization plane.
+ dcomplex **tqcps;
+ //! \brief Variation of unitary radiation vector. QUESTION: correct?
+ double *duk;
+ //! \brief Cluster extinction cross-section components referred to scattering plane.
+ double **cextlr;
+ //! \brief Cluster extinction cross-section components referred to meridional plane.
+ double **cext;
+ //! \brief Cluster Mueller Transformation Matrix components referred to scattering plane.
+ double **cmullr;
+ //! \brief Cluster Mueller Transformation Matrix components referred to meridional plane.
+ double **cmul;
+ //! \brief Geometric asymmetry parameter components.
+ double *gapv;
+ //! \brief Radiation extinction torque components.
+ double *tqev;
+ //! \brief Radiation scattering torque components.
+ double *tqsv;
+ //! \brief Incident unitary vector components.
+ double *u;
+ //! \brief Scattered unitary vector components.
+ double *us;
+ //! \brief Normal unitary vector components.
+ double *un;
+ //! \brief Normal scattered unitary vector components.
+ double *uns;
+ //! \brief Incident unitary vector components on polarization plane.
+ double *up;
+ //! \brief Scattered unitary vector components on polarization plane.
+ double *ups;
+ //! \brief Mean unitary vector components normal to polarization plane.
+ double *unmp;
+ //! \brief Mean scattered unitary vector components normal to polarization plane.
+ double *unsmp;
+ //! \brief Mean incident unitary vector components on polarization plane.
+ double *upmp;
+ //! \brief Mean scattered unitary vector components on polarization plane.
+ double *upsmp;
+ //! \brief Scattering angle.
+ double scan;
+ //! \brief Control parameter on incidence direction referred to meridional plane.
+ double cfmp;
+ //! \brief Control parameter on scattering direction referred to meridional plane.
+ double sfmp;
+ //! \brief Control parameter on incidence direction referred to scattering plane.
+ double cfsp;
+ //! \brief Control parameter on scattering direction referred to scattering plane.
+ double sfsp;
+ //! \brief SQSFI = XI^-2
+ double sqsfi;
+ //! \brief Vectorized scattering coefficient matrix.
+ dcomplex *am_vector;
+ //! \brief Scattering coefficient matrix.
+ dcomplex **am;
+ //! \brief Argument of harmonic functions. QUESTION: correct?
+ dcomplex arg;
+ //! \brief Vacuum magnitude of wave vector.
+ double vk;
+ //! \brief Wave number.
+ double wn;
+ //! \brief Normalization scale. QUESTION: correct?
+ double xip;
+ //! \brief Number of scales (wavelengths) to be computed.
+ int number_of_scales;
+ //! \brief Size of the block of scales handled by the current process.
+ int xiblock;
+ //! \brief Index of the first scale handled by the current process.
+ int firstxi;
+ //! \brief Index of the last scale handled by the current process.
+ int lastxi;
+ //! \brief ID of the GPU used by one MPI process.
+ int proc_device;
+ //! \brief Refinement mode selction flag.
+ int refinemode;
+ //! \brief Maximum number of refinement iterations.
+ int maxrefiters;
+ //! \brief Required accuracy level.
+ double accuracygoal;
+
+ /*! \brief `ClusterIterationData` default instance constructor.
+ *
+ * \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
+ * \param sconf: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` object.
+ * \param mpidata: `mixMPI *` Pointer to a `mixMPI` object.
+ * \param device_count: `const int` Number of offload devices available on the system.
+ */
+ ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count);
+
+ /*! \brief `ClusterIterationData` copy constructor.
+ *
+ * \param rhs: `const ClusterIterationData &` Reference to the `ClusterIterationData` object to be copied.
+ */
+ ClusterIterationData(const ClusterIterationData& rhs);
+
+#ifdef MPI_VERSION
+ /*! \brief `ClusterIterationData` MPI constructor.
+ *
+ * \param mpidata: `const mixMPI *` Pointer to a `mixMPI` instance.
+ * \param device_count: `const int` Number of offload devices available on the system.
+ */
+ ClusterIterationData(const mixMPI *mpidata, const int device_count);
+
+ /*! \brief Broadcast over MPI the ClusterIterationData instance from MPI process 0 to all others.
+ *
+ * When using MPI, the initial ClusterIterationData instance created by MPI process 0
+ * needs to be replicated on all other processes. This function sends it using
+ * MPI broadcast calls. The MPI broadcast calls in this function must match those
+ * in the constructor using the mixMPI pointer.
+ *
+ * \param mpidata: `mixMPI *` Pointer to the mpi structure used to do the MPI broadcast.
+ */
+ void mpibcast(const mixMPI *mpidata);
+#endif // MPI_VERSION
+
+ /*! \brief `ClusterIterationData` instance destroyer.
+ */
+ ~ClusterIterationData();
+};
+// >>> END OF DEFINITION OF ClusterIterationData CLASS <<<
+
+#endif // INCLUDE_ITERATION_DATA_H_
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index e3b7eea555653d1186c28044ec36f177c90a87f7..cfdec141721873bb3ecddaaa80614fefa6a87959 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -61,537 +61,6 @@ mixMPI::mixMPI(const mixMPI& rhs) {
mixMPI::~mixMPI() {
}
-ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count) {
- c1 = new ParticleDescriptorCluster(gconf, sconf);
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
- gaps = new double[c1->nsph]();
- tqev = new double[3]();
- tqsv = new double[3]();
- tqse = new double*[2];
- tqspe = new dcomplex*[2];
- tqss = new double*[2];
- tqsps = new dcomplex*[2];
- tqce = new double*[2];
- tqcpe = new dcomplex*[2];
- tqcs = new double*[2];
- tqcps = new dcomplex*[2];
- for (int ti = 0; ti < 2; ti++) {
- tqse[ti] = new double[c1->nsph]();
- tqspe[ti] = new dcomplex[c1->nsph]();
- tqss[ti] = new double[c1->nsph]();
- tqsps[ti] = new dcomplex[c1->nsph]();
- tqce[ti] = new double[3]();
- tqcpe[ti] = new dcomplex[3]();
- tqcs[ti] = new double[3]();
- tqcps[ti] = new dcomplex[3]();
- }
- gapv = new double[3]();
- gapp = new dcomplex*[3];
- gappm = new dcomplex*[3];
- gap = new double*[3];
- gapm = new double*[3];
- for (int gi = 0; gi < 3; gi++) {
- gapp[gi] = new dcomplex[2]();
- gappm[gi] = new dcomplex[2]();
- gap[gi] = new double[2]();
- gapm[gi] = new double[2]();
- }
- u = new double[3]();
- us = new double[3]();
- un = new double[3]();
- uns = new double[3]();
- up = new double[3]();
- ups = new double[3]();
- unmp = new double[3]();
- unsmp = new double[3]();
- upmp = new double[3]();
- upsmp = new double[3]();
- argi = new double[1]();
- args = new double[1]();
- duk = new double[3]();
- cextlr = new double*[4];
- cext = new double*[4];
- cmullr = new double*[4];;
- cmul = new double*[4];
- for (int ci = 0; ci < 4; ci++) {
- cextlr[ci] = new double[4]();
- cext[ci] = new double[4]();
- cmullr[ci] = new double[4]();
- cmul[ci] = new double[4]();
- }
- zpv = new double***[c1->lm];
- for (int zi = 0; zi < c1->lm; zi++) {
- zpv[zi] = new double**[3];
- for (int zj = 0; zj < 3; zj++) {
- zpv[zi][zj] = new double*[2];
- for (int zk = 0; zk < 2; zk++) {
- zpv[zi][zj][zk] = new double[2]();
- }
- }
- }
- am_vector = new dcomplex[ndit * ndit]();
- am = new dcomplex*[ndit];
- for (int ai = 0; ai < ndit; ai++) {
- am[ai] = (am_vector + ai * ndit);
- }
-
- arg = 0.0 + 0.0 * I;
- // These are suspect initializations
- scan = 0.0;
- cfmp = 0.0;
- sfmp = 0.0;
- cfsp = 0.0;
- sfsp = 0.0;
- // End of suspect initializations
- wn = sconf->wp / 3.0e8;
- xip = sconf->xip;
- sqsfi = 1.0;
- vk = 0.0;
- number_of_scales = sconf->number_of_scales;
- xiblock = (int) ceil(((double) (sconf->number_of_scales-1))/((double) mpidata->nprocs));
- lastxi = ((mpidata->rank+1) * xiblock)+1;
- firstxi = lastxi-xiblock+1;
- if (lastxi > sconf->number_of_scales) lastxi = sconf->number_of_scales;
-
-#ifdef USE_MAGMA
- proc_device = mpidata->rank % device_count;
-#else
- proc_device = 0;
-#endif
-
- // In the first iteration, if refinement is enabled, determine the number of refinement iterations required to arrive at the target accuracy (if achievable in a reasonable number of iterations)
- refinemode = 2;
- // maxrefiters and accuracygoal should be configurable and preferably set somewhere else
- maxrefiters = 20;
- accuracygoal = 1e-6;
-}
-
-ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
- c1 = new ParticleDescriptorCluster(reinterpret_cast<ParticleDescriptorCluster &>(*(rhs.c1)));
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
- gaps = new double[c1->nsph]();
- for (int gi = 0; gi < c1->nsph; gi++) gaps[gi] = rhs.gaps[gi];
- tqev = new double[3]();
- tqsv = new double[3]();
- for (int ti = 0; ti < 3; ti++) {
- tqev[ti] = rhs.tqev[ti];
- tqsv[ti] = rhs.tqsv[ti];
- }
- tqse = new double*[2];
- tqspe = new dcomplex*[2];
- tqss = new double*[2];
- tqsps = new dcomplex*[2];
- tqce = new double*[2];
- tqcpe = new dcomplex*[2];
- tqcs = new double*[2];
- tqcps = new dcomplex*[2];
- for (int ti = 0; ti < 2; ti++) {
- tqse[ti] = new double[c1->nsph]();
- tqspe[ti] = new dcomplex[c1->nsph]();
- tqss[ti] = new double[c1->nsph]();
- tqsps[ti] = new dcomplex[c1->nsph]();
- for (int tj = 0; tj < c1->nsph; tj++) {
- tqse[ti][tj] = rhs.tqse[ti][tj];
- tqspe[ti][tj] = rhs.tqspe[ti][tj];
- tqss[ti][tj] = rhs.tqss[ti][tj];
- tqsps[ti][tj] = rhs.tqsps[ti][tj];
- }
- tqce[ti] = new double[3]();
- tqcpe[ti] = new dcomplex[3]();
- tqcs[ti] = new double[3]();
- tqcps[ti] = new dcomplex[3]();
- for (int tj = 0; tj < 3; tj++) {
- tqce[ti][tj] = rhs.tqce[ti][tj];
- tqcpe[ti][tj] = rhs.tqcpe[ti][tj];
- tqcs[ti][tj] = rhs.tqcs[ti][tj];
- tqcps[ti][tj] = rhs.tqcps[ti][tj];
- }
- }
- gapv = new double[3]();
- gapp = new dcomplex*[3];
- gappm = new dcomplex*[3];
- gap = new double*[3];
- gapm = new double*[3];
- for (int gi = 0; gi < 3; gi++) {
- gapv[gi] = rhs.gapv[gi];
- gapp[gi] = new dcomplex[2]();
- gappm[gi] = new dcomplex[2]();
- gap[gi] = new double[2]();
- gapm[gi] = new double[2]();
- for (int gj = 0; gj < 2; gj++) {
- gapp[gi][gj] = rhs.gapp[gi][gj];
- gappm[gi][gj] = rhs.gappm[gi][gj];
- gap[gi][gj] = rhs.gap[gi][gj];
- gapm[gi][gj] = rhs.gapm[gi][gj];
- }
- }
- u = new double[3]();
- us = new double[3]();
- un = new double[3]();
- uns = new double[3]();
- up = new double[3]();
- ups = new double[3]();
- unmp = new double[3]();
- unsmp = new double[3]();
- upmp = new double[3]();
- upsmp = new double[3]();
- duk = new double[3]();
- for (int ui = 0; ui < 3; ui++) {
- u[ui] = rhs.u[ui];
- us[ui] = rhs.us[ui];
- un[ui] = rhs.un[ui];
- uns[ui] = rhs.uns[ui];
- up[ui] = rhs.up[ui];
- ups[ui] = rhs.ups[ui];
- unmp[ui] = rhs.unmp[ui];
- unsmp[ui] = rhs.unsmp[ui];
- upmp[ui] = rhs.upmp[ui];
- upsmp[ui] = rhs.upsmp[ui];
- duk[ui] = rhs.duk[ui];
- }
- argi = new double[1]();
- args = new double[1]();
- argi[0] = rhs.argi[0];
- args[0] = rhs.args[0];
- cextlr = new double*[4];
- cext = new double*[4];
- cmullr = new double*[4];;
- cmul = new double*[4];
- for (int ci = 0; ci < 4; ci++) {
- cextlr[ci] = new double[4]();
- cext[ci] = new double[4]();
- cmullr[ci] = new double[4]();
- cmul[ci] = new double[4]();
- for (int cj = 0; cj < 4; cj++) {
- cextlr[ci][cj] = rhs.cextlr[ci][cj];
- cext[ci][cj] = rhs.cext[ci][cj];
- cmullr[ci][cj] = rhs.cmullr[ci][cj];
- cmul[ci][cj] = rhs.cmul[ci][cj];
- }
- }
- zpv = new double***[c1->lm];
- for (int zi = 0; zi < c1->lm; zi++) {
- zpv[zi] = new double**[3];
- for (int zj = 0; zj < 3; zj++) {
- zpv[zi][zj] = new double*[2];
- for (int zk = 0; zk < 2; zk++) {
- zpv[zi][zj][zk] = new double[2]();
- zpv[zi][zj][zk][0] = rhs.zpv[zi][zj][zk][0];
- zpv[zi][zj][zk][1] = rhs.zpv[zi][zj][zk][1];
- }
- }
- }
- am_vector = new dcomplex[ndit * ndit]();
- for (np_int ai = 0; ai < ndit * ndit; ai++) am_vector[ai] = rhs.am_vector[ai];
- am = new dcomplex*[ndit];
- for (np_int ai = 0; ai < ndit; ai++) {
- am[ai] = (am_vector + ai * ndit);
- }
-
- arg = rhs.arg;
- // These are suspect initializations
- scan = rhs.scan;
- cfmp = rhs.cfmp;
- sfmp = rhs.sfmp;
- cfsp = rhs.cfsp;
- sfsp = rhs.sfsp;
- // End of suspect initializations
- wn = rhs.wn;
- xip = rhs.xip;
- sqsfi = rhs.sqsfi;
- vk = rhs.vk;
- firstxi = rhs.firstxi;
- lastxi = rhs.lastxi;
- xiblock = rhs.xiblock;
- number_of_scales = rhs.number_of_scales;
-
- proc_device = rhs.proc_device;
- refinemode = rhs.refinemode;
- maxrefiters = rhs.maxrefiters;
- accuracygoal = rhs.accuracygoal;
-}
-
-#ifdef MPI_VERSION
-ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int device_count) {
- c1 = new ParticleDescriptorCluster(mpidata);
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
- gaps = new double[c1->nsph]();
- MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- tqev = new double[3]();
- tqsv = new double[3]();
- MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- tqse = new double*[2];
- tqspe = new dcomplex*[2];
- tqss = new double*[2];
- tqsps = new dcomplex*[2];
- tqce = new double*[2];
- tqcpe = new dcomplex*[2];
- tqcs = new double*[2];
- tqcps = new dcomplex*[2];
- for (int ti = 0; ti < 2; ti++) {
- tqse[ti] = new double[c1->nsph]();
- tqspe[ti] = new dcomplex[c1->nsph]();
- tqss[ti] = new double[c1->nsph]();
- tqsps[ti] = new dcomplex[c1->nsph]();
- MPI_Bcast(tqse[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqspe[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqss[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqsps[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- tqce[ti] = new double[3]();
- tqcpe[ti] = new dcomplex[3]();
- tqcs[ti] = new double[3]();
- tqcps[ti] = new dcomplex[3]();
- MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- }
- gapv = new double[3]();
- gapp = new dcomplex*[3];
- gappm = new dcomplex*[3];
- gap = new double*[3];
- gapm = new double*[3];
- MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- for (int gi = 0; gi < 3; gi++) {
- gapp[gi] = new dcomplex[2]();
- gappm[gi] = new dcomplex[2]();
- gap[gi] = new double[2]();
- gapm[gi] = new double[2]();
- MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- }
- u = new double[3]();
- us = new double[3]();
- un = new double[3]();
- uns = new double[3]();
- up = new double[3]();
- ups = new double[3]();
- unmp = new double[3]();
- unsmp = new double[3]();
- upmp = new double[3]();
- upsmp = new double[3]();
- duk = new double[3]();
- MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- argi = new double[1]();
- args = new double[1]();
- MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- cextlr = new double*[4];
- cext = new double*[4];
- cmullr = new double*[4];;
- cmul = new double*[4];
- for (int ci = 0; ci < 4; ci++) {
- cextlr[ci] = new double[4]();
- cext[ci] = new double[4]();
- cmullr[ci] = new double[4]();
- cmul[ci] = new double[4]();
- MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- }
- zpv = new double***[c1->lm];
- for (int zi = 0; zi < c1->lm; zi++) {
- zpv[zi] = new double**[3];
- for (int zj = 0; zj < 3; zj++) {
- zpv[zi][zj] = new double*[2];
- for (int zk = 0; zk < 2; zk++) {
- zpv[zi][zj][zk] = new double[2]();
- MPI_Bcast(zpv[zi][zj][zk], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- }
- }
- }
- am_vector = new dcomplex[ndit * ndit]();
- am = new dcomplex*[ndit];
- for (np_int ai = 0; ai < ndit; ai++) {
- am[ai] = (am_vector + ai * ndit);
- MPI_Bcast(am[ai], ndit, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- }
- MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
- lastxi = ((mpidata->rank+1) * xiblock)+1;
- firstxi = lastxi-xiblock+1;
- if (lastxi > number_of_scales) lastxi = number_of_scales;
-
-#ifdef USE_MAGMA
- proc_device = mpidata->rank % device_count;
-#else
- proc_device = 0;
-#endif
- MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
-}
-
-void ClusterIterationData::mpibcast(const mixMPI *mpidata) {
- c1->mpibcast(mpidata);
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
- MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- for (int ti = 0; ti < 2; ti++) {
- MPI_Bcast(tqse[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqspe[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqss[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqsps[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- }
- MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- for (int gi = 0; gi < 3; gi++) {
- MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- }
- MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- for (int ci = 0; ci < 4; ci++) {
- MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- }
- for (int zi = 0; zi < c1->lm; zi++) {
- for (int zj = 0; zj < 3; zj++) {
- for (int zk = 0; zk < 2; zk++) {
- MPI_Bcast(zpv[zi][zj][zk], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- }
- }
- }
- // since MPI expects an int argument for the number of elements to transfer in one go, transfer am one row at a time
- for (int ai = 0; ai < ndit; ai++) {
- MPI_Bcast(am[ai], ndit, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- }
- MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
-}
-#endif
-
-ClusterIterationData::~ClusterIterationData() {
- const int nsph = c1->nsph;
- delete[] am_vector;
- delete[] am;
- for (int zi = c1->lm - 1; zi > -1; zi--) {
- for (int zj = 2; zj > -1; zj--) {
- delete[] zpv[zi][zj][1];
- delete[] zpv[zi][zj][0];
- delete[] zpv[zi][zj];
- }
- delete[] zpv[zi];
- }
- delete[] zpv;
- delete c1;
- delete[] gaps;
- for (int ti = 1; ti > -1; ti--) {
- delete[] tqse[ti];
- delete[] tqss[ti];
- delete[] tqspe[ti];
- delete[] tqsps[ti];
- delete[] tqce[ti];
- delete[] tqcpe[ti];
- delete[] tqcs[ti];
- delete[] tqcps[ti];
- }
- delete[] tqse;
- delete[] tqss;
- delete[] tqspe;
- delete[] tqsps;
- delete[] tqce;
- delete[] tqcpe;
- delete[] tqcs;
- delete[] tqcps;
- delete[] tqev;
- delete[] tqsv;
- for (int gi = 2; gi > -1; gi--) {
- delete[] gapp[gi];
- delete[] gappm[gi];
- delete[] gap[gi];
- delete[] gapm[gi];
- }
- delete[] gapp;
- delete[] gappm;
- delete[] gap;
- delete[] gapm;
- delete[] gapv;
- delete[] u;
- delete[] us;
- delete[] un;
- delete[] uns;
- delete[] up;
- delete[] ups;
- delete[] unmp;
- delete[] unsmp;
- delete[] upmp;
- delete[] upsmp;
- delete[] argi;
- delete[] args;
- delete[] duk;
- for (int ci = 3; ci > -1; ci--) {
- delete[] cextlr[ci];
- delete[] cext[ci];
- delete[] cmullr[ci];
- delete[] cmul[ci];
- }
- delete[] cextlr;
- delete[] cext;
- delete[] cmullr;
- delete[] cmul;
-}
-
// >>> ParticleDescriptor class implementation. <<< //
ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
_class_type = BASE_TYPE;