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