- for (mwSize u=0; u<j; u++)
- dotPhiGram2 += phi[u*m*k+mm*k+r] * Gram2[j*p*k+u*k+r];
- for (mwSize u=j+1; u<p; u++)
- dotPhiGram2 += phi[u*m*k+mm*k+r] * Gram2[j*p*k+u*k+r];
- //S(j,r,mm)=-rho(mm,mm,r)*ps2(j,mm,r)+sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))
- // +sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r)));
- S[j*m*k+mm*k+r] = -rho[mm*m*k+mm*k+r] * ps2[j*m*k+mm*k+r] + dotPhiGram2;
- if (fabs(S[j*m*k+mm*k+r]) <= n*lambda*pow(pi[r],gamma))
- phi[j*m*k+mm*k+r] = 0;
- else if (S[j*m*k+mm*k+r] > n*lambda*pow(pi[r],gamma))
- phi[j*m*k+mm*k+r] = (n*lambda*pow(pi[r],gamma) - S[j*m*k+mm*k+r])
- / Gram2[j*p*k+j*k+r];
+ 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)
+ {
+ phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow_pir_gamma - S[ai(j,mm,r,p,m,k)])
+ / Gram2[ai(j,j,r,p,p,k)];
+ }