Add pointwise matrix inversion.

1 files changed, 21 insertions, 0 deletions

diff --git a/utils.py b/utils.py
new file mode 100644
index 0000000..b225bbe
--- /dev/null
+++ b/utils.py
@@ -0,0 +1,21 @@
+import numpy as np
+
+def matrix_invert(mat):
+    """
+    Pointwise matrix inversion.
+
+    Given an array that stores many square matrices (e.g. the values of the
+    metric tensor at spacetime points), compute the array of corresponding
+    inverses.
+
+    :param array_like mat: N-D array with N>2 and each mat[i0, ..., :, :] a
+                           square matrix to be inverted.
+    :return: Array (same shape as mat) of inverses.
+    """
+    oldshape = mat.shape
+    newshape = oldshape[:2] + (np.product(mat.shape[2:]),)
+
+    mat = mat.reshape(newshape).transpose((2, 0, 1))
+    inv = np.linalg.inv(mat)
+    return inv.transpose((1, 2, 0)).reshape(oldshape)
+