Set Vinv = inverse of V, where both V and Vinv are real-symmetric matrices.*/
| 39 | |
| 40 | /* Set Vinv = inverse of V, where both V and Vinv are real-symmetric matrices.*/ |
| 41 | void sym_matrix_invert(symmetric_matrix *Vinv, const symmetric_matrix *V) { |
| 42 | double m00 = V->m00, m11 = V->m11, m22 = V->m22; |
| 43 | double m01 = V->m01, m02 = V->m02, m12 = V->m12; |
| 44 | |
| 45 | if (m01 == 0.0 && m02 == 0.0 && m12 == 0.0) { |
| 46 | /* optimize common case of a diagonal matrix: */ |
| 47 | Vinv->m00 = 1.0 / m00; |
| 48 | Vinv->m11 = 1.0 / m11; |
| 49 | Vinv->m22 = 1.0 / m22; |
| 50 | Vinv->m01 = Vinv->m02 = Vinv->m12 = 0.0; |
| 51 | } |
| 52 | else { |
| 53 | double detinv; |
| 54 | |
| 55 | /* compute the determinant: */ |
| 56 | detinv = m00 * m11 * m22 - m02 * m11 * m02 + 2.0 * m01 * m12 * m02 - m01 * m01 * m22 - |
| 57 | m12 * m12 * m00; |
| 58 | |
| 59 | if (detinv == 0.0) meep::abort("singular 3x3 matrix"); |
| 60 | |
| 61 | detinv = 1.0 / detinv; |
| 62 | |
| 63 | Vinv->m00 = detinv * (m11 * m22 - m12 * m12); |
| 64 | Vinv->m11 = detinv * (m00 * m22 - m02 * m02); |
| 65 | Vinv->m22 = detinv * (m11 * m00 - m01 * m01); |
| 66 | |
| 67 | Vinv->m02 = detinv * (m01 * m12 - m11 * m02); |
| 68 | Vinv->m01 = detinv * (m12 * m02 - m01 * m22); |
| 69 | Vinv->m12 = detinv * (m01 * m02 - m00 * m12); |
| 70 | } |
| 71 | } |
| 72 | |
| 73 | /* Returns whether or not V is positive-definite. */ |
| 74 | int sym_matrix_positive_definite(symmetric_matrix *V) { |
no test coverage detected