Function to multiply 3x3 matrices gives wrong answer for middle column only

Question

While teaching myself c, I thought it would be good practice to write a function which multiplies two 3x3 matrices and then make it more general. The function seems to calculate the correct result for the first and last columns but not the middle one. In addition, each value down the middle column is out by 3 more than the last.

For example:

[1 2 3]   [23  4  6]
[4 5 6] * [ 2 35  0]
[7 8 9]   [14  2 43]

The answer I receive is:

[ 69  80 135]
[190 273 282]
[303 326 429]

The actual answer should be:

[ 69  83 135]
[190 279 282]
[303 335 429]

Isolating the middle columns for clarity:

Received  Expected
   [ 80]  [ 83]
   [273]  [279]
   [326]  [335]

My code is as follows:

#include <stdio.h>

typedef struct mat_3x3
{
    double values [3][3];
} mat_3x3;



void SetMatrix(mat_3x3 * matrix, double vals[3][3])
{
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            (matrix->values)[i][j] = vals[i][j];
        }
    }
    putchar('\n');
}



mat_3x3 MatrixMultiply(mat_3x3 * m1, mat_3x3 * m2)
{
    mat_3x3 result;
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            double temp = 0;
            for (int k = 0; k < 3; k++)
            {               
                temp += ((m1->values)[i][k] * (m2->values)[k][j]);
            }
            (result.values)[i][j] = temp;
        }
    }
    return result;
}



void PrintMatrix(mat_3x3 * matrix)
{
    putchar('\n');
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 3; j++)
        {
            printf("%lf  ", (matrix->values)[i][j]);
        }
        putchar('\n');
    }
    putchar('\n');
}



int main()
{
    mat_3x3 m1;
    mat_3x3 * pm1 = &m1;
    
    mat_3x3 m2;
    mat_3x3 * pm2 = &m2;
    
    double vals[3][3] = {
        {1,2,3},
        {4,7,6},
        {7,8,9}
    };
    
    double vals2[3][3] = {
        {23,4,6},
        {2,35,0},
        {14,2,43}
    };
    
    SetMatrix(pm1, vals);
    SetMatrix(pm2, vals2);
    
    printf("\nm1:");
    PrintMatrix(pm1);
    printf("\nm2:");
    PrintMatrix(pm2);
        
    mat_3x3 m3 = MatrixMultiply(pm1, pm2);
    mat_3x3 * pm3 = &m3;
    printf("\nm3 = m1 * m2");
    PrintMatrix(pm3);
    
}

Have been working on this for a while now comparing it against other simple examples and can't find the problem, so help would be appreciated!

Also if I've done anything atrocious syntax wise etc, I'm open to any criticism on how it's written as well.

Answer 1

While teaching myself c, I thought it would be good practice to write a function which multiplies two 3x3 matrices and then make it more general. The function seems to calculate the correct result for the first and last columns but not the middle one. In addition, each value down the middle column is out by 3 more than the last.

In practice, when coding in C, you should take care of the following issues:

• refer to a good C reference website and read a good C programming book, such as Modern C

• floating point numbers are not mathematical real numbers, see floating-point-gui.de for much more. For example, addition is associative in math, but not on a computer using IEEE-754.

• we all make bugs (e.g. buffer overflows or undefined behavior). So you need to learn how to use a debugger. I recommend GDB. But you need to learn how to use it and spend a few hours reading documentation. Tools like valgrind are also useful (to hunt memory leaks) as soon as you use C dynamic memory allocation.

• recent compilers can be helpful. I recommend GCC. You should invoke it with all warnings and debug info, e.g. gcc -Wall -Wextra -g. Be sure to spend some time in reading the documentation of your compiler. You might later consider using static program analysis tools such as Frama-C or the Clang analyzer or (for precision analysis) Fluctuat or CADNA

• consider having a matrix abstract data type like here. You would then easily generalize your code to "arbitrary" N*M matrixes.

• later, for benchmarking purposes, you will want to use an optimizing compiler. If you use GCC, you could compile your code using gcc -Wall -Wextra -g -O3 but then you could have surprising optimizations, see e.g. this draft report.

• in some cases, you could need arbitrary-precision arithmetic. Consider then using specialized libraries such as GMPlib.

• Most computers today are multi-core. You could want to use Pthreads or MPI to take advantage of that with concurrent programming.

• many open source libraries exist for scientific computations. Look at least for inspiration on github and gitlab and see also this list. You could be interested by GNU GSL and study its source code since it is free software (and later improve it).

If you want to make serious scientific computations, you might consider switching (for expressiveness) to functional languages such as Ocaml. If you care about making a lot of iterative computing (like in finite element methods) you might switch to OpenCL or OpenACC.

Be aware that scientific computation is a very difficult field.

Expect to spend a decade in learning it.

---

I'm open to any criticism on how it's written as well.

 mat_3x3 MatrixMultiply(mat_3x3 * m1, mat_3x3 * m2)

is unusual. Why don't you return a pointer (to a fresh memory zone obtained with malloc and correctly initialized) ? That is likely to be faster (a pointer is usually 8 bytes, a 3x3 matrix takes 72 bytes to be copied) and enable you to code things like MatrixMultiply(MatrixMultiply(M1, M2), MatrixAdd(M2, M3)). Of course, garbage collection (read the GC handbook, consider using Boehm GC) then becomes an issue. If you used Ocaml, the system GC would be very helpful.