- 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];
+ 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;
+ if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow(pi[r],gamma))
+ phi[ai(j,mm,r,p,m,k)] = 0;
+ else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow(pi[r],gamma))
+ phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow(pi[r],gamma) - S[ai(j,mm,r,p,m,k)])
+ / Gram2[ai(j,j,r,p,p,k)];