Improve nested loop for bio-statistics calculation I am doing a bio-statistics calculation and the following code works. However, can someone help to improve the messy nested loop?

for(int i=0; i&lt;NN; i++) {
    for (int j=0; j&lt;NN; j++) {
        if (i != j){
            thirdlayer = 0;
            for (int k=0; k&lt;NN; k++) {
                fourthlayer = 0;
                for (int l=0; l&lt;NN; l++) {
                    fourthlayer =  fourthlayer + V[j*NN+l]*V[NN+l]*J[k*NN+l];
                }
                thirdlayer = thirdlayer + V[k]*V[i*NN+k]*fourthlayer;
            }
            if(pi_cod[j] != 0)
                Transitions[i*NN +j] =  sqrt(pi_cod[i]*pi_cod[1]/(pi_cod[0]*pi_cod[j]))*Q[i*NN +j]*thirdlayer/Padt;
        }
    }
}
Using small l for the index is bad because it looks like the digit 1.  It's better to use large L.

Instead of j*NN, it's better to use a cached index that increments by for:

int NN2 = NN*NN
for(int iNN=k; iNN &lt; (NN2+k); iNN+=NN) {
  thirdlayer = thirdlayer + V[k]*V[iNN]*fourthlayer;
}


This could be a bit faster.

Another hook - more use pointer as array sintax:
for 3layer better get a row vector in wich 4layer for process:

int* VVj = &amp;(V[j*NN]);
int* VNN = &amp;(V[NN]);
for (int k=0; k&lt;NN; k++) {
                int* JNNk = &amp;(J[k*NN]);
                fourthlayer = 0;
                for (int l=0; l&lt;NN; l++) {
                    fourthlayer =  fourthlayer + VVj[l]*VNN[l]*JNNk[l];
                }
                thirdlayer = thirdlayer + V[k]*V[i*NN+k]*fourthlayer;
            }


good compiler do it for you byself, but in such decomposition may better see data dependents, and it is a bit simpler and short

also you can deploy from 4layer V[j*NN+l]*V[NN+l] into stanalone vector that can be prepared in outer j cycle.

Instead of division (xxx)/Padt (better=faster), use multiplication *(1/Padt), or move out from the last cycle:

double thp = thirdlayer/Padt;
for(int i=0; i&lt;NN; i++) { 
    if(i != j &amp;&amp; pi_cod[j] != 0)
      Transitions[i*NN +j] =  sqrt(pi_cod[i]*pi_cod[1]/(pi_cod[0]*pi_cod[j]))*Q[i*NN +j]*thp;
  }


Instead of a conditional calculation, it's better to use a conditional assignment since it can better optimized for x86:

double thp = thirdlayer/Padt;
for(int i=0; i&lt;NN; i++) {
    int pcodj = (pi_cod[j] != 0)?pi_cod[j]: 1;
    double transition = sqrt(pi_cod[i]*pi_cod[1]/(pi_cod[0]*pi_cod[j]))*Q[i*NN +j]*thp;
    if(i != j &amp;&amp; pi_cod[j] != 0)
      Transitions[i*NN +j] = transition;
  }


If it's rare misses for assignment, so penalty for calculation could be negligible.