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.

Your formulas are confused: if `U` is `MxM` and `V` is `NxN` then `Vck*Ulc` doesn't make sense because for the first term the `c` range is `N` and for the second term it should be instead `M` and given that there are two names `N` and `M` one should assume that the two values may be different. — 6502, Mar 04 '14 at 20:33
I addition to the inconsistent matrix dimensions, you should check where you set and reset you sum counter, namely around the `c` loop. Set `thesum` to zero in the `l` loop or better yet, make `thesum` a local variable in the `l` loop scope and initialise it to zero. Assignment to `result` could be placed better after the `c` loop. — M Oehm, Mar 05 '14 at 08:07
@M Oehm:I had some mistakes indeed.Thank you for your suggestions. — George, Mar 05 '14 at 19:50

score 1 · Accepted Answer · answered Mar 04 '14 at 20:36

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).

multiplication of 3 matrices in a sum

1 Answers1

Linked