Commit b8b059d6 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement ZTM subroutine and fix function GHIT

parent 5d960f80

src/cluster/cluster.cpp

+39 −17
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("../../test_data/cluster/DEDFB");

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

	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("../../test_data/cluster/DCLU");

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

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

		// Shortcuts to variables stored in configuration objects

		int nsph = gconf->number_of_spheres;
@@ -91,6 +91,8 @@ void cluster() { 
		C2 *c2 = new C2(nsph, max_ici, npnt, npntts);

		complex<double> **am = new complex<double>*[mxndm];

		for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm];

		const int ndi = c4->nsph * c4->nlim;

		C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);

		// End of global variables for CLU

		fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");

		fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
@@ -225,11 +227,6 @@ void cluster() { 
								c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,

								c1, c2, jer, lcalc, arg

						);

						//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;
@@ -238,19 +235,43 @@ void cluster() { 
					if (jer != 0) break;

				} // i132 loop

				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);

				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());

				if (jer != 0) break; // jxi488 loop: goes to memory clean

				ztm(am, c1, c1ao, c4, c6, c9);

				if (sconf->idfc >= 0) {

					if (jxi488 == gconf->jwtm) {

						int is = 1;

						fstream ttms_file;

						ttms_file.open("c_TTMS", ios::out | ios::binary);

						if (ttms_file.is_open()) {

							ttms_file.write(reinterpret_cast<char *>(&is), sizeof(int));

							ttms_file.write(reinterpret_cast<char *>(&lm), sizeof(int));

							ttms_file.write(reinterpret_cast<char *>(&vk), sizeof(double));

							ttms_file.write(reinterpret_cast<char *>(&exri), sizeof(double));

							int nlemt = 2 * c4->nlem;

							for (int ami = 0; ami < nlemt; ami++) {

								for (int amj = 0; amj < nlemt; amj++) {

									complex<double> value = c1ao->am0m[ami][amj];

									double rval = value.real();

									double ival = value.imag();

									ttms_file.write(reinterpret_cast<char *>(&rval), sizeof(double));

									ttms_file.write(reinterpret_cast<char *>(&ival), sizeof(double));

								}

							}

							ttms_file.close();

						} else { // Could not open TM file. Should never occur.

							printf("ERROR: failed to open TTMS file.\n");

							break;

						}

					}

				}

				// label 156: continue from here



				if (jer != 0) break;

			} // jxi488 loop

			tppoan.close();

		} else { // In case TPPOAN could not be opened. Should never happen.
@@ -263,6 +284,7 @@ void cluster() { 
		delete c3;

		delete c4;

		delete c6;

		delete c9;

		for (int zi = c4->lm - 1; zi > -1; zi--) {

			for (int zj = 2; zj > -1; zj--) {

				delete[] zpv[zi][zj][1];

src/include/clu_subs.h

+103 −48
Original line number Diff line number Diff line
@@ -413,7 +413,7 @@ std::complex<double> ghit( 
						int ny = l3 * l3 + lt54;

						double aors = 1.0 * (l3 + lt54);

						double f3j = (c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;

						cfun = (c1ao->vh[nbhj + lt54 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;

						cfun = (c1ao->vj0[nbhj + lt54 - 1] * c1ao->vyj0[nby + ny - 1]) * f3j;

						csum += cfun;

						jsn *= -1;

						lt54 += 2;
@@ -431,7 +431,7 @@ std::complex<double> ghit( 
							int ny = l3 * l3  + lt60 + m1mm2;

							double aors = 1.0 * (l3 + lt60);

							double f3j = (c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;

							cfun = (c1ao->vh[nbhj + lt60 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;

							cfun = (c1ao->vj0[nbhj + lt60 - 1] * c1ao->vyj0[nby + ny - 1]) * f3j;

							csum += cfun;

						}

						// label 60
@@ -469,33 +469,6 @@ std::complex<double> ghit( 
void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {

	std::complex<double> dm, de, cgh, cgk;

	const std::complex<double> cc0(0.0, 0.0);

	// DEBUG CODE

	/*double srxx = 0.0, sryy = 0.0, srzz = 0.0;

	for (int di = 0; di < c4->nsph; di++) {

		srxx += c1->rxx[di]; sryy += c1->ryy[di]; srzz += c1->rzz[di];

	}

	printf("DEBUG: in CMS srxx = %lE, sryy = %lE, srzz = %lE\n", srxx, sryy, srzz);

	int sind3j = 0;

	for (int di = 0; di < 9; di++) {

		for (int dj = 0; dj < 8; dj++) sind3j += c1ao->ind3j[di][dj];

	}

	printf("DEBUG: in CMS sind3j = %d\n", sind3j);

	std::complex<double> sv3j0(0.0, 0.0);

	for (int di = 0; di < c4->nv3j; di++) sv3j0 += c1ao->v3j0[di];

	printf("DEBUG: in CMS sv3j0 = (%lE,%lE)\n", sv3j0.real(), sv3j0.imag());

	std::complex<double> svh(0.0, 0.0);

	int dsize = (c4->nsph * c4->nsph - 1) * c4->litpo;

	for (int di = 0; di < dsize; di++) svh += c1ao->vh[di];

	printf("DEBUG: in CMS svh = (%lE,%lE)\n", svh.real(), svh.imag());

	std::complex<double> svyhj(0.0, 0.0);

	dsize = (c4->nsph * c4->nsph - 1) * c4->litpos;

	for (int di = 0; di < dsize; di++) svyhj += c1ao->vyhj[di];

	printf("DEBUG: in CMS svyhj = (%lE,%lE)\n", svyhj.real(), svyhj.imag());

	std::complex<double> svyj0(0.0, 0.0);

	dsize = c4->nsph * c4->lmtpos;

	for (int di = 0; di < dsize; di++) svyj0 += c1ao->vyj0[di];

	printf("DEBUG: in CMS svyj0 = (%lE,%lE)\n", svyj0.real(), svyj0.imag());

	// END DEBUG */

	int ndi = c4->nsph * c4->nlim;

	int nbl = 0;

	int nsphmo = c4->nsph - 1;
@@ -531,18 +504,6 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) { 
							int i2e = in2 + ilm2e;

							int j2 = in1 + ilm2;

							int j2e = in1 + ilm2e;

							bool make_break = false;

							if (i1 == 5 && i2 == 487) make_break = true;

							if (i1 == 5 && i2e == 487) make_break = true;

							if (i1e == 5 && i2 == 487) make_break = true;

							if (i1e == 5 && i2e == 487) make_break = true;

							if (j1 == 5 && j2 == 487) make_break = true;

							if (j1 == 5 && j2e == 487) make_break = true;

							if (j1e == 5 && j2 == 487) make_break = true;

							if (j1e == 5 && j2e == 487) make_break = true;

							if (make_break) {

								printf("DEBUG: this is a diverging place.\n");

							}

							cgh = ghit(0, 0, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);

							cgk = ghit(0, 1, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);

							am[i1 - 1][i2 - 1] = cgh;
@@ -585,9 +546,6 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) { 
			} // im1 loop

		} // l1 loop

	} // n1 loop

	double srac3j = 0.0;

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

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

}



/*! \brief C++ porting of HJV
@@ -1013,10 +971,6 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) { 
			for (int iv38 = 1; iv38 <= c4->litpos; iv38++) {

				c1ao->vyhj[iv38 + ivy - 1] = dconjg(ylm[iv38 - 1]);

			} // iv38 loop

			//std::complex<double> svyhj(0.0, 0.0);

			//int dsize = (c4->nsph * c4->nsph - 1) * c4->litpos;

			//for (int di = 0; di < dsize; di++) svyhj += c1ao->vyhj[di];

			//printf("DEBUG: in STR svyhj = (%lE,%lE)\n", svyhj.real(), svyhj.imag());

			ivy += c4->litpos;

		} // ns40 loop

	} // nf40 loop
@@ -1040,6 +994,107 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) { 
	delete[] ylm;

}



/*! \brief C++ porting of ZTM

 *

 * \param am: Matrix of complex.

 * \param c1: `C1 *` Pointer to a C1 instance.

 * \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.

 * \param c4: `C4 *` Pointer to a C4 structure.

 * \param c6: `C6 *` Pointer to a C6 instance.

 * \param c9: `C9 *` Pointer to a C9 instance.

 */

void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {

	std::complex<double> gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;

	const std::complex<double> cc0(0.0, 0.0);

	int ndi = c4->nsph * c4->nlim;

	int i2 = 0;

	for (int n2 = 1; n2 <= c4->nsph; n2++) { // GPU portable?

		for (int l2 = 1; l2 <= c4->li; l2++) {

			int l2tpo = l2 + l2 + 1;

			int m2 = -l2 - 1;

			for (int im2 = 1; im2 <= l2tpo; im2++) {

				m2++;

				i2++;

				int i3 = 0;

				for (int l3 = 1; l3 <= c4->le; l3++) {

					int l3tpo = l3 + l3 + 1;

					int m3 = -l3 - 1;

					for (int im3 = 1; im3 <= l3tpo; im3++) {

						m3++;

						i3++;

						c9->gis[i2 - 1][i3 - 1] = ghit(2, 0, n2, l2, m2, l3, m3, c1, c1ao, c4, c6);

						c9->gls[i2 - 1][i3 - 1] = ghit(2, 1, n2, l2, m2, l3, m3, c1, c1ao, c4, c6);

					} // im3 loop

				} // l3 loop

			} // im2 loop

		} // l2 loop

	} // n2 loop

	/* // DEBUG CODE

	std::complex<double> dbgtst(0.0, 0.0);

	for (int di = 0; di < ndi; di++) {

		for (int dj = 0; dj < c4->nlem; dj++) dbgtst += c9->gis[di][dj];

	}

	printf("DEBUG: in ZTM init sgis = (%lE, %lE)\n", dbgtst.real(), dbgtst.imag());

	dbgtst = std::complex<double>(0.0, 0.0);

	for (int di = 0; di < ndi; di++) {

		for (int dj = 0; dj < c4->nlem; dj++) dbgtst += c9->gls[di][dj];

	}

	printf("DEBUG: in ZTM init sgls = (%lE, %lE)\n", dbgtst.real(), dbgtst.imag());

	// END DEBUG */

	for (int i1 = 1; i1 <= ndi; i1++) { // GPU portable?

		int i1e = i1 + ndi;

		for (int i3 = 1; i3 <= c4->nlem; i3++) {

			int i3e = i3 + c4->nlem;

			sum1 = cc0;

			sum2 = cc0;

			sum3 = cc0;

			sum4 = cc0;

			for (int i2 = 1; i2 <= ndi; i2++) {

				int i2e = i2 + ndi;

				gie = c9->gis[i2 - 1][i3 - 1];

				gle = c9->gls[i2 - 1][i3 - 1];

				a1 = am[i1 - 1][i2 - 1];

				a2 = am[i1 - 1][i2e - 1];

				a3 = am[i1e - 1][i2 - 1];

				a4 = am[i1e - 1][i2e - 1];

				sum1 += (a1 * gie + a2 * gle);

				sum2 += (a1 * gle + a2 * gie);

				sum3 += (a3 * gie + a4 * gle);

				sum4 += (a3 * gle + a4 * gie);

			} // i2 loop

			c9->sam[i1 - 1][i3 - 1] = sum1;

			c9->sam[i1 - 1][i3e - 1] = sum2;

			c9->sam[i1e - 1][i3 - 1] = sum3;

			c9->sam[i1e - 1][i3e - 1] = sum4;

		} // i3 loop

	} // i1 loop

	for (int i1 = 1; i1 <= ndi; i1++) {

		for (int i0 = 1; i0 <= c4->nlem; i0++) {

			c9->gis[i1 - 1][i0 - 1] = dconjg(c9->gis[i1 - 1][i0 - 1]);

			c9->gls[i1 - 1][i0 - 1] = dconjg(c9->gls[i1 - 1][i0 - 1]);

		} // i0 loop

	} // i1 loop

	int nlemt = c4->nlem + c4->nlem;

	for (int i0 = 1; i0 <= c4->nlem; i0++) {

		int i0e = i0 + c4->nlem;

		for (int i3 = 1; i3 <= nlemt; i3++) {

			sum1 = cc0;

			sum2 = cc0;

			for (int i1 = 1; i1 <= ndi; i1 ++) {

				int i1e = i1 + ndi;

				a1 = c9->sam[i1 - 1][i3 - 1];

				a2 = c9->sam[i1e - 1][i3 - 1];

				gie = c9->gis[i1 - 1][i0 - 1];

				gle = c9->gls[i1 - 1][i0 - 1];

				sum1 += (a1 * gie + a2 * gle);

				sum2 += (a1 * gle + a2 * gie);

			} // i1 loop

			c1ao->am0m[i0 - 1][i3 - 1] = -sum1;

			c1ao->am0m[i0e - 1][i3 - 1] = -sum2;

		} // i3 loop

	} // i0 loop

}



/*! \brief Write a matrix to a log file (debug function).

 *

 *  \param mat: Matrix of complex.