Loading src/libnptm/clu_subs.cpp +1 −127 Original line number Diff line number Diff line Loading @@ -821,9 +821,7 @@ dcomplex ghit( int lmaxpo = l2 + l1mp + 1; int i3j0in = c1->ind3j[l1mp][l2 - 1]; int ilin = (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ? 0 : -1; // if (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ilin = 0; int isn = (m1 % 2 != 0) ? -1 : 1; // if (m1 % 2 != 0) isn *= -1; if (lminpo % 2 == 0) { isn *= -1; if (l2 > l1mp) isn *= -1; Loading Loading @@ -1557,7 +1555,7 @@ void r3j000(int j2, int j3, double *rac3j) { } } void r3jjr(int j2, int j3, int m2, int m3, double *rac3j) { void r3jjr(int j2, int j3, int m2, int m3, double rac3j[]) { int jmx = j3 + j2; int jdf = j3 - j2; int m1 = -m2 - m3; Loading Loading @@ -1673,130 +1671,6 @@ void r3jjr(int j2, int j3, int m2, int m3, double *rac3j) { } } // #ifdef USE_TARGET_OFFLOAD // #pragma omp begin declare target device_type(any) // #endif void r3jjr_d(int j2, int j3, int m2, int m3, double *rac3j) { int jmx = j3 + j2; int jdf = j3 - j2; int m1 = -m2 - m3; int abs_jdf = (jdf >= 0) ? jdf : -jdf; int abs_m1 = (m1 >= 0) ? m1 : -m1; int jmn = (abs_jdf > abs_m1) ? abs_jdf : abs_m1; int njmo = jmx - jmn; int jf = jmx + jmx + 1; int isn = 1; if ((jdf + m1) % 2 != 0) isn = -1; if (njmo <= 0) { double sj = 1.0 * jf; double cnr = (1.0 / sqrt(sj)) * isn; rac3j[0] = cnr; } else { // label 15 double sjt = 1.0; double sjr = 1.0 * jf; int jsmpos = (jmx + 1) * (jmx + 1); int jdfs = jdf * jdf; int m1s = m1 * m1; int mdf = m3 - m2; int idjc = m1 * (j3 * (j3 + 1) - j2 * (j2 +1)); int j1 = jmx; int j1s = j1 * j1; int j1po = j1 + 1; double ccj = 1.0 * (j1s - jdfs) * (j1s - m1s); double cj = sqrt(ccj * (jsmpos - j1s)); double dj = 1.0 * jf * (j1 * j1po * mdf + idjc); if (njmo <= 1) { rac3j[0] = -dj / (cj * j1po); double sj = sjr + (rac3j[0] * rac3j[0]) * (jf - 2); double cnr = (1.0 / sqrt(sj)) * isn; rac3j[1] = cnr; rac3j[0] *= cnr; } else { // label 20 double cjp = 0.0; int nj = njmo + 1; int nmat = (nj + 1) / 2; rac3j[nj - 1] = 1.0; rac3j[njmo - 1] = -dj / (cj * j1po); if (nmat != njmo) { int nbr = njmo - nmat; for (int ibr45 = 1; ibr45 <= nbr; ibr45++) { int irr = nj - ibr45; jf -= 2; j1--; j1s = j1 * j1; j1po = j1 + 1; cjp = cj; ccj = 1.0 * (j1s - jdfs) * (j1s - m1s); cj = sqrt(ccj * (jsmpos - j1s)); sjt = rac3j[irr - 1] * rac3j[irr - 1]; dj = 1.0 * jf * (j1 * j1po * mdf + idjc); rac3j[irr - 2] = -(rac3j[irr - 1] * dj + rac3j[irr] * cjp * j1) / (cj * j1po); sjr += (sjt * jf); } // ibr45 loop } // label 50 double osjt = sjt; sjt = rac3j[nmat - 1] * rac3j[nmat - 1]; if (sjt >= osjt) { sjr += (sjt * (jf - 2)); } else { // label 55 nmat++; } // label 60 double racmat = rac3j[nmat - 1]; rac3j[0] = 1.0; jf = jmn + jmn + 1; double sjl = 1.0 * jf; j1 = jmn; if (j1 != 0) { j1po = j1 + 1; int j1pos = j1po * j1po; double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s); cjp = sqrt(ccjp * (jsmpos - j1pos)); dj = 1.0 * jf * (j1 * j1po * mdf + idjc); rac3j[1] = - dj / (cjp * j1); } else { // label 62 cjp = sqrt(1.0 * (jsmpos - 1)); dj = 1.0 * mdf; rac3j[1] = -dj / cjp; } // label 63 int nmatmo = nmat - 1; if (nmatmo >= 2) { for (int irl70 = 2; irl70 <= nmatmo; irl70++) { jf += 2; j1++; j1po = j1 + 1; int j1pos = j1po * j1po; cj = cjp; double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s); cjp = sqrt(ccjp * (jsmpos - j1pos)); sjt = rac3j[irl70 - 1] * rac3j[irl70 - 1]; dj = 1.0 * jf * (j1 * j1po * mdf + idjc); rac3j[irl70] = -( rac3j[irl70 - 1] * dj + rac3j[irl70 - 2] * cj * j1po ) / (cjp * j1); sjl += (sjt * jf); } } // label 75 double ratrac = racmat / rac3j[nmat - 1]; double rats = ratrac * ratrac; double sj = sjr + sjl * rats; rac3j[nmat - 1] = racmat; double cnr = (1.0 / sqrt(sj)) * isn; for (int irr80 = nmat; irr80 <= nj; irr80++) rac3j[irr80 - 1] *= cnr; double cnl = cnr * ratrac; for (int irl85 = 1; irl85 <= nmatmo; irl85++) rac3j[irl85 - 1] *= cnl; } } } // #ifdef USE_TARGET_OFFLOAD // #pragma omp end declare target // #endif void r3jmr(int j1, int j2, int j3, int m1, double *rac3j) { int mmx = (j2 < j3 - m1) ? j2 : j3 - m1; int mmn = (-j2 > -(j3 + m1)) ? -j2 : -(j3 + m1); Loading Loading
src/libnptm/clu_subs.cpp +1 −127 Original line number Diff line number Diff line Loading @@ -821,9 +821,7 @@ dcomplex ghit( int lmaxpo = l2 + l1mp + 1; int i3j0in = c1->ind3j[l1mp][l2 - 1]; int ilin = (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ? 0 : -1; // if (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ilin = 0; int isn = (m1 % 2 != 0) ? -1 : 1; // if (m1 % 2 != 0) isn *= -1; if (lminpo % 2 == 0) { isn *= -1; if (l2 > l1mp) isn *= -1; Loading Loading @@ -1557,7 +1555,7 @@ void r3j000(int j2, int j3, double *rac3j) { } } void r3jjr(int j2, int j3, int m2, int m3, double *rac3j) { void r3jjr(int j2, int j3, int m2, int m3, double rac3j[]) { int jmx = j3 + j2; int jdf = j3 - j2; int m1 = -m2 - m3; Loading Loading @@ -1673,130 +1671,6 @@ void r3jjr(int j2, int j3, int m2, int m3, double *rac3j) { } } // #ifdef USE_TARGET_OFFLOAD // #pragma omp begin declare target device_type(any) // #endif void r3jjr_d(int j2, int j3, int m2, int m3, double *rac3j) { int jmx = j3 + j2; int jdf = j3 - j2; int m1 = -m2 - m3; int abs_jdf = (jdf >= 0) ? jdf : -jdf; int abs_m1 = (m1 >= 0) ? m1 : -m1; int jmn = (abs_jdf > abs_m1) ? abs_jdf : abs_m1; int njmo = jmx - jmn; int jf = jmx + jmx + 1; int isn = 1; if ((jdf + m1) % 2 != 0) isn = -1; if (njmo <= 0) { double sj = 1.0 * jf; double cnr = (1.0 / sqrt(sj)) * isn; rac3j[0] = cnr; } else { // label 15 double sjt = 1.0; double sjr = 1.0 * jf; int jsmpos = (jmx + 1) * (jmx + 1); int jdfs = jdf * jdf; int m1s = m1 * m1; int mdf = m3 - m2; int idjc = m1 * (j3 * (j3 + 1) - j2 * (j2 +1)); int j1 = jmx; int j1s = j1 * j1; int j1po = j1 + 1; double ccj = 1.0 * (j1s - jdfs) * (j1s - m1s); double cj = sqrt(ccj * (jsmpos - j1s)); double dj = 1.0 * jf * (j1 * j1po * mdf + idjc); if (njmo <= 1) { rac3j[0] = -dj / (cj * j1po); double sj = sjr + (rac3j[0] * rac3j[0]) * (jf - 2); double cnr = (1.0 / sqrt(sj)) * isn; rac3j[1] = cnr; rac3j[0] *= cnr; } else { // label 20 double cjp = 0.0; int nj = njmo + 1; int nmat = (nj + 1) / 2; rac3j[nj - 1] = 1.0; rac3j[njmo - 1] = -dj / (cj * j1po); if (nmat != njmo) { int nbr = njmo - nmat; for (int ibr45 = 1; ibr45 <= nbr; ibr45++) { int irr = nj - ibr45; jf -= 2; j1--; j1s = j1 * j1; j1po = j1 + 1; cjp = cj; ccj = 1.0 * (j1s - jdfs) * (j1s - m1s); cj = sqrt(ccj * (jsmpos - j1s)); sjt = rac3j[irr - 1] * rac3j[irr - 1]; dj = 1.0 * jf * (j1 * j1po * mdf + idjc); rac3j[irr - 2] = -(rac3j[irr - 1] * dj + rac3j[irr] * cjp * j1) / (cj * j1po); sjr += (sjt * jf); } // ibr45 loop } // label 50 double osjt = sjt; sjt = rac3j[nmat - 1] * rac3j[nmat - 1]; if (sjt >= osjt) { sjr += (sjt * (jf - 2)); } else { // label 55 nmat++; } // label 60 double racmat = rac3j[nmat - 1]; rac3j[0] = 1.0; jf = jmn + jmn + 1; double sjl = 1.0 * jf; j1 = jmn; if (j1 != 0) { j1po = j1 + 1; int j1pos = j1po * j1po; double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s); cjp = sqrt(ccjp * (jsmpos - j1pos)); dj = 1.0 * jf * (j1 * j1po * mdf + idjc); rac3j[1] = - dj / (cjp * j1); } else { // label 62 cjp = sqrt(1.0 * (jsmpos - 1)); dj = 1.0 * mdf; rac3j[1] = -dj / cjp; } // label 63 int nmatmo = nmat - 1; if (nmatmo >= 2) { for (int irl70 = 2; irl70 <= nmatmo; irl70++) { jf += 2; j1++; j1po = j1 + 1; int j1pos = j1po * j1po; cj = cjp; double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s); cjp = sqrt(ccjp * (jsmpos - j1pos)); sjt = rac3j[irl70 - 1] * rac3j[irl70 - 1]; dj = 1.0 * jf * (j1 * j1po * mdf + idjc); rac3j[irl70] = -( rac3j[irl70 - 1] * dj + rac3j[irl70 - 2] * cj * j1po ) / (cjp * j1); sjl += (sjt * jf); } } // label 75 double ratrac = racmat / rac3j[nmat - 1]; double rats = ratrac * ratrac; double sj = sjr + sjl * rats; rac3j[nmat - 1] = racmat; double cnr = (1.0 / sqrt(sj)) * isn; for (int irr80 = nmat; irr80 <= nj; irr80++) rac3j[irr80 - 1] *= cnr; double cnl = cnr * ratrac; for (int irl85 = 1; irl85 <= nmatmo; irl85++) rac3j[irl85 - 1] *= cnl; } } } // #ifdef USE_TARGET_OFFLOAD // #pragma omp end declare target // #endif void r3jmr(int j1, int j2, int j3, int m1, double *rac3j) { int mmx = (j2 < j3 - m1) ? j2 : j3 - m1; int mmn = (-j2 > -(j3 + m1)) ? -j2 : -(j3 + m1); Loading