Commit 5d960f80 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Set up execution of cluster case

parent 907c8beb

src/cluster/Makefile

+26 −1
Original line number Diff line number Diff line
@@ -2,8 +2,12 @@ BUILDDIR=../../build/cluster 
FC=gfortran

FCFLAGS=-std=legacy -O3

LFLAGS=

LFLAGS=

CXX=g++

CXXFLAGS=-O0 -ggdb -pg -coverage

CXXLFLAGS=



all: clu edfb

all: clu edfb np_cluster



clu: clu.o

	$(FC) $(FCFLAGS) -o $(BUILDDIR)/clu $(BUILDDIR)/clu.o $(LFLAGS)
@@ -11,6 +15,27 @@ clu: clu.o 
edfb: edfb.o

	$(FC) $(FCFLAGS) -o $(BUILDDIR)/edfb $(BUILDDIR)/edfb.o $(LFLAGS)



np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o

	$(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o



$(BUILDDIR)/np_cluster.o:

	$(CXX) $(CXXFLAGS) -c ../np_cluster.cpp -o $(BUILDDIR)/np_cluster.o



$(BUILDDIR)/Commons.o:

	$(CXX) $(CXXFLAGS) -c ../libnptm/Commons.cpp -o $(BUILDDIR)/Commons.o



$(BUILDDIR)/Configuration.o:

	$(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o



$(BUILDDIR)/Parsers.o:

	$(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o



$(BUILDDIR)/cluster.o:

	$(CXX) $(CXXFLAGS) -c cluster.cpp -o $(BUILDDIR)/cluster.o



$(BUILDDIR)/sphere.o:

	$(CXX) $(CXXFLAGS) -c ../sphere/sphere.cpp -o $(BUILDDIR)/sphere.o



clean:

	rm -f $(BUILDDIR)/*.o

src/cluster/cluster.cpp

+42 −36
Original line number Diff line number Diff line
@@ -19,13 +19,13 @@ using namespace std; 
//! \brief C++ implementation of CLU

void cluster() {

	printf("INFO: making legacy configuration ...\n");

	ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("DEDFB_clu");

	ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/cluster/DEDFB");

	conf->write_formatted("c_OEDFB_clu");

	conf->write_binary("c_TEDF_clu");

	delete conf;

	printf("INFO: reading binary configuration ...\n");

	ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_clu");

	GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("DCLU");

	GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/cluster/DCLU");

	if (sconf->number_of_spheres == gconf->number_of_spheres) {

		// Shortcuts to variables stored in configuration objects

		int nsph = gconf->number_of_spheres;
@@ -61,43 +61,27 @@ void cluster() { 
			c1->ryy[c1i] = gconf->sph_y[c1i];

			c1->rzz[c1i] = gconf->sph_z[c1i];

			c1->ros[c1i] = sconf->radii_of_spheres[c1i];

			int iogi = c1->iog[c1i];

			if (iogi >= c1i + 1) {

				double gcss = pi * c1->ros[c1i] * c1->ros[c1i];

				c1->gcsv[c1i] = gcss;

				int nsh = c1->nshl[c1i];

				for (int j16 = 1; j16 <= nsh; j16++) {

					c1->rc[c1i][j16- 1] = sconf->rcf[c1i][j16] * c1->ros[c1i];

				}

				c3->gcs += c1->gcsv[c1i - 1];

			}

		}

		C4 *c4 = new C4;

		c4->li = gconf->li;

		c4->le = gconf->le;

		c4->lm = lm;

		c4->lmpo = c4->lm + 1;

		c4->litpo = 2 * gconf->li + 1;

		c4->lm = (c4->li > c4-> le) ? c4->li : c4->le;

		c4->nv3j = (c4->lm * (c4->lm  + 1) * (2 * c4->lm + 7)) / 6;

		c4->nsph = nsph;

		// The following is needed to initialize C1_AddOns

		c4->litpo = c4->li + c4->li + 1;

		c4->litpos = c4->litpo * c4->litpo;

		c4->lmtpo = gconf->li + gconf->le + 1;

		c4->lmtpo = c4->li + c4->le + 1;

		c4->lmtpos = c4->lmtpo * c4->lmtpo;

		c4->nlim = c4->li * (c4->li + 2);

		c4->nlem = c4->le * (c4->le + 2);

		c4->nsph = nsph;

		C6 *c6 = new C6(c4->lmtpo);

		c4->lmpo = c4->lm + 1;

		C1_AddOns *c1ao = new C1_AddOns(c4);

		// End of add-ons initialization

		C6 *c6 = new C6(c4->lmtpo);

		FILE *output = fopen("c_OCLU", "w");

		double ****zpv = new double***[c4->lm];

		for (int zi = 0; zi < c4->lm; zi++) {

			zpv[zi] = new double**[3];

			for (int zj = 0; zj < 3; zj++) {

				zpv[zi][zj] = new double*[2];

				zpv[zi][zj][0] = new double[2];

				zpv[zi][zj][1] = new double[2];

			}

		}

		int jer = 0, lcalc = 0;

		complex<double> arg = complex<double>(0.0, 0.0);

		complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);

		int max_ici = 0;

		for (int insh = 0; insh < nsph; insh++) {

			int nsh = sconf->nshl_vec[insh];
@@ -162,7 +146,19 @@ void cluster() { 
		int nks = nths * nphs;

		int nkks = nk * nks;

		double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;

		str(c1, c1ao, c3, c4, c6);

		str(sconf->rcf, c1, c1ao, c3, c4, c6);

		double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]

		for (int zi = 0; zi < c4->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] = 0.0;

					zpv[zi][zj][zk][1] = 0.0;

				}

			}

		}

		thdps(c4->lm, zpv);

		double exri = sqrt(sconf->exdc);

		double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
@@ -200,10 +196,7 @@ void cluster() { 
					sqsfi = 1.0 / (xi * xi);

					fprintf(output, "  XI=%15.7lE\n", xi);

				}

				// Would call HJV(EXRI, VKARG, JER, LCALC, ARG)

				hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);

				//printf("INFO: calculation went up to %d and jer = %d\n", lcalc, jer);

				//printf("INFO: arg = (%lE,%lE)\n", arg.real(), arg.imag());

				if (jer != 0) {

					fprintf(output, "  STOP IN HJV\n");

					break; // jxi488 loop: goes to memory cleaning and  return
@@ -218,7 +211,6 @@ void cluster() { 
					} else {

						int nsh = sconf->nshl_vec[i132 - 1];

						int ici = (nsh + 1) / 2;

						int size_dc0 = (nsh % 2 == 0) ? ici + 1 : ici;

						if (sconf->idfc == 0) {

							for (int ic = 0; ic < ici; ic++)

								c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
@@ -233,8 +225,11 @@ void cluster() { 
								c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,

								c1, c2, jer, lcalc, arg

						);

						//printf("INFO: DME returned jer = %d , lcalc = %d and arg = (%lE, %lE)\n",

						//		jer, lcalc, arg.real(), arg.imag());

						//for (int idl = 1; idl <= 8; idl++) {

						//	printf("DEBUG: RMI( %d , %d ) = (%lE, %lE)\n", idl, i132, c1->rmi[idl - 1][i132 - 1].real(), c1->rmi[idl - 1][i132 - 1].imag());

						//	printf("DEBUG: REI( %d , %d ) = (%lE, %lE)\n", idl, i132, c1->rei[idl - 1][i132 - 1].real(), c1->rei[idl - 1][i132 - 1].imag());

						//}

						//printf("DEBUG: DME - jer = %d, lcalc = %d, arg = (%lE,%lE)\n", jer, lcalc, arg.real(), arg.imag());

						if (jer != 0) {

							fprintf(output, "  STOP IN DME\n");

							break;
@@ -242,7 +237,18 @@ void cluster() { 
					}

					if (jer != 0) break;

				} // i132 loop

				// Would call CMS(AM)

				cms(am, c1, c1ao, c4, c6);

				if (jxi488 == 1) logmat(am, 960, 960);

				ccsam = summat(am, 960, 960);

				printf("DEBUG: after CMS CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag());

				double srac3j = 0.0;

				for (int di = 0; di < c4->lmtpo; di++) srac3j += c6->rac3j[di];

				printf("DEBUG: after CMS SRAC3J = %lE\n", srac3j);

				//int ndit = 2 * nsph * c4->nlim;

				//lucin(am, mxndm, ndit, jer);

				//ccsam = summat(am, 960, 960);

				//printf("DEBUG: after LUCIN CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag());

				//printf("DEBUG: after LUCIN, JER = %d\n", jer);



				if (jer != 0) break;

			} // jxi488 loop

src/include/Commons.h

+6 −4
Original line number Diff line number Diff line
@@ -163,11 +163,11 @@ struct C4 { 
	int litpo;

	//! \brief QUESTION: definition?

	int litpos;

	//! \brief QUESTION: definition?

	//! \brief Maximum field expansion order plus one. QUESTION: correct?

	int lmpo;

	//! \brief QUESTION: definition?

	//! \brief Twice maximum field expansion order plus one. QUESTION: correct?

	int lmtpo;

	//! \brief QUESTION: definition?

	//! \brief Square of `lmtpo`.

	int lmtpos;

	//! \brief QUESTION: definition?

	int li;
@@ -181,6 +181,8 @@ struct C4 { 
	int lm;

	//! \brief Number of spheres.

	int nsph;

	//! \brief QUESTION: definition?

	int nv3j;

};



/*! \brief Vectors and matrices that are specific to cluster C1 blocks.
@@ -192,7 +194,7 @@ protected: 
	int nsph;

	//! \brief QUESTION: definition?

	int nlemt;

	//! \brief QUESTION: definition?

	//! \brief Maximum expansion order plus one. QUESTION: correct?

	int lmpo;



	/*! \brief Allocate the necessary common vectors depending on configuration.

src/include/clu_subs.h

+404 −186

File changed.

Preview size limit exceeded, changes collapsed.

src/include/sph_subs.h

+11 −11
Original line number Diff line number Diff line
@@ -1046,14 +1046,13 @@ label47: 
 * \param zpv: `double ****`

 */

void thdps(int lm, double ****zpv) {

	// WARNING: unclear nested loop in THDPS

	// The optimized interpretation was adopted here.

	for (int l10 = 1; l10 <= lm; l10++) {

		for (int ilmp = 1; ilmp <= 3; ilmp++) {

			zpv[l10 - 1][ilmp - 1][0][0] = 0.0;

			zpv[l10 - 1][ilmp - 1][0][1] = 0.0;

			zpv[l10 - 1][ilmp - 1][1][0] = 0.0;

			zpv[l10 - 1][ilmp - 1][1][1] = 0.0;

	//double szpv = 0.0;

	for (int l10 = 0; l10 < lm; l10++) {

		for (int ilmp = 0; ilmp < 3; ilmp++) {

			zpv[l10][ilmp][0][0] = 0.0;

			zpv[l10][ilmp][0][1] = 0.0;

			zpv[l10][ilmp][1][0] = 0.0;

			zpv[l10][ilmp][1][1] = 0.0;

		}

	}

	for (int l15 = 1; l15 <= lm; l15++) {
@@ -1061,6 +1060,7 @@ void thdps(int lm, double ****zpv) { 
		double zp = -1.0 / sqrt(xd);

		zpv[l15 - 1][1][0][1] = zp;

		zpv[l15 - 1][1][1][0] = zp;

		//szpv += 2.0 * zp;

	}

	if (lm != 1) {

		for (int l20 = 2; l20 <= lm; l20++) {
@@ -1069,6 +1069,7 @@ void thdps(int lm, double ****zpv) { 
			double zp = sqrt(xn / xd);

			zpv[l20 - 1][0][0][0] = zp;

			zpv[l20 - 1][0][1][1] = zp;

			//szpv += 2.0 * zp;

		}

		int lmmo = lm - 1;

		for (int l25 = 1; l25 <= lmmo; l25++) {
@@ -1077,8 +1078,10 @@ void thdps(int lm, double ****zpv) { 
			double zp = -1.0 * sqrt(xn / xd);

			zpv[l25 - 1][2][0][0] = zp;

			zpv[l25 - 1][2][1][1] = zp;

			//szpv += 2.0 * zp;

		}

	}

	//printf("DEBUG: szpv = %lE\n", szpv);

}



/*! \brief C++ porting of UPVMP
@@ -1412,9 +1415,7 @@ void dme( 
	}

	for (int j43 = 1; j43 <= lipt; j43++) {

		fb[j43 - 1] = rfj[j43 - 1];

		//printf("DEBUG: fb[%d] = (%lE,%lE)\n", j43, fb[j43 - 1].real(), fb[j43 - 1].imag());

		fn[j43 - 1] = rfn[j43 - 1];

		//printf("DEBUG: fn[%d] = (%lE,%lE)\n", j43, fn[j43 - 1].real(), fn[j43 - 1].imag());

	}

	if (nsh <= 1) {

		cri = c2->dc0[0] / exdc;
@@ -1429,7 +1430,6 @@ void dme( 
			ccnb = fn[lpo - 1] * dfbi;

			ccnc = fbi[lpo - 1] * dfb;

			ccnd = fb[lpo - 1] * dfbi;

			// CLUSTER terminated in this loop with li = 8.

			c1->rmi[l60 - 1][i - 1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);

			c1->rei[l60 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);

		}