title: Inaccuracies in matrix multiplication- C
body: <p>I'm doing a DS project in C.
  I've been struggling with a matrix multiplication function:</p>
  <pre><code>void mat_mult(double **mat1, double **mat2, double **res, int n) {
      int i, j, k;
      
      for (i = 0; i &lt; n; i++) {
          
          for (j = 0; j &lt; n; j++) {
              assert(res[i][j] == 0);
              for (k = 0; k &lt; n; k++)
                  res[i][j] += mat1[i][k] * mat2[k][j];;
          }
      }
  } 
  </code></pre>
  <p>The matrices I'm trying to multiply are:</p>
  <pre><code>0.00000000, 0.00007291, 0.00000000 
  0.00007291, 0.00000000, 0.00000001 
  0.00000000, 0.00000001, 0.00000000
  </code></pre>
  <p>and:</p>
  <pre><code>117.11258762, 0.00000000, 0.00000000 
  0.00000000, 117.10123871, 0.00000000 
  0.00000000, 0.00000000, 8087.59220061
  </code></pre>
  <p>The result is:</p>
  <pre><code>0.00000000, 0.00853872, 0.00000007 
  0.00853790, 0.00000000, 0.00000172 
  0.00000467, 0.00011897, 0.00000000 
  </code></pre>
  <p>And not:</p>
  <pre><code>0.00000000 0.00853868 0.00000000
  0.00853785 0.00000000 0.00000117
  0.00000000 0.00008088 0.00000000
  </code></pre>
  <p>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?</p>
  <p>***Edit- generating code ***</p>
  <pre><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&lt;n; i++)
  {
      a= (double*)calloc(dim,sizeof(double));
      assert(a!=NULL);
      for(j=0; j&lt;dim; j++)
      {
          a[j]= rand() %50;
      }
      vecs[i]= a;
      
  }
  wam_res= wam(vecs,dim,n);
  printf(&quot;%s&quot;,&quot;wam:\n&quot;);
  print_matrix(wam_res, n);
  ddg_res= ddg(wam_res, n);
  printf(&quot;%s&quot;,&quot;ddg:\n&quot;);
  print_sym_mat(ddg_res, n); 
  
  lnorm_res= lnorm(ddg_res, wam_res, n);
  printf(&quot;%s&quot;,&quot;lnorm:\n&quot;); }
  
  
  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&lt; n; i++)
  {
      res[i]= (double*)calloc(n, sizeof(double));
  
      /*For vector i, calculating wij for all 0&lt;=j&lt;n*/
      for(j=0; j&lt;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&lt;n; i++)
  {
      res[i]= (double*)calloc(n, sizeof(double));
      for(j=0; j&lt;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&lt;n; i++)
  {
      res[i]= (double*)calloc(n, sizeof(double));
      tmp[i]= (double*)calloc(n, sizeof(double));
  }
  
  mat_mult(ddg,wam, res, n);
  printf(&quot;%s&quot;, &quot;step 1 results: \n&quot;);
  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(&quot;%s&quot;, &quot;step 2 results: \n&quot;);
  print_sym_mat(res,n); /* errors happened here too */
  mat_subI(res,n); 
  printf(&quot;%s&quot;, &quot;step 3 results: \n&quot;);
  print_sym_mat(res,n);
  return res; }
  </code></pre>
  <p>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</p>