Add triangular solve

3 files changed, 68 insertions, 25 deletions

diff --git a/examples/linear_system.cpp b/examples/linear_system.cpp
new file mode 100644
index 0000000..184f576
--- /dev/null
+++ b/examples/linear_system.cpp
@@ -0,0 +1,36 @@
+#include "matrix.hpp"
+#include "norms.hpp"
+#include "triangular_solve.hpp"
+#include "vector.hpp"
+
+#include <iomanip>
+#include <iostream>
+
+int main() {
+    const linalg::Matrix basis = linalg::Matrix::identity(3);
+    const linalg::Matrix upper{
+        {4.0, -2.0, 1.0},
+        {0.0, 3.0, 5.0},
+        {0.0, 0.0, -2.0}
+    };
+    const linalg::Vector expected{2.0, -1.0, 3.0};
+    const linalg::Vector rhs = upper * expected;
+    const linalg::Vector x = linalg::backward_substitution(upper, rhs);
+
+    std::cout << "Week 3 triangular solve demo\n";
+    std::cout << "Identity matrix diagonal: ";
+    for (std::size_t i = 0; i < basis.rows(); ++i) {
+        std::cout << basis(i, i) << (i + 1 == basis.rows() ? '\n' : ' ');
+    }
+
+    std::cout << "Recovered solution x: ";
+    for (std::size_t i = 0; i < x.size(); ++i) {
+        std::cout << std::fixed << std::setprecision(2) << x[i]
+                  << (i + 1 == x.size() ? '\n' : ' ');
+    }
+
+    const linalg::Vector residual = (upper * x) - rhs;
+    std::cout << "Residual 2-norm = " << linalg::norm2(residual) << '\n';
+
+    return 0;
+}
diff --git a/examples/matmul.cpp b/examples/matmul.cpp
new file mode 100644
index 0000000..c6245f3
--- /dev/null
+++ b/examples/matmul.cpp
@@ -0,0 +1,32 @@
+#include "matrix.hpp"
+
+#include <iomanip>
+#include <iostream>
+
+int main() {
+    const linalg::Matrix a{
+        {1.0, 2.0, 3.0},
+        {4.0, 5.0, 6.0}
+    };
+    const linalg::Matrix b{
+        {7.0, 8.0},
+        {9.0, 10.0},
+        {11.0, 12.0}
+    };
+
+    const linalg::Matrix c = a * b;
+
+    std::cout << "Matrix multiplication\n";
+    std::cout << "A is " << a.rows() << " x " << a.cols() << '\n';
+    std::cout << "B is " << b.rows() << " x " << b.cols() << '\n';
+    std::cout << "C = A * B:\n";
+
+    for (std::size_t i = 0; i < c.rows(); ++i) {
+        for (std::size_t j = 0; j < c.cols(); ++j) {
+            std::cout << std::fixed << std::setprecision(2) << c(i, j)
+                      << (j + 1 == c.cols() ? '\n' : ' ');
+        }
+    }
+
+    return 0;
+}
diff --git a/examples/solve_linear_system.cpp b/examples/solve_linear_system.cpp
deleted file mode 100644
index af45094..0000000
--- a/examples/solve_linear_system.cpp
+++ /dev/null
@@ -1,25 +0,0 @@
-#include "matrix.hpp"
-#include "vector.hpp"
-
-#include <iomanip>
-#include <iostream>
-
-int main() {
-    const linalg::Matrix basis = linalg::Matrix::identity(3);
-    const linalg::Vector x{1.0, -2.0, 0.5};
-
-    std::cout << "Identity matrix diagonal: ";
-    for (std::size_t i = 0; i < basis.rows(); ++i) {
-        std::cout << basis(i, i) << (i + 1 == basis.rows() ? '\n' : ' ');
-    }
-
-    std::cout << "Vector x: ";
-    for (std::size_t i = 0; i < x.size(); ++i) {
-        std::cout << std::fixed << std::setprecision(2) << x[i]
-                  << (i + 1 == x.size() ? '\n' : ' ');
-    }
-
-    std::cout << "x dot x = " << linalg::dot(x, x) << '\n';
-
-    return 0;
-}