multiplication of 3 matrices in a sum

Question

I want to compute:

result = SUM (c=0,N) { V_ck * U_lc * S_c }

But my 2D matrices are indexed as 1D and stored as column major.

I am trying:

float *A,*B;

int M,N;


A = (float *) malloc( M * N * sizeof(float) );
B = (float *) malloc( M * sizeof(float) );


float *S = (float *)malloc( N * sizeof( float) ); 
float *U = (float *)malloc( M * M * sizeof( float) ); 
float *V = (float *)malloc( N * N  * sizeof( float) ); 

float * result = (float *) malloc( M * N * sizeof(float) );

float thesum=0;

for (int k=0; k < M ; k++) {

      for (int l=0 ; l < N; l++) {

        for ( int c=0; c < N; c++) {

        thesum += V[ k+c*N ] * S[c] * U[l*M + k];
        result[ k+l*M ]=thesum;

                }
        }

        }

I have one big mistake ,I think in the above,because I need another loop ? in order to do properly the multiplication first and then use the:

for ( int c=0; c < N; c++)

Loop for the summation,right?

And then, do I have to create an array which will hold the multiplication values and then use this array to hold the summation values?

If I were using 2D notation , I would simple use U[l][k] and so on.

But now,I am confused about how to apply the appropriate indices.

And I want someone to explain me how should I proceed with this.

Answer 1

If I were using 2D notation , I would simple use U[l][k] and so on.

So, add that layer of abstraction - don't let everything else get complicated. You have:

A = (float *) malloc( M * N * sizeof(float) );

At a minimum, you can use:

float& at(float* p, int rows, int col, int row) { return p[rows * col + row]; }

(reorder arguments to taste)

Then you can say:

at(A, M, col, row)

(Or similar - I wouldn't swear I got all the rows/column names right - but IMHO you should have used Rows and Columns instead of M and N so I'm not going to bust a gut over it.)

If you want to get a little fancier, in C++ you can wrap the allocations in a class that stores the pointer and #rows/columns, then overloads float& operator()(int col, int row) and const float& operator()(int col, int row) const (or just float operator()(int col, int row) const if you don't care about ability to take the address of the array entry).