rotate A by a unitary (real) rotation matrix R: RAR = transpose(R) * A * R */
| 15 | RAR = transpose(R) * A * R |
| 16 | */ |
| 17 | void sym_matrix_rotate(symmetric_matrix *RAR, const symmetric_matrix *A_, const double R[3][3]) { |
| 18 | int i, j; |
| 19 | double A[3][3], AR[3][3]; |
| 20 | A[0][0] = A_->m00; |
| 21 | A[1][1] = A_->m11; |
| 22 | A[2][2] = A_->m22; |
| 23 | A[0][1] = A[1][0] = A_->m01; |
| 24 | A[0][2] = A[2][0] = A_->m02; |
| 25 | A[1][2] = A[2][1] = A_->m12; |
| 26 | for (i = 0; i < 3; ++i) |
| 27 | for (j = 0; j < 3; ++j) |
| 28 | AR[i][j] = A[i][0] * R[0][j] + A[i][1] * R[1][j] + A[i][2] * R[2][j]; |
| 29 | for (i = 0; i < 3; ++i) |
| 30 | for (j = i; j < 3; ++j) |
| 31 | A[i][j] = R[0][i] * AR[0][j] + R[1][i] * AR[1][j] + R[2][i] * AR[2][j]; |
| 32 | RAR->m00 = A[0][0]; |
| 33 | RAR->m11 = A[1][1]; |
| 34 | RAR->m22 = A[2][2]; |
| 35 | RAR->m01 = A[0][1]; |
| 36 | RAR->m02 = A[0][2]; |
| 37 | RAR->m12 = A[1][2]; |
| 38 | } |
| 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) { |