From b8b059d6844682b444e3d1c892f03d4439216fe7 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Mon, 27 Nov 2023 16:31:47 +0100
Subject: [PATCH] Implement ZTM subroutine and fix function GHIT
---
src/cluster/cluster.cpp | 56 ++++++++++-----
src/include/clu_subs.h | 151 +++++++++++++++++++++++++++-------------
2 files changed, 142 insertions(+), 65 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 03038fc2..ad983afc 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -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];
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index dd1ee920..c10aecb1 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -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.
--
GitLab