Inaccuracies in matrix multiplication- C

Question

I'm doing a DS project in C. I've been struggling with a matrix multiplication function:

void mat_mult(double **mat1, double **mat2, double **res, int n) {
    int i, j, k;
    
    for (i = 0; i < n; i++) {
        
        for (j = 0; j < n; j++) {
            assert(res[i][j] == 0);
            for (k = 0; k < n; k++)
                res[i][j] += mat1[i][k] * mat2[k][j];;
        }
    }
} 

The matrices I'm trying to multiply are:

0.00000000, 0.00007291, 0.00000000 
0.00007291, 0.00000000, 0.00000001 
0.00000000, 0.00000001, 0.00000000

and:

117.11258762, 0.00000000, 0.00000000 
0.00000000, 117.10123871, 0.00000000 
0.00000000, 0.00000000, 8087.59220061

The result is:

0.00000000, 0.00853872, 0.00000007 
0.00853790, 0.00000000, 0.00000172 
0.00000467, 0.00011897, 0.00000000 

And not:

0.00000000 0.00853868 0.00000000
0.00853785 0.00000000 0.00000117
0.00000000 0.00008088 0.00000000

So, the inaccuracies are very small but could be significant. Does anyone have an idea on how to solve or what to look for in trying to solve this?

***Edit- generating code ***

int main(){
int i,j;
int dim=3;
int n=3;
double* a;
double** wam_res, **ddg_res, **lnorm_res, **test_res;
double** vecs;
vecs= (double**)calloc(n, sizeof(double*));
assert(vecs!=NULL);

for(i=0; i<n; i++)
{
    a= (double*)calloc(dim,sizeof(double));
    assert(a!=NULL);
    for(j=0; j<dim; j++)
    {
        a[j]= rand() %50;
    }
    vecs[i]= a;
    
}
wam_res= wam(vecs,dim,n);
printf("%s","wam:\n");
print_matrix(wam_res, n);
ddg_res= ddg(wam_res, n);
printf("%s","ddg:\n");
print_sym_mat(ddg_res, n); 

lnorm_res= lnorm(ddg_res, wam_res, n);
printf("%s","lnorm:\n"); }


double** wam(double** vecs, int dim, int n) {
 /*Calculating wam for the given set of vectors- n vectors of dimension dim */
double** res;
int i,j, curr_index=0;
double norm;
res= (double**)calloc(n, sizeof(double*));
for(i=0; i< n; i++)
{
    res[i]= (double*)calloc(n, sizeof(double));

    /*For vector i, calculating wij for all 0<=j<n*/
    for(j=0; j<n; j++)
    {
        if(i!=j)

        {
            
            norm= l2_norm2vec(*(vecs+i),*(vecs+j),dim); /* function for l2 norm */
            
            res[i][j]= exp((-1*norm/2));
            
            curr_index++;
        }
        else
        {
            /* no self loops */
            res[i][j]=0;
        }
    }
}
return res; }

double** ddg(double** wam, int n) {
int i, j; 
double** res;
double d=0; 
res= (double**)calloc(n*n, sizeof(double*)); /* a n*n array with just 0s */
for(i=0; i<n; i++)
{
    res[i]= (double*)calloc(n, sizeof(double));
    for(j=0; j<n; j++)
    {
        d+=wam[i][j];
    }
    res[i][i]= 1/sqrt(d);
    d=0;

}
return res; }

double** lnorm(double** ddg, double** wam,  int n) {
double** res, **tmp;
int i;
res= (double**)calloc(n, sizeof(double));
tmp= (double**)calloc(n, sizeof(double));
for(i=0; i<n; i++)
{
    res[i]= (double*)calloc(n, sizeof(double));
    tmp[i]= (double*)calloc(n, sizeof(double));
}

mat_mult(ddg,wam, res, n);
printf("%s", "step 1 results: \n");
print_sym_mat(res,n); /* errors happened here */
mat_mult(copyMatrix(res,n), ddg,tmp,n);  /*copyMatrix- simple loop to allocate an identical matrix */
res= copyMatrix(tmp,n);
printf("%s", "step 2 results: \n");
print_sym_mat(res,n); /* errors happened here too */
mat_subI(res,n); 
printf("%s", "step 3 results: \n");
print_sym_mat(res,n);
return res; }

