Skip to content
Snippets Groups Projects
Commit dea86173 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Initiate porting of cluster to C++

parent 7ad971ab
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/cluster/cluster.cpp 0 → 100644
+ 279
− 0
View file @ dea86173
#include <cstdio>
#include <fstream>
#include <string>
#include <complex>
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
#ifndef INCLUDE_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif
using namespace std;
/*
* >>> WARNING: works only for IDFC >= 0, as the original code <<<
*
*/
//! \brief C++ implementation of CLU
void cluster() {
printf("INFO: making legacy configuration ...\n");
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("DCLU");
if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects
int nsph = gconf->number_of_spheres;
int mxndm = gconf->mxndm;
int inpol = gconf->in_pol;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
int iavm = gconf->iavm;
int isam = gconf->meridional_type;
int nxi = sconf->number_of_scales;
double th = gconf->in_theta_start;
double thstp = gconf->in_theta_step;
double thlst = gconf->in_theta_end;
double ths = gconf->sc_theta_start;
double thsstp = gconf->sc_theta_step;
double thslst = gconf->sc_theta_end;
double ph = gconf->in_phi_start;
double phstp = gconf->in_phi_step;
double phlst = gconf->in_phi_end;
double phs = gconf->sc_phi_start;
double phsstp = gconf->sc_phi_step;
double phslst = gconf->sc_phi_end;
double wp = sconf->wp;
// Global variables for CLU
double pi = 2.0 * acos(0.0);
int lm = gconf->l_max;
if (gconf->li > lm) lm = gconf->li;
if (gconf->le > lm) lm = gconf->le;
C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec);
C3 *c3 = new C3();
for (int c1i = 0; c1i < nsph; c1i++) {
c1->rxx[c1i] = gconf->sph_x[c1i];
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->litpos = c4->litpo * c4->litpo;
c4->lmtpo = gconf->li + gconf->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);
C1_AddOns *c1ao = new C1_AddOns(c4);
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);
int max_ici = 0;
for (int insh = 0; insh < nsph; insh++) {
int nsh = sconf->nshl_vec[insh];
int ici = (nsh + 1) / 2;
if (ici > max_ici) max_ici = ici;
}
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];
// 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",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
gconf->sph_x[ri], gconf->sph_y[ri], gconf->sph_z[ri]
);
fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
th, thstp, thlst, ths, thsstp, thslst
);
fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
ph, phstp, phlst, phs, phsstp, phslst
);
fprintf(output, " READ(IR,*)JWTM\n");
fprintf(output, " %5d\n", gconf->jwtm);
fprintf(output, " READ(ITIN)NSPHT\n");
fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
fprintf(output, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
fprintf(output, " READ(ITIN)(XIV(I),I=1,NXI)\n");
fprintf(output, " READ(ITIN)NSHL(I),ROS(I)\n");
fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
fprintf(output, " \n");
double small = 1.0e-3;
int nth = 0, nph = 0;
if (thstp != 0.0) nth = int((thlst - th) / thstp + small);
nth++;
if (phstp != 0.0) nph = int((phlst - ph) / phstp + small);
nph++;
int nths = 0, nphs = 0;
double thsca = 0.0;
if (isam > 1) {
nths = 1;
thsca = ths - th;
} else { // ISAM <= 1
if (thsstp == 0.0) nths = 0;
else nths = int ((thslst - ths) / thsstp + small);
nths++;
}
if (isam >= 1) {
nphs = 1;
} else {
if (phsstp == 0.0) nphs = 0;
else nphs = int((phslst - phs) / phsstp + small);
nphs++;
}
int nk = nth * nph;
int nks = nths * nphs;
int nkks = nk * nks;
double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
str(c1, c1ao, c3, c4, c6);
thdps(c4->lm, zpv);
double exri = sqrt(sconf->exdc);
double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream tppoan;
tppoan.open("c_TPPOAN", ios::out | ios::binary);
if (tppoan.is_open()) {
tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nxi), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nth), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
int nlemt = c4->nlem + c4->nlem;
double wn = wp / 3.0e8;
double sqsfi = 1.0;
if (sconf->idfc < 0) {
vk = sconf->xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
int jaw = 1;
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->scale_vec[jxi488 - 1];
double vkarg = 0.0;
if (sconf->idfc >= 0) {
vk = xi * wn;
vkarg = vk;
fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", vk, xi);
} else {
vkarg = xi * vk;
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
}
for (int i132 = 1; i132 <= nsph; i132++) {
int iogi = c1->iog[i132 - 1];
if (iogi != i132) {
for (int l123 = 1; l123 <= gconf->li; l123++) {
c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
} // l123 loop
} 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];
} else {
if (jxi488 == 1) {
for (int ic = 0; ic < ici; ic++)
c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0];
}
}
if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
dme(
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());
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
break;
}
}
if (jer != 0) break;
} // i132 loop
// Would call CMS(AM)
if (jer != 0) break;
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
printf("ERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
delete c1;
delete c1ao;
delete c3;
delete c4;
delete c6;
for (int zi = c4->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;
for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
delete[] am;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
);
}
delete sconf;
delete gconf;
printf("Done.\n");
}
src/include/clu_subs.h 0 → 100644
+ 740
− 0
View file @ dea86173
/*! \file clu_subs.h
*
* \brief C++ porting of CLU functions and subroutines.
*
* Remember that FORTRAN passes arguments by reference, so, every time we use
* a subroutine call, we need to add a referencing layer to the C++ variable.
* All the functions defined below need to be properly documented and ported
* to C++.
*
* Currently, only basic documenting information about functions and parameter
* types are given, to avoid doxygen warning messages.
*/
#ifndef INCLUDE_COMMONS_H_
#include "Commons.h"
#endif
#ifndef INCLUDE_CLU_SUBS_H_
#define INCLUDE_CLU_SUBS_H_
#include <complex>
// >>> DECLARATION OF SPH_SUBS <<<
extern std::complex<double> dconjg(std::complex<double> value);
extern void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
);
extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
extern void rbf(int n, double x, int &nm, double sj[]);
extern void rnf(int n, double x, int &nm, double sy[]);
extern void thdps(int lm, double ****zpv);
// >>> END OF SPH_SUBS DECLARATION <<<
/*! \brief C++ porting of CGEV
*
* \param ipamo: `int`
* \param mu: `int`
* \param l: `int`
* \param m: `int`
*/
double cgev(int ipamo, int mu, int l, int m) {
double result = 0.0;
double xd = 0.0, xn = 0.0;
if (ipamo == 0) {
if (m != 0 || mu != 0) { // label 10
if (mu != 0) {
xd = 2.0 * l * (l + 1);
if (mu <= 0) {
xn = 1.0 * (l + m) * (l - m + 1);
result = sqrt(xn / xd);
} else { // label 15
xn = 1.0 * (l - m) * (l + m + 1);
result = -sqrt(xn / xd);
}
} else { // label 20
xd = 1.0 * (l + 1) * l;
xn = -1.0 * m;
result = xn / sqrt(xd);
}
}
} else { // label 30
xd = 2.0 * l * (l * 2 - 1);
if (mu < 0) { // label 35. CHECK: is clu.f code line 2466 a switch?
xn = 1.0 * (l - 1 + m) * (l + m);
} else if (mu == 0) { // label 40
xn = 2.0 * (l - m) * (l + m);
} else { // mu > 0, label 45
xn = 1.0 * (l - 1 - m) * (l - m);
}
result = sqrt(xn / xd);
}
return result;
}
/*! \brief C++ porting of R3JJR
*
* \param j2: `int`
* \param j3: `int`
* \param m2: `int`
* \param m3: `int`
* \param c6: `C6 *`
*/
void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
}
/*! \brief C++ porting of GHIT
*
* \param ihi: `int`
* \param ipamo: `int`
* \param nbl: `int`
* \param l1: `int`
* \param m1: `int`
* \param l2: `int`
* \param m2: `int`
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
* \param c6: `C6 *`
*/
std::complex<double> ghit(
int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
C1_AddOns *c1ao, C4 *c4, C6 *c6
) {
/* NBL identifies transfer vector going from N2 to N1;
* IHI=0 for Hankel, IHI=1 for Bessel, IHI=2 for Bessel from origin;
* depending on IHI, IPAM=0 gives H or I, IPAM= 1 gives K or L. */
const std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
const std::complex<double> uim = std::complex<double>(0.0, 1.0);
std::complex<double> result = cc0;
if (ihi == 2) {
if (!(c1->rxx[nbl - 1] != 0.0 || c1->ryy[nbl - 1] != 0.0 || c1->rzz[nbl - 1] != 0.0)) {
if (ipamo != 0) return result;
if (l1 == l2 && m1 == m2) result = std::complex<double>(1.0, 0.0);
return result;
}
}
int l1mp = l1 - ipamo;
int l1po = l1 + 1;
int m1mm2 = m1 - m2;
int m1mm2m = (m1mm2 > 0) ? m1mm2 + 1 : -m1mm2 + 1;
int lminpo = (l2 - l1mp > 0) ? l2 - l1mp + 1 : l1mp - l2 + 1;
int lmaxpo = l2 + l1mp + 1;
int i3j0in = c1ao->ind3j[l1mp][l2 - 1];
int ilin = -1;
if (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ilin = 0;
int isn = 1;
if (m1 % 2 != 0) isn *= -1;
if (lminpo % 2 == 0) {
isn *= -1;
if (l2 > l1mp) isn *= -1;
}
int nblmo = nbl - 1;
if (ihi != 2) {
int nbhj = nblmo * c4->litpo;
int nby = nblmo * c4->litpos;
if (ihi != 1) {
for (int jm24 = 1; jm24 <= 3; jm24++) {
std::complex<double> csum = cc0;
int mu = jm24 - 2;
int mupm1 = mu + m1;
int mupm2 = mu + m2;
if (mupm1 < -l1mp || mupm1 > l1mp || mupm2 < -l2 || mupm2 > l2)
continue; //jm24 loop
int jsn = -isn;
if (mu == 0) jsn = isn;
double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
int i3j0 = i3j0in;
if (mupm1 == 0 && mupm2 == 0) {
int lt14 = lminpo;
while (lt14 <= lmaxpo) {
i3j0 += 1;
int l3 = lt14 - 1;
int ny = l3 * l3 + lt14;
double aors = 1.0 * (l3 + lt14);
double f3j = (
c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] *
sqrt(aors)
) * jsn;
std::complex<double> cfun = (
c1ao->vh[nbhj + lt14 - 1] * c1ao->vyhj[nby + ny - 1]
) * f3j;
csum += cfun;
jsn *= - 1;
lt14 += 2;
}
} else { // label 16
// Would call R3JJR(L1MP, L2, -MUPM1, MUPM2)
int il = ilin;
int lt20 = lminpo;
while (lt20 <= lmaxpo) {
i3j0 += 1;
if (m1mm2m <= lt20) {
il = il + 2;
int l3 = lt20 - 1;
int ny = l3 * l3 + lt20 + m1mm2;
double aors = 1.0 * (l3 + lt20);
double f3j = (
c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
) * jsn;
std::complex<double> cfun = (
c1ao->vh[nbhj + lt20 - 1] * c1ao->vyhj[nby + ny - 1]
) * f3j;
csum += cfun;
}
jsn *= -1;
lt20 += 2;
}
}
csum *= cr;
result += csum;
} // jm24 loop
} else { // label 30, IHI = 1
for (int jm44 = 1; jm44 <= 3; jm44 ++) {
std::complex<double> csum = cc0;
int mu = jm44 - 2;
int mupm1 = mu + m1;
int mupm2 = mu + m2;
if (mupm1 >= -l1mp && mupm1 <= l1mp && mupm2 >= -l2 && mupm2 >= l2) {
int jsn = -isn;
if (mu == 0) jsn = isn;
double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
int i3j0 = i3j0in;
if (mupm1 == 0 && mupm2 == 0) {
int lt34 = lminpo;
while (lt34 <= lmaxpo) {
i3j0++;
int l3 = lt34 - 1;
int ny = l3 * l3 + lt34;
double aors = 1.0 * (l3 + lt34);
std::complex<double> f3j = (
c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
) * jsn;
std::complex<double> cfun = (
c1ao->vj[nbhj + lt34 - 1] * c1ao->vyhj[nby + ny - 1]
) * f3j;
csum += cfun;
jsn *= -1;
lt34 += 2;
}
} else { // label 36
// Would call R3JJR
int il = ilin;
int lt40 = lminpo;
while (lt40 <= lmaxpo) {
i3j0++;
if (m1mm2m <= lt40) {
il += 2;
int l3 = lt40 - 1;
int ny = l3 * l3 + lt40 + m1mm2;
double aors = 1.0 * (l3 + lt40);
std::complex<double> f3j = (
c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
) * jsn;
std::complex<double> cfun = (
c1ao->vj[nbhj + lt40 - 1] * c1ao->vyhj[nby + ny - 1]
) * f3j;
csum += cfun;
}
jsn *= -1;
lt40 += 2;
}
}
// label 42
csum *= cr;
}
result += csum;
} // jm44 loop
}
} else { // label 50, IHI = 2
int nbhj = nblmo * c4->lmtpo;
int nby = nblmo * c4->lmtpos;
for (int jm64 =1; jm64 <= 3; jm64++) {
std::complex<double> csum = cc0;
int mu = jm64 - 2;
int mupm1 = mu + m1;
int mupm2 = mu + m2;
if (mupm1 >= -l1mp && mupm1 < l1mp && mupm2 >= -l2 && mupm2 <= l2) {
int jsn = - isn;
if (mu == 0) jsn = isn;
double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
int i3j0 = i3j0in;
if (mupm1 == 0 && mupm2 == 0) {
int lt54 = lminpo;
while (lt54 <= lmaxpo) {
i3j0++;
int l3 = lt54 - 1;
int ny = l3 * l3 + lt54;
double aors = 1.0 * (l3 + lt54);
std::complex<double> f3j = (
c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
) * jsn;
std::complex<double> cfun = (
c1ao->vj0[nbhj + lt54 - 1] * c1ao->vyj0[nby + ny - 1]
) * f3j;
csum += cfun;
jsn *= -1;
lt54 += 2;
}
} else { // label 56
// Would call R3JJR(L1MP, L2, -MUPM1, MUPM2)
int il = ilin;
int lt60 = lminpo;
while (lt60 <= lmaxpo) {
i3j0++;
if (m1mm2m <= lt60) {
il += 2;
int l3 = lt60 - 1;
int ny = l3 * l3 + lt60 + m1mm2;
double aors = 1.0 * (l3 + lt60);
std::complex<double> f3j = (
c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
) * jsn;
std::complex<double> cfun = (
c1ao->vj0[nbhj + lt60 - 1] * c1ao->vyj0[nby + ny - 1]
) * f3j;
csum += cfun;
}
jsn *= -1;
lt60 += 2;
}
}
// label 62
csum *= cr;
}
result += csum;
} // jm64 loop
}
// label 70
const double four_pi = acos(0.0) * 8.0;
if (ipamo != 1) {
double cr = sqrt(four_pi * (l1 + l1po) * (l2 + l2 + 1));
result *= cr;
} else {
double cr = sqrt(four_pi * (l1 + l1mp) * (l1 + l1po) * (l2 + l2 + 1) / l1po);
result *= (cr * uim);
}
return result;
}
/*! \brief C++ porting of CMS
*
* \param am: Matrix of complex.
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
* \param c6: `C6 *`
*/
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 = std::complex<double>(0.0, 0.0);
int ndi = c4->nsph * c4->nlim;
int nbl = 0;
int nsphmo = c4->nsph - 1;
for (int n1 = 1; n1 <= nsphmo; n1++) { // GPU portable?
int in1 = (n1 - 1) * c4->nlim;
int n1po = n1 + 1;
for (int n2 = n1po; n2 <= c4->nsph; n2++) {
int in2 = (n2 - 1) * c4->nlim;
nbl++;
for (int l1 = 1; l1 <= c4->li; l1++) {
int l1po = l1 + 1;
int il1 = l1po * l1;
int l1tpo = l1po + l1;
for (int im1 =1; im1 <= l1tpo; il1++) {
int m1 = im1 - l1po;
int ilm1 = il1 + m1;
int ilm1e = ilm1 + ndi;
int i1 = in1 + ilm1;
int i1e = in1 + ilm1e;
int j1 = in2 + ilm1;
int j1e = in2 + ilm1e;
for (int l2 =1; l2 <= c4->li; l2++) {
int l2po = l2 + 1;
int il2 = l2po * l2;
int l2tpo = l2po + l2;
int ish = ((l2 + l1) % 2 == 0) ? 1 : -1;
int isk = -ish;
for (int im2 = 1; im2 <= l2tpo; im2++) {
int m2 = im2 - l2po;
int ilm2 = il2 + m2;
int ilm2e = ilm2 + ndi;
int i2 = in2 + ilm2;
int i2e = in2 + ilm2e;
int j2 = in1 + ilm2;
int j2e = in1 + ilm2e;
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;
am[i1 - 1][i2e - 1] = cgk;
am[i1e - 1][i2 - 1] = cgk;
am[i1e - 1][i2e - 1] = cgh;
am[j1 - 1][j2 - 1] = cgh * (1.0 * ish);
am[j1 - 1][j2e - 1] = cgk * (1.0 * isk);
am[j1e - 1][j2 - 1] = cgk * (1.0 * isk);
am[j1e - 1][j2e - 1] = cgh * (1.0 * ish);
}
}
} // im1 loop
} // l1 loop
} // n2 loop
} // n1 loop
for (int n1 = 1; n1 <= c4->nsph; n1++) { // GPU portable?
int in1 = (n1 - 1) * c4->nlim;
for (int l1 = 1; l1 <= c4->li; l1++) {
dm = c1->rmi[l1 - 1][n1 - 1];
de = c1->rei[l1 - 1][n1 - 1];
int l1po = l1 + 1;
int il1 = l1po * l1;
int l1tpo = l1po + l1;
for (int im1 = 1; im1 <= l1tpo; im1++) {
int m1 = im1 - l1po;
int ilm1 = il1 + m1;
int i1 = in1 + ilm1;
int i1e = i1 + ndi;
for (int ilm2 = 1; ilm2 <= c4->nlim; ilm2++) {
int i2 = in1 + ilm2;
int i2e = i2 + ndi;
am[i1 - 1][i2 - 1] = cc0;
am[i1 - 1][i2e - 1] = cc0;
am[i1e - 1][i2 - 1] = cc0;
am[i1e - 1][i2e - 1] = cc0;
}
am[i1 - 1][i1 - 1] = dm;
am[i1e - 1][i1e - 1] = de;
} // im1 loop
} // l1 loop
} // n1 loop
}
/*! \brief C++ porting of HJV
*
* \param exri: `double`
* \param vk: `double`
* \param jer: `int &`
* \param lcalc: `int &`
* \param arg: `complex\<double\> &`
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
*/
void hjv(
double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
C1 *c1, C1_AddOns *c1ao, C4 *c4
) {
int nsphmo = c4->nsph - 1;
int lit = c4->li + c4->li;
int lmt = c4->lm + c4->lm;
const int rfj_size = (lit > lmt) ? lit : lmt;
const int rfn_size = c4->litpo;
double *rfj, *rfn;
rfj = new double[rfj_size];
rfn = new double[rfn_size];
jer = 0;
int ivhb = 0;
for (int nf40 = 1; nf40 <= nsphmo; nf40++) { // GPU portable?
int nfpo = nf40 + 1;
for (int ns40 = nfpo; ns40 <= c4->nsph; ns40++) {
double rx = c1->rxx[nf40 - 1] - c1->rxx[ns40 - 1];
double ry = c1->ryy[nf40 - 1] - c1->ryy[ns40 - 1];
double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
double rr = sqrt(rx * rx + ry * ry + rz * rz);
double rarg = rr * vk * exri;
arg = std::complex<double>(rarg, 0.0);
rbf(lit, rarg, lcalc, rfj);
if (lcalc < lit) {
jer = 1;
delete[] rfj;
delete[] rfn;
return;
}
rnf(lit, rarg, lcalc, rfn);
if (lcalc < lit) {
jer = 2;
delete[] rfj;
delete[] rfn;
return;
}
for (int lpo38 = 1; lpo38 <= c4->litpo; lpo38++) {
c1ao->vh[lpo38 + ivhb - 1] = std::complex<double>(rfj[lpo38 - 1], rfn[lpo38 - 1]);
}
ivhb += c4->litpo;
} // ns40 loop
} // nf40 loop
ivhb = 0;
for (int nf50 = 1; nf50 <= c4->nsph; nf50++) {
double rx = c1->rxx[nf50 - 1];
double ry = c1->ryy[nf50 - 1];
double rz = c1->rzz[nf50 - 1];
if (!(rx == 0.0 && ry == 0.0 && rz == 0.0)) {
double rr = sqrt(rx * rx + ry * ry + rz * rz);
double rarg = rr * vk * exri;
rbf(lmt, rarg, lcalc, rfj);
if (lcalc < lmt) {
jer = 3;
delete[] rfj;
delete[] rfn;
return;
}
for (int lpo47 = 1; lpo47 <= c4->lmtpo; lpo47++) {
c1ao->vj0[lpo47 + ivhb - 1] = rfj[lpo47 - 1];
}
}
ivhb += c4->lmtpo;
} // nf50 loop
delete[] rfj;
delete[] rfn;
}
/*! \brief C++ porting of POLAR
*
* \param x: `double`
* \param y: `double`
* \param z: `double`
* \param r: `double &`
* \param cth: `double &`
* \param sth: `double &`
* \param cph: `double &`
* \param sph: `double &`
*/
void polar(
double x, double y, double z, double &r, double &cth, double &sth,
double &cph, double &sph
) {
bool onx = (y == 0.0);
bool ony = (x == 0.0);
bool onz = (onx && ony);
double rho;
if (!(onz || onx || ony)) {
rho = sqrt(x * x + y * y);
cph = x / rho;
sph = y / rho;
} else {
if (onz) {
cph = 1.0;
sph = 0.0;
} else if (onx) {
rho = x;
cph = 1.0;
if (x < 0.0) {
rho *= -1.0;
cph *= -1.0;
}
sph = 0.0;
} else if (ony) {
rho = y;
sph = 1.0;
if (y < 0.0) {
rho *= -1.0;
sph *= -1.0;
}
cph = 0.0;
}
}
if (z == 0.0) {
if (!onz) {
r = rho;
cth = 0.0;
sth = 1.0;
return;
} else {
r = 0.0;
cth = 1.0;
sth = 0.0;
return;
}
} else {
if (!onz) {
r = sqrt(rho + z * z);
cth = z / r;
sth = rho / r;
} else {
r = z;
cth = 1.0;
sth = 0.0;
if (z < 0.0) {
r *= -1.0;
cth *= -1.0;
}
}
}
}
/*! \brief C++ porting of R3J000
*
* \param j2: `int`
* \param j3: `int`
* \param c6: `C6 *` Pointer to a C6 instance.
*/
void r3j000(int j2, int j3, C6 *c6) {
int jmx = j3 + j2;
if (jmx <= 0) {
c6->rac3j[0] = 1.0;
//printf("DEBUG: RAC3J(1) = %lE\n", c6->rac3j[0]);
return;
}
int jmn = j3 - j2;
if (jmn < 0) jmn *= -1;
int njmo = (jmx - jmn) / 2;
int jf = jmx + jmx + 1;
int isn = 1;
if (jmn % 2 != 0) isn = -1;
if (njmo <= 0) {
double sj = 1.0 * jf;
double cnr = (1 / sqrt(sj)) * isn;
c6->rac3j[0] = cnr;
//printf("DEBUG: RAC3J(1) = %lE\n", c6->rac3j[0]);
return;
}
double sjr = jf;
int jmxpos = (jmx + 1) * (jmx + 1);
int jmns = jmn * jmn;
int j1mo = jmx - 1;
int j1s = jmx * jmx; // a slight optimization was applied here
double cj = sqrt(1.0 * (jmxpos - j1s) * (j1s - jmns));
int j1mos = j1mo * j1mo;
double cjmo = sqrt(1.0 * (jmxpos - j1mos) * (j1mos - jmns));
if (njmo <= 1) {
c6->rac3j[0] = -cj / cjmo;
double sj = sjr + (c6->rac3j[0] * c6->rac3j[0]) * (jf - 4);
double cnr = (1.0 / sqrt(sj)) * isn;
c6->rac3j[1] = cnr;
c6->rac3j[0] *= cnr;
//printf("DEBUG: RAC3J(1) = %lE\n", c6->rac3j[0]);
//printf("DEBUG: RAC3J(2) = %lE\n", c6->rac3j[1]);
return;
}
int nj = njmo + 1;
int nmat = (nj + 1) / 2;
c6->rac3j[nj - 1] = 1.0;
c6->rac3j[njmo - 1] = -cj / cjmo;
if (nmat != njmo) {
int nbr = njmo - nmat;
for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
int irr = nj - ibr45;
jf = jf -4;
j1mo = j1mo - 2;
j1s = (j1mo + 1) * (j1mo + 1);
cj = sqrt(1.0 * (jmxpos - j1s) * (j1s - jmns));
j1mos = j1mo * j1mo;
cjmo = sqrt(1.0 * (jmxpos - j1mos) * (j1s - jmns));
c6->rac3j[irr - 2] = c6->rac3j[irr - 1] * (-cj / cjmo);
sjr = sjr + (c6->rac3j[irr - 1] * c6->rac3j[irr - 1]) * jf;
}
}
double racmat = c6->rac3j[nmat - 1];
sjr = sjr + (racmat * racmat) * (jf - 4);
c6->rac3j[0] = 1.0;
jf = jmn + jmn + 1;
double sjl = 1.0 * jf;
int j1pt = jmn + 2;
int j1pos = (j1pt - 1) * (j1pt - 1);
double cjpo = sqrt(1.0 * (jmxpos - j1pos) * (j1pos - jmns));
int j1pts = j1pt * j1pt;
double cjpt = sqrt(1.0 * (jmxpos - j1pts) * (j1pts - jmns));
c6->rac3j[1] = -cjpo / cjpt;
int nmatmo = nmat - 1;
if (nmatmo >= 2) {
for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
jf = jf + 4;
j1pt = j1pt + 2;
j1pos = (j1pt - 1) * (j1pt - 1);
cjpo = sqrt(1.0 * (jmxpos - j1pos) * (j1pos - jmns));
j1pts = j1pt * j1pt;
cjpt = sqrt(1.0 * (jmxpos - j1pts) * (j1pts - jmns));
c6->rac3j[irl70] = c6->rac3j[irl70 - 1] * (-cjpo / cjpt);
sjl = sjl + (c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1]) * jf;
}
}
double ratrac = racmat / c6->rac3j[nmat - 1];
double rats = ratrac * ratrac;
double sj = sjr + sjl * rats;
c6->rac3j[nmat - 1] = racmat;
double cnr = (1.0 / sqrt(sj)) * isn;
for (int irr80 = nmat; irr80 <= nj; irr80++) {
c6->rac3j[irr80 - 1] *= cnr;
}
double cnl = cnr * ratrac;
for (int irl85 = 1; irl85 <= nmatmo; irl85++) {
c6->rac3j[irl85 - 1] *= cnl;
}
}
/*! \brief C++ porting of STR
*
* \param c1: `C1 *` Pointer to a C1 instance.
* \param c1ao: `C1_AddOns *` Pointer to a C1_AddOns instance.
* \param c3: `C3 *` Pointer to a C3 instance.
* \param c4: `C4 *` Pointer to a C4 structure.
* \param c6: `C6 *` Pointer to a C6 instance.
*/
void str(C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
std::complex<double> *ylm;
c3 = new C3(); // this makes up for GCS = 0.0D0
int i = 0;
for (int l1po = 1; l1po <= c4->lmpo; l1po++) {
int l1 = l1po - 1;
for (int l2 = 1; l2 <= c4->lm; l2++) {
r3j000(l1, l2, c6);
c1ao->ind3j[l1po - 1][l2 - 1] = i;
int iabs = l2 - l1;
if (iabs < 0) iabs *= -1;
int lmnpo = iabs + 1;
int lmxpo = l2 + l1 + 1;
int il = 0;
int lpo = lmnpo;
while (lpo <= lmxpo) {
i++;
il++;
c1ao->v3j0[i - 1] = c6->rac3j[il - 1];
lpo += 2;
}
}
} // l1po loop
int nsphmo = c4->nsph - 1;
int lit = c4->li + c4-> li;
int lmt = c4->li + c4->le;
int size40 = nsphmo * c4->nsph * c4->litpos;
int size50 = c4->nsph * c4->lmtpos;
int ylm_size = (size40 > size50) ? size40 : size50;
ylm = new std::complex<double>[ylm_size];
int ivy = 0;
for (int nf40 = 1; nf40 <= nsphmo; nf40++) {
int nfpo = nf40 + 1;
for (int ns40 = nfpo; ns40 <= c4->nsph; ns40++) {
double rx = c1->rxx[nf40 - 1] - c1->rxx[ns40 - 1];
double ry = c1->ryy[nf40 - 1] - c1->ryy[ns40 - 1];
double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
double rr, crth, srth, crph, srph;
polar(rx, ry, rz, rr, crth, srth, crph, srph);
//printf("DEBUG: crth = %lE\n", crth);
//printf("DEBUG: srth = %lE\n", srth);
//printf("DEBUG: crph = %lE\n", crph);
//printf("DEBUG: srph = %lE\n", srph);
sphar(crth, srth, crph, srph, lit, ylm);
for (int iv38 = 1; iv38 <= c4->litpos; iv38++) {
c1ao->vyhj[iv38 + ivy - 1] = dconjg(ylm[iv38 - 1]);
} // iv38 loop
ivy += c4->litpos;
} // ns40 loop
} // nf40 loop
ivy = 0;
for (int nf50 = 1; nf50 <= c4->nsph; nf50++) {
double rx = c1->rxx[nf50 - 1];
double ry = c1->ryy[nf50 - 1];
double rz = c1->rzz[nf50 - 1];
double rr, crth, srth, crph, srph;
if (rx == 0.0 && ry == 0.0 && rz == 0.0) continue;
polar(rx, ry, rz, rr, crth, srth, crph, srph);
sphar(crth, srth, crph, srph, lmt, ylm);
for (int iv48 = 1; iv48 <= c4->lmtpos; iv48++) {
c1ao->vyj0[iv48 + ivy - 1] = dconjg(ylm[iv48 - 1]);
} // iv48 loop
ivy += c4->lmtpos;
} // nf50 loop
}
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment