title: Symmetric Matrix Inversion in C using CBLAS/LAPACK
body: <p>I am writing an algorithm in C that requires Matrix and Vector multiplications. I have a matrix <strong><em>Q</em></strong> (W x W) which is created by multiplying the transpose of a vector <strong><em>J</em></strong>(1 x W) with itself and adding Identity matrix <strong><em>I</em></strong>, scaled using scalar <strong><em>a</em></strong>.  </p>
  
  <p>Q = [(J^T) * J + aI].</p>
  
  <p>I then have to multiply the <strong>inverse of Q</strong> with <strong>vector G</strong> to get vector <strong>M</strong>.</p>
  
  <p>M = (Q^(-1)) * G.</p>
  
  <p>I am using <em>cblas</em> and <em>clapack</em> to develop my algorithm. When matrix <strong>Q</strong> is populated using random numbers (type float) and inverted using the routines <em>sgetrf_</em> and <em>sgetri_</em> , the calculated inverse is <strong>correct</strong>.</p>
  
  <p><strong>But when matrix Q is symmetrical</strong>, which is the case when you multiply (J^T) x J, the calculated <strong>inverse is wrong!!</strong>.</p>
  
  <p>I am aware of the row-major (in C) and column-major (in FORTRAN) format of arrays while calling <em>lapack</em> routines from C, but for a symmetrical matrix this should not be a problem as A^T = A.</p>
  
  <p>I have attached my C function code for matrix inversion below.</p>
  
  <p>I am sure there is a better way to solve this. Can anyone help me with this?</p>
  
  <p>A solution using cblas would be great...</p>
  
  <p>Thanks.</p>
  
  <pre><code>void InverseMatrix_R(float *Matrix, int W)
  {   
      int     LDA = W;
      int     IPIV[W];
      int     ERR_INFO;
      int     LWORK = W * W;
      float   Workspace[LWORK];
  
      // - Compute the LU factorization of a M by N matrix A
      sgetrf_(&amp;W, &amp;W, Matrix, &amp;LDA, IPIV, &amp;ERR_INFO);
  
      // - Generate inverse of the matrix given its LU decompsotion
      sgetri_(&amp;W, Matrix, &amp;LDA, IPIV, Workspace, &amp;LWORK, &amp;ERR_INFO);
  
      // - Display the Inverted matrix
      PrintMatrix(Matrix, W, W);
  
  }
  
  void PrintMatrix(float* Matrix, int row, int colm)
  {
      int i,k;
  
      for (i =0; i &lt; row; i++) 
      {
          for (k = 0; k &lt; colm; k++) 
          {
              printf("%g, ",Matrix[i*colm + k]);
          }
  
          printf("\n");
      }
  }
  </code></pre>