Matrix multiplication

[dukebball]

Hey, so I have to make a program that has to multiply two 3x3 matrices (among other things). I used classes to do the project, but I just can't seem to get this multiplication to work. Here is my code so far:

void matrix::mult(const matrix& u)
{
    double c[9];
    double a[9];
    for (int i=0; i<rows; i++)
    {
        for (int j=0; j<columns; j++)
        {
            for(int r=0; r<columns; r++)
            {
                c[i*columns+j]+=a[i*columns+r]*u.data[r*columns+j];
            }
        }
    }
}

Any advice would be great, and please let me know if there's anything I should clarify. Thank You.

[Lindley]

You don't seem to be initializing the c array prior to using it, and I have no idea what you think the a array contains, because it appears to be uninitialized as well from that code.

[deightm]

void matrix::mult(const matrix& u)
{
unsigned char i=0,j=0,k=0;

double c[9];
for(i=0;i<9;i++) c[i] = 0;

for(i=0;i<3;i++)
for(j=0;j<3);j++)
for(k=0;k<3;k++)
c[i*3+j] += this->data[i*columns+j] * u.data[k*columns+j];
}

This should set c = "this matrix" * u, assuming that "data" is double[9] array which stores matrix value;

[bitshifter420]

Rule #1) Use single dimensional array for speed and clarity.
Rule #2) Dont use a loop, even if your compiler can unroll them.
Rule #3) Use SIMD for this type of operation (do 4 muls at once)

Tomorrow when i get to my other PC i will have some code to share.
I have recently written an entire matrix library in assembly language for SSE2.
This library was developed by me and a friend and each procedure
was carefully coded and clocked using high resolution uops counters.
If you dont like assembly language then use intrinsics instead.

Edit:
Ooh, i found a copy of my C version right here on this HDD

void MultiplyMatrices(float dst[16], const float a[16], const float b[16])
{
   dst[0]  = a[0] * b[0]  + a[4] * b[1]  + a[8]  * b[2]  + a[12] * b[3];
   dst[1]  = a[1] * b[0]  + a[5] * b[1]  + a[9]  * b[2]  + a[13] * b[3];
   dst[2]  = a[2] * b[0]  + a[6] * b[1]  + a[10] * b[2]  + a[14] * b[3];
   dst[3]  = a[3] * b[0]  + a[7] * b[1]  + a[11] * b[2]  + a[15] * b[3];

   dst[4]  = a[0] * b[4]  + a[4] * b[5]  + a[8]  * b[6]  + a[12] * b[7];
   dst[5]  = a[1] * b[4]  + a[5] * b[5]  + a[9]  * b[6]  + a[13] * b[7];
   dst[6]  = a[2] * b[4]  + a[6] * b[5]  + a[10] * b[6]  + a[14] * b[7];
   dst[7]  = a[3] * b[4]  + a[7] * b[5]  + a[11] * b[6]  + a[15] * b[7];

   dst[8]  = a[0] * b[8]  + a[4] * b[9]  + a[8]  * b[10] + a[12] * b[11];
   dst[9]  = a[1] * b[8]  + a[5] * b[9]  + a[9]  * b[10] + a[13] * b[11];
   dst[10] = a[2] * b[8]  + a[6] * b[9]  + a[10] * b[10] + a[14] * b[11];
   dst[11] = a[3] * b[8]  + a[7] * b[9]  + a[11] * b[10] + a[15] * b[11];

   dst[12] = a[0] * b[12] + a[4] * b[13] + a[8]  * b[14] + a[12] * b[15];
   dst[13] = a[1] * b[12] + a[5] * b[13] + a[9]  * b[14] + a[13] * b[15];
   dst[14] = a[2] * b[12] + a[6] * b[13] + a[10] * b[14] + a[14] * b[15];
   dst[15] = a[3] * b[12] + a[7] * b[13] + a[11] * b[14] + a[15] * b[15];
}

If you want the SSE2 assembly code version just let me know...

[potatoCode]

never mind.

[Lindley]

See, if I needed that degree of optimization, I'd just use one of the optimized BLAS implementations.