Added triangular linear solvers

drreynolds · Jan 25, 2021 · 50de86c · 50de86c
1 parent 4015fa6
commit 50de86c
Show file tree

Hide file tree

Showing 6 changed files with 290 additions and 0 deletions.
diff --git a/TriangularSolves/BackwardSub.m b/TriangularSolves/BackwardSub.m
@@ -0,0 +1,45 @@
+function [x] = BackwardSub(U,y)
+  % usage: [x] = BackwardSub(U,y)
+  %
+  % Row-oriented backward substitution to solve the upper-triangular
+  % linear system
+  %       U x = y
+  % This function does not ensure that U is upper-triangular.  It does,
+  % however, attempt to catch the case where U is singular.
+  %
+  % Inputs:
+  %   U - square n-by-n matrix (assumed upper triangular)
+  %   y - right-hand side vector (n-by-1)
+  %
+  % Outputs:
+  %   x - solution vector (n-by-1)
+  %
+  % Daniel R. Reynolds
+  % SMU Mathematics
+  % Math 4315
+
+  % check inputs
+  [m,n] = size(U);
+  if (m ~= n)
+    error('BackwardSub error: matrix must be square')
+  end
+  [p,q] = size(y);
+  if ((p ~= n) || (q ~= 1))
+    error('BackwardSub error: right-hand side vector has incorrect dimensions')
+  end
+  if (min(abs(diag(U))) < 100*eps)
+    error('BackwardSub error: matrix is [close to] singular')
+  end
+
+  % create output vector
+  x = y;
+
+  % perform forward-subsitution algorithm
+  for i=n:-1:1
+    for j=i+1:n
+      x(i) = x(i) - U(i,j)*x(j);
+    end
+    x(i) = x(i)/U(i,i);
+  end
+
+  return
diff --git a/TriangularSolves/BackwardSub.py b/TriangularSolves/BackwardSub.py
@@ -0,0 +1,49 @@
+# BackwardSub.py
+#
+# Daniel R. Reynolds
+# SMU Mathematics
+# Math 4315
+
+# imports
+import numpy
+
+def BackwardSub(U,y):
+    """
+    usage: x = BackwardSub(U,y)
+
+    Row-oriented backward substitution to solve the upper-triangular
+    linear system
+          U x = y
+    This function does not ensure that U is upper-triangular.  It does,
+    however, attempt to catch the case where U is singular.
+
+    Inputs:
+      U - square n-by-n matrix (assumed upper triangular)
+      y - right-hand side vector (n-by-1)
+
+    Outputs:
+      x - solution vector (n-by-1)
+    """
+
+    # check inputs
+    m, n = numpy.shape(U)
+    if (m != n):
+        raise ValueError("BackwardSub error: matrix must be square")
+    p = numpy.size(y)
+    if (p != n):
+        raise ValueError("BackwardSub error: right-hand side vector has incorrect dimensions")
+    if (numpy.min(numpy.abs(numpy.diag(U))) < 100*numpy.finfo(float).eps):
+        raise ValueError("BackwardSub error: matrix is [close to] singular")
+
+    # create output vector
+    x = y.copy()
+
+    # perform forward-subsitution algorithm
+    for i in range(n-1,-1,-1):
+        for j in range(i+1,n):
+            x[i] -= U[i,j]*x[j]
+        x[i] /= U[i,i]
+
+    return x
+
+# end function
diff --git a/TriangularSolves/ForwardSub.m b/TriangularSolves/ForwardSub.m
@@ -0,0 +1,45 @@
+function [y] = ForwardSub(L,b)
+  % usage: [y] = ForwardSub(L,b)
+  %
+  % Row-oriented forward substitution to solve the lower-triangular
+  % linear system
+  %       L y = b
+  % This function does not ensure that L is lower-triangular.  It does,
+  % however, attempt to catch the case where L is singular.
+  %
+  % Inputs:
+  %   L - square n-by-n matrix (assumed lower triangular)
+  %   b - right-hand side vector (n-by-1)
+  %
+  % Outputs:
+  %   y - solution vector (n-by-1)
+  %
+  % Daniel R. Reynolds
+  % SMU Mathematics
+  % Math 4315
+
+  % check inputs
+  [m,n] = size(L);
+  if (m ~= n)
+    error('ForwardSub error: matrix must be square')
+  end
+  [p,q] = size(b);
+  if ((p ~= n) || (q ~= 1))
+    error('ForwardSub error: right-hand side vector has incorrect dimensions')
+  end
+  if (min(abs(diag(L))) < 100*eps)
+    error('ForwardSub error: matrix is [close to] singular')
+  end
+
+  % create output vector
+  y = b;
+
+  % perform forward-subsitution algorithm
+  for i=1:n
+    for j=1:i-1
+      y(i) = y(i) - L(i,j)*y(j);
+    end
+    y(i) = y(i)/L(i,i);
+  end
+
+  return
diff --git a/TriangularSolves/ForwardSub.py b/TriangularSolves/ForwardSub.py
@@ -0,0 +1,49 @@
+# ForwardSub.py
+#
+# Daniel R. Reynolds
+# SMU Mathematics
+# Math 4315
+
+# imports
+import numpy
+
+def ForwardSub(L,b):
+    """
+    usage: y = ForwardSub(L,b)
+
+    Row-oriented forward substitution to solve the lower-triangular
+    linear system
+          L y = b
+    This function does not ensure that L is lower-triangular.  It does,
+    however, attempt to catch the case where L is singular.
+
+    Inputs:
+      L - square n-by-n matrix (assumed lower triangular)
+      b - right-hand side vector (n-by-1)
+
+    Outputs:
+      y - solution vector (n-by-1)
+    """
+
+    # check inputs
+    m, n = numpy.shape(L)
+    if (m != n):
+        raise ValueError("ForwardSub error: matrix must be square")
+    p = numpy.size(b)
+    if (p != n):
+        raise ValueError("ForwardSub error: right-hand side vector has incorrect dimensions")
+    if (numpy.min(numpy.abs(numpy.diag(L))) < 100*numpy.finfo(float).eps):
+        raise ValueError("ForwardSub error: matrix is [close to] singular")
+
+    # create output vector
+    y = b.copy()
+
+    # perform forward-subsitution algorithm
+    for i in range(n):
+        for j in range(i):
+            y[i] -= L[i,j]*y[j]
+        y[i] /= L[i,i]
+
+    return y
+
+# end function
diff --git a/TriangularSolves/TestTriSolves.m b/TriangularSolves/TestTriSolves.m
@@ -0,0 +1,48 @@
+% Script to test ForwardSub and BackwardSub on a variety of linear systems.
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% Math 4315
+clear
+
+% set matrix sizes for tests
+nvals = [50, 100, 200, 400];
+
+% full-rank square matrix tests
+for n = nvals
+
+   fprintf('Testing with full-rank triangular matrices of dimension %i\n',n);
+
+   % create the matrices
+   L = tril(rand(n,n) + 2*eye(n));
+   U = triu(rand(n,n) + 2*eye(n));
+
+   % create solution vector
+   x = rand(n,1);
+
+   % solve triangular linear systems
+   x_bs = BackwardSub(U, U*x);
+   x_fs = ForwardSub(L, L*x);
+
+   % output results
+   fprintf('   BackwardSub error = %g\n', norm(x-x_bs));
+   fprintf('   ForwardSub error  = %g\n', norm(x-x_fs));
+
+end
+
+% ensure that rank-deficient case fails
+n = 100;
+fprintf('Testing with rank-deficient triangular matrices (should fail)\n');
+L = tril(rand(n,n));  L(n-3,n-3) = eps;
+U = triu(rand(n,n));  U(n-5,n-5) = eps;
+x = rand(n,1);
+try
+  x_bs = BackwardSub(U, U*x);
+catch
+  fprintf('   rank deficiency correctly caught by BackwardSub\n');
+end
+try
+  x_fs = ForwardSub(L, L*x);
+catch
+  fprintf('   rank deficiency correctly caught by ForwardSub\n');
+end
diff --git a/TriangularSolves/TestTriSolves.py b/TriangularSolves/TestTriSolves.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+#
+# Script to test ForwardSub and BackwardSub on a variety of linear systems.
+#
+# Daniel R. Reynolds
+# SMU Mathematics
+# Math 4315
+
+# imports
+import numpy
+from ForwardSub  import ForwardSub
+from BackwardSub import BackwardSub
+
+# set matrix sizes for tests
+nvals = [50, 100, 200, 400]
+
+# full-rank square matrix tests
+for n in nvals:
+
+    print("Testing with full-rank triangular matrices of dimension ", n)
+
+    # create the matrices
+    L = numpy.tril(numpy.random.rand(n,n) + 2*numpy.eye(n))
+    U = numpy.triu(numpy.random.rand(n,n) + 2*numpy.eye(n))
+
+    # create solution vector
+    x = numpy.random.rand(n)
+
+    # solve triangular linear systems
+    x_bs = BackwardSub(U, U@x)
+    x_fs = ForwardSub(L, L@x)
+
+    # output results
+    print("   BackwardSub error = ", numpy.linalg.norm(x-x_bs))
+    print("   ForwardSub error  = ", numpy.linalg.norm(x-x_fs))
+
+# ensure that rank-deficient case fails
+n = 100
+print("Testing with rank-deficient triangular matrices (should fail)")
+L = numpy.tril(numpy.random.rand(n,n))
+L[n-4,n-4] = numpy.finfo(float).eps
+U = numpy.triu(numpy.random.rand(n,n))
+U[n-6,n-6] = numpy.finfo(float).eps
+x = numpy.random.rand(n)
+try:
+  x_bs = BackwardSub(U, U@x)
+except:
+  print("   rank deficiency correctly caught by BackwardSub")
+try:
+  x_fs = ForwardSub(L, L@x)
+except:
+  print("   rank deficiency correctly caught by ForwardSub")
+
+# end of script