Matrix Multiplication Error in C

Question

I have been working on some code to implement a discrete cosine transform using C. This is the code I am using for my matrix multiplication.

float Yt[8][8][4];
float Y8[8][8][4]={{{0.980785,-0.555570,0.555570 ...}}};
float A[8][8]={{0.707107,0.707107 ... }};

for(k=0;k<4;++k){
        for(i=0;i<8;++i){
            for(j=0;j<8;++j){
            Yt[i][j][k]=0;
            }}}

for(k=0;k<4;++k){
    for(i=0;i<8;++i){
         for(j=0;j<8;++j){
            for(l=0;l<8;++l){
                Yt[i][j][k]+=A[i][l]*Y8[l][j][k];
                }}}}

Y8 is the same 8x8 matrix, repeated 4 times. My code successfully compiles and runs. However, the output that I receive is obviously incorrect. Yt[i][j][k] for k=0,1,2,3 should all be the same value, but I end up with the correct values for k=0,1 and different values for k=2,3.

Can anyone see why I am getting incorrect values for the last two matrices?

Answer 1

@chux makes a good call about initializing they value of Yt before using it. You might do something like this:

for(k=0;k<4;++k){
  for(i=0;i<8;++i){
    for(j=0;j<8;++j){
      double temp = 0.0;
      for(l=0;l<8;++l){
        temp += A[i][l]*Y8[l][j][k];
      }
      Yt[i][j][k] = temp;
    }
  }
}

This would allow the compiler to use a register for the accumulation, then access Yt (with its three indexing operations) just once.

Answer 2

Initialize Yt[i][j][k].

I suspect this would be done as below

for(j=0;j<8;++j){
  Yt[i][j][k] = 0.0;
  for(l=0;l<8;++l){
    Yt[i][j][k] += A[i][l]*Y8[l][j][k];
  }
}

... or in its declaration. Something like double Yt[8][8][4] = { 0.0 };