title: numerically stable inverse of a 2x2 matrix
body: <p>In a numerical solver I am working on in C, I need to invert a 2x2 matrix and it then gets multiplied on the right side by another matrix:</p>
  
  <pre><code>C = B . inv(A)
  </code></pre>
  
  <p>I have been using the following definition of an inverted 2x2 matrix:</p>
  
  <pre><code>a = A[0][0];
  b = A[0][1];
  c = A[1][0];
  d = A[1][1];
  invA[0][0] = d/(a*d-b*c);
  invA[0][1] = -b/(a*d-b*c);
  invA[1][0] = -c/(a*d-b*c);
  invA[1][1] = a/(a*d-b*c);
  </code></pre>
  
  <p>In the first few iterations of my solver this seems to give the correct answers, however, after a few steps things start to grow and eventually explode.</p>
  
  <p>Now, comparing to an implementation using SciPy, I found that the same math does not explode.  The only difference I can find is that the SciPy code uses <code>scipy.linalg.inv()</code>, which internally uses LAPACK internally to perform the inversion.</p>
  
  <p>When I replace the call to <code>inv()</code> with the above calculations the Python version does explode, so I'm pretty sure this is the problem.  Small differences in the calculations are creeping in, which leads me to believe it is a numerical problem--not entirely surprising for an inversion operation.</p>
  
  <p>I am using double-precision floats (64-bit), hoping that numerical issues would not be a problem, but apparently that is not the case.</p>
  
  <p>But: I would like to solve this in my C code without needing to call out to a library like LAPACK, because the whole reason for porting it to pure C is to get it running on a target system.  Moreover, I'd like to understand the problem, not just call out to a black box.  Eventually I'd like to it run with single-precision too, if possible.</p>
  
  <p>So, my question is, for such a small matrix, is there a numerically more stable way to calculate the inverse of A?</p>
  
  <p>Thanks.</p>
  
  <p>Edit:  Currently trying to figure out if I can just <a href="http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/" rel="nofollow">avoid the inversion</a> by solving for <code>C</code>.</p>