The main function that calls the matrix multiplication is lnorm, that appears at the end of the code. It receives two vector calculations for different graph representations

Answer 1

Inaccuracies are not to blame here: the second matrix is a diagonal matrix so the matrix multiplication involves a single multiplication per element of the resulting matrix. Your code outputs only 8 decimal places, whereas the values have more data.

Here is a modified version of the code that compiles and outputs the posted results, exhibiting the issue: the values in the matrices are rounded to 8 decimal places, which is a significant difference for some of the elements:

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double **wam(double **vecs, int dim, int n);
double **ddg(double **wam, int n);
double **lnorm(double **ddg, double **wam, int n);

double l2_norm2vec(double *vec1, double *vec2, int n) {
    double res = 0;
    for (int i = 0; i < n; i++) {
        double d = vec2[i] - vec1[i];
        res += d * d;
    }
    return sqrt(res);
}

#define print_sym_mat print_matrix

double **copyMatrix(double **res, int n) {
    double **mat = malloc(sizeof(*mat) * n);
    for (int i = 0; i < n; i++) {
        mat[i] = malloc(sizeof(*mat[i]) * n);
        for (int j = 0; j < n; j++) {
            mat[i][j] = res[i][j];
        }
    }
    return mat;
}

double **mat_subI(double **mat, int n) {
    for (int i = 0; i < n; i++) {
        mat[i][i] -= 1;
    }
    return mat;
}

void mat_mult(double **mat1, double **mat2, double **res, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            res[i][j] = 0;
            for (int k = 0; k < n; k++)
                res[i][j] += mat1[i][k] * mat2[k][j];
        }
    }
}

