[BF] Fix bug in multiple_fast_inv.

VirgileFritsch · VirgileFritsch · commit d915f892eb06 · 2013-01-18T14:57:59.000+01:00
Some computations were done outside the loop, but they depend
on values changing at each iteration.
diff --git a/nipy/algorithms/statistics/utils.py b/nipy/algorithms/statistics/utils.py
@@ -43,32 +43,31 @@ def multiple_fast_inv(a):
     from scipy.linalg import calc_lwork
     from scipy.linalg.lapack import  get_lapack_funcs
     a1, n = a[0], a.shape[0]
-    getrf, getri = get_lapack_funcs(('getrf','getri'), (a1,))
-    if getrf.module_name[:7] == 'clapack' != getri.module_name[:7]:
-        # ATLAS 3.2.1 has getrf but not getri.
-        for i in range(n):
+    getrf, getri = get_lapack_funcs(('getrf', 'getri'), (a1,))
+    for i in range(n):
+        if getrf.module_name[:7] == 'clapack' != getri.module_name[:7]:
+            # ATLAS 3.2.1 has getrf but not getri.
             lu, piv, info = getrf(np.transpose(a[i]), rowmajor=0,
                                   overwrite_a=True)
             a[i] = np.transpose(lu)
-    else:
-        for i in range(n):
+        else:
             a[i], piv, info = getrf(a[i], overwrite_a=True)
-    if info == 0:
-        if getri.module_name[:7] == 'flapack':
-            lwork = calc_lwork.getri(getri.prefix, a1.shape[0])
-            lwork = lwork[1]
-            # XXX: the following line fixes curious SEGFAULT when
-            # benchmarking 500x500 matrix inverse. This seems to
-            # be a bug in LAPACK ?getri routine because if lwork is
-            # minimal (when using lwork[0] instead of lwork[1]) then
-            # all tests pass. Further investigation is required if
-            # more such SEGFAULTs occur.
-            lwork = int(1.01 * lwork)
-            for i in range(n):
+        if info == 0:
+            if getri.module_name[:7] == 'flapack':
+                lwork = calc_lwork.getri(getri.prefix, a1.shape[0])
+                lwork = lwork[1]
+                # XXX: the following line fixes curious SEGFAULT when
+                # benchmarking 500x500 matrix inverse. This seems to
+                # be a bug in LAPACK ?getri routine because if lwork is
+                # minimal (when using lwork[0] instead of lwork[1]) then
+                # all tests pass. Further investigation is required if
+                # more such SEGFAULTs occur.
+                lwork = int(1.01 * lwork)
                 a[i], _ = getri(a[i], piv, lwork=lwork, overwrite_lu=1)
-        else: # clapack
-            for i in range(n):
+            else:  # clapack
                 a[i], _ = getri(a[i], piv, overwrite_lu=1)
+        else:
+            raise ValueError('Matrix LU decomposition failed')
     return a
 
 
