unify svd computations

FIAT/alfeld_sorokina.py

@@ -41,10 +41,8 @@ def AlfeldSorokinaSpace(ref_el, degree):
 
     if len(rows) > 0:
         dual_mat = numpy.vstack(rows)
-        _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True)
-        tol = sig[0] * 1E-10
-        num_sv = len([s for s in sig if abs(s) > tol])
-        coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(-1, 0))
+        nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True)
+        coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0))
     return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs)
 
 

FIAT/christiansen_hu.py

@@ -28,10 +28,8 @@ def ChristiansenHuSpace(ref_el, degree, reduced=False):
     tab = C0.tabulate(Q.get_points(), 1)
     divC0 = sum(tab[alpha][:, alpha.index(1), :] for alpha in tab if sum(alpha) == 1)
 
-    _, sig, vt = numpy.linalg.svd(divC0.T, full_matrices=True)
-    tol = sig[0] * 1E-10
-    num_sv = len([s for s in sig if abs(s) > tol])
-    coeffs = numpy.tensordot(vt[num_sv:], C0.get_coeffs(), axes=(-1, 0))
+    nsp = polynomial_set.spanning_basis(divC0.T, nullspace=True)
+    coeffs = numpy.tensordot(nsp, C0.get_coeffs(), axes=(-1, 0))
 
     verts = numpy.array(ref_complex.get_vertices())
     WT = verts[-1]

FIAT/macro.py

@@ -429,10 +429,7 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs):
 
         if len(rows) > 0:
             dual_mat = numpy.vstack(rows)
-            _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True)
-            tol = sig[0] * 1E-10
-            num_sv = len([s for s in sig if abs(s) > tol])
-            coeffs = vt[num_sv:]
+            coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True)
         else:
             coeffs = numpy.eye(expansion_set.get_num_members(degree))
 
@@ -445,10 +442,8 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs):
             jumps = expansion_set.tabulate_jumps(degree, points, order=vorder)
             for r in range(sorder+1, vorder+1):
                 dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T)
-                _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True)
-                tol = sig[0] * 1E-10
-                num_sv = len([s for s in sig if abs(s) > tol])
-                coeffs = numpy.dot(vt[num_sv:], coeffs)
+                nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True)
+                coeffs = numpy.dot(nsp, coeffs)
 
         if shape != ():
             m, n = coeffs.shape
@@ -499,8 +494,7 @@ def __init__(self, ref_el, degree, order=0, **kwargs):
 
         if len(rows) > 0:
             dual_mat = numpy.vstack(rows)
-            _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True)
-            num_sv = len([s for s in sig if abs(s) > 1.e-10])
-            coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(1, 0))
+            nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True)
+            coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0))
 
         super().__init__(ref_el, degree, degree, expansion_set, coeffs)

FIAT/polynomial_set.py

@@ -162,28 +162,25 @@ def form_matrix_product(mats, alpha):
     return result
 
 
+def spanning_basis(A, nullspace=False, rtol=1e-10):
+    """Construct a basis that spans the rows of A via SVD.
+    """
+    Aflat = A.reshape(A.shape[0], -1)
+    u, sig, vt = numpy.linalg.svd(Aflat, full_matrices=True)
+    atol = rtol * (sig[0] + 1)
+    num_sv = len([s for s in sig if abs(s) > atol])
+    basis = vt[num_sv:] if nullspace else vt[:num_sv]
+    return numpy.reshape(basis, (-1, *A.shape[1:]))
+
+
 def polynomial_set_union_normalized(A, B):
     """Given polynomial sets A and B, constructs a new polynomial set
     whose span is the same as that of span(A) union span(B).  It may
     not contain any of the same members of the set, as we construct a
     span via SVD.
     """
-    new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs))
-    func_shape = new_coeffs.shape[1:]
-    if len(func_shape) == 1:
-        (u, sig, vt) = numpy.linalg.svd(new_coeffs)
-        num_sv = len([s for s in sig if abs(s) > 1.e-10])
-        coeffs = vt[:num_sv]
-    else:
-        new_shape0 = new_coeffs.shape[0]
-        new_shape1 = numpy.prod(func_shape)
-        newshape = (new_shape0, new_shape1)
-        nc = numpy.reshape(new_coeffs, newshape)
-        (u, sig, vt) = numpy.linalg.svd(nc, 1)
-        num_sv = len([s for s in sig if abs(s) > 1.e-10])
-
-        coeffs = numpy.reshape(vt[:num_sv], (num_sv,) + func_shape)
-
+    new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0)
+    coeffs = spanning_basis(new_coeffs)
     return PolynomialSet(A.get_reference_element(),
                          A.get_degree(),
                          A.get_embedded_degree(),