void print_matrix(double **mat, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%19.15f,", mat[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

int main() {
    int i, j;
    int dim = 3;
    int n = 3;
    double *a;
    double **wam_res, **ddg_res, **lnorm_res;
    double **vecs;
    vecs = (double**)calloc(n, sizeof(double*));
    assert(vecs != NULL);

    for (i = 0; i < n; i++) {
        a = (double*)calloc(dim, sizeof(double));
        assert(a != NULL);
        for (j = 0; j < dim; j++) {
            a[j] = rand() % 50;
        }
        vecs[i] = a;
    }
    wam_res = wam(vecs, dim, n);
    printf("%s", "wam:\n");
    print_matrix(wam_res, n);
    ddg_res = ddg(wam_res, n);
    printf("%s", "ddg:\n");
    print_sym_mat(ddg_res, n);

    lnorm_res = lnorm(ddg_res, wam_res, n);
    printf("%s", "lnorm:\n");
    print_matrix(lnorm_res, n);
    return 0;
}

double **wam(double **vecs, int dim, int n) {
    /* Calculating wam for the given set of vectors- n vectors of dimension dim */
    double **res;
    int i, j, curr_index = 0;
    double norm;
    res = (double **)calloc(n, sizeof(double*));
    for (i = 0; i < n; i++) {
        res[i] = (double*)calloc(n, sizeof(double));

        /*For vector i, calculating wij for all 0<=j<n*/
        for (j = 0; j < n; j++) {
            if (i != j) {
                norm = l2_norm2vec(*(vecs+i), *(vecs+j), dim); /* function for l2 norm */
                res[i][j] = exp((-1 * norm / 2));
                curr_index++;
            } else {
                /* no self loops */
                res[i][j] = 0;
            }
        }
    }
    return res;
}

double **ddg(double **wam, int n) {
    int i, j;
    double **res;
    double d = 0;
    res = (double**)calloc(n*n, sizeof(double*)); /* a n*n array with just 0s */
    for (i = 0; i < n; i++) {
        res[i] = (double*)calloc(n, sizeof(double));
        for (j = 0; j < n; j++)  {
            d += wam[i][j];
        }
        res[i][i] = 1 / sqrt(d);
        d = 0;
    }
    return res;
}

double **lnorm(double **ddg, double **wam, int n) {
    double **res, **tmp;
    int i;
    res = (double**)calloc(n, sizeof(double));
    tmp = (double**)calloc(n, sizeof(double));
    for (i = 0; i < n; i++) {
        res[i] = (double*)calloc(n, sizeof(double));
        tmp[i] = (double*)calloc(n, sizeof(double));
    }

    mat_mult(ddg,wam, res, n);
    printf("%s", "step 1 results: \n");
    print_sym_mat(res, n); /* errors happened here */
    mat_mult(copyMatrix(res, n), ddg, tmp, n);  /*copyMatrix- simple loop to allocate an identical matrix */
    res = copyMatrix(tmp, n);
    printf("%s", "step 2 results: \n");
    print_sym_mat(res, n); /* errors happened here too */
    mat_subI(res, n);
    printf("%s", "step 3 results: \n");
    print_sym_mat(res, n);
    return res;
}

Output with 8 decimal places:

wam:
  0.00000000,  0.00007291,  0.00000000,
  0.00007291,  0.00000000,  0.00000001,
  0.00000000,  0.00000001,  0.00000000,

ddg:
117.11258762,  0.00000000,  0.00000000,
  0.00000000,117.10123871,  0.00000000,
  0.00000000,  0.00000000,8087.59220061,

step 1 results:
  0.00000000,  0.00853872,  0.00000007,
  0.00853790,  0.00000000,  0.00000172,
  0.00000467,  0.00011897,  0.00000000,

step 2 results:
  0.00000000,  0.99989517,  0.00054713,
  0.99989517,  0.00000000,  0.01393205,
  0.00054713,  0.01393205,  0.00000000,

step 3 results:
 -1.00000000,  0.99989517,  0.00054713,
  0.99989517, -1.00000000,  0.01393205,
  0.00054713,  0.01393205, -1.00000000,

lnorm:
 -1.00000000,  0.99989517,  0.00054713,
  0.99989517, -1.00000000,  0.01393205,
  0.00054713,  0.01393205, -1.00000000,

With more places, it appears the coefficients in wam are quite different from 0.00000001:

wam:
  0.000000000000000,  0.000072910387337,  0.000000000577653,
  0.000072910387337,  0.000000000000000,  0.000000014710728,
  0.000000000577653,  0.000000014710728,  0.000000000000000,

ddg:
117.112587619035523,  0.000000000000000,  0.000000000000000,
  0.000000000000000,117.101238706028283,  0.000000000000000,
  0.000000000000000,  0.000000000000000,8087.592200608456551,

step 1 results:
  0.000000000000000,  0.008538724125301,  0.000000067650464,
  0.008537896671658,  0.000000000000000,  0.000001722644500,
  0.000004671823707,  0.000118974371006,  0.000000000000000,

step 2 results:
  0.000000000000000,  0.999895172041842,  0.000547129363250,
  0.999895172041842,  0.000000000000000,  0.013932046219111,
  0.000547129363250,  0.013932046219111,  0.000000000000000,

step 3 results:
 -1.000000000000000,  0.999895172041842,  0.000547129363250,
  0.999895172041842, -1.000000000000000,  0.013932046219111,
  0.000547129363250,  0.013932046219111, -1.000000000000000,

lnorm:
 -1.000000000000000,  0.999895172041842,  0.000547129363250,
  0.999895172041842, -1.000000000000000,  0.013932046219111,
  0.000547129363250,  0.013932046219111, -1.000000000000000,