- //sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r])
- Real dotPhiGram2 = 0.0;
- for (int u=0; u<j; u++)
- dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
- //sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r])
- for (int u=j+1; u<p; u++)
- dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
- //S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] +
- // (if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) +
- // (if(j<p) sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r]) else 0)
- S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)] + dotPhiGram2;
- Real pow_pir_gamma = pow(pi[r],gamma);
- if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow_pir_gamma)
- phi[ai(j,mm,r,p,m,k)] = 0;
- else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow_pir_gamma)
+ //sum(phi[-j,mm,r] * Gram2[j,-j,r])
+ Real phiDotGram2 = 0.;
+ for (int u=0; u<p; u++)
+ {
+ if (u != j)
+ phiDotGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
+ }
+ //S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r])
+ S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)]
+ + phiDotGram2;
+ Real pirPowGamma = pow(pi[r],gamma);
+ if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pirPowGamma)
+ phi[ai(j,mm,r,p,m,k)] = 0.;
+ else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pirPowGamma)