05/30 - Fixed compilation, implemented tridiag dot product, tridiag s…

…olve needs debugging

:wq

@@ -0,0 +1,58 @@
+#include "../../include/tuvx/linalg/linalg.h"
+
+namespace tuvx {
+namespace linalg {
+
+template <typename T>
+std::vector<T> dot(tuvx::linalg::trid_mat<T> A, std::vector<T> x) {
+  int size = x.size();
+  std::vector<T> v(size);
+  v[0] = A.mdiag[0] * x[0] + A.udiag[0] * x[1];
+
+  int i = 0;
+  for (i = 1; i < x.size() - 1; i++) {
+    v[i] =
+        A.mdiag[i] * x[i] + A.udiag[i] * x[i + 1] + A.ldiag[i - 1] * x[i - 1];
+  }
+  v[i] = A.mdiag[i] * x[i] + A.ldiag[i - 1] * x[i - 1];
+  return v;
+}
+
+/* solve tridiagonal system */
+template <typename T>
+std::vector<T> tridiag_solve(trid_mat<T> A, std::vector<T> b) {
+  T temp = A.udiag[0] / A.mdiag[0];
+  std::vector<T> x(b.size());
+
+  // solve the first lowest diag
+  x[0] = b[0] / A.mdiag[0];
+
+  // forward iteration
+  for (int i = 1; i < b.size(); i++) {
+    temp = A.udiag[i - 1] / (A.mdiag[i] - A.ldiag[i - 1] * temp);
+    x[i] =
+        (b[i] - A.ldiag[i] * x[i - 1]) / (A.mdiag[i] - A.ldiag[i - 1] * temp);
+  }
+
+  // backward iteration
+  for (int i = b.size(); i > 0; i--) {
+    temp = (A.mdiag[i] - (A.ldiag[i - 1] / temp)) / A.ldiag[i - 1];
+    b[i] = b[i] - temp * b[i + 1];
+  }
+  return x;
+}
+/*
+    cp(1) = c(1) / b(1)
+    u(1) = r(1) / b(1)
+    do i = 2, size( b )
+      denom = 1.0 / ( b(i) - a(i) * cp(i-1) )
+      cp(i) = c(i) * denom
+      u(i) = ( r(i) - a(i) * u(i-1) ) * denom
+    end do
+    do i = size( b ) - 1, 1, -1
+      u(i) = u(i) - cp(i) * u(i+1)
+    end do
+
+ */
+} // namespace linalg
+} // namespace tuvx

include/tuvx/linalg/linalg.h

@@ -1,21 +1,22 @@
 #pragma once
 
-#include <vector>
 #include <iostream>
-namespace tuvx{
-namespace linalg{
+#include <vector>
+namespace tuvx {
+namespace linalg {
 
+// tridiag matrix
+template <typename T> struct trid_mat {
+  int size;
+  std::vector<T> udiag; // upper diagonal
+  std::vector<T> ldiag; // lower diagonal
+  std::vector<T> mdiag; // main diagonal
+};
 
-  // tridiag matrix
-  template <typename T>
-  struct trid_mat {
-    int size;
-    std::vector<T> udiag; // upper diagonal 
-    std::vector<T> ldiag; // lower diagonal
-    std::vector<T> mdiag; // main diagonal
-  };
+template <typename T>
+std::vector<T> tridiag_solve(trid_mat<T> A, std::vector<T> b);
 
-  template <typename T>
-  std::vector<T> tridiag_solve(trid_mat<T> A, std::vector<T> b); 
-} 
-}
+template <typename T>
+std::vector<T> dot(tuvx::linalg::trid_mat<T>, std::vector<T>);
+} // namespace linalg
+} // namespace tuvx

src/linalg/linalg.cpp

@@ -0,0 +1,58 @@
+#include "../../include/tuvx/linalg/linalg.h"
+
+namespace tuvx {
+namespace linalg {
+
+template <typename T>
+std::vector<T> dot(tuvx::linalg::trid_mat<T> A, std::vector<T> x) {
+  int size = x.size();
+  std::vector<T> v(size);
+  v[0] = A.mdiag[0] * x[0] + A.udiag[0] * x[1];
+
+  int i = 0;
+  for (i = 1; i < x.size() - 1; i++) {
+    v[i] =
+        A.mdiag[i] * x[i] + A.udiag[i] * x[i + 1] + A.ldiag[i - 1] * x[i - 1];
+  }
+  v[i] = A.mdiag[i] * x[i] + A.ldiag[i - 1] * x[i - 1];
+  return v;
+}
+
+/* solve tridiagonal system */
+template <typename T>
+std::vector<T> tridiag_solve(trid_mat<T> A, std::vector<T> b) {
+  T temp = A.udiag[0] / A.mdiag[0];
+  std::vector<T> x(b.size());
+
+  // solve the first lowest diag
+  x[0] = b[0] / A.mdiag[0];
+
+  // forward iteration
+  for (int i = 1; i < b.size(); i++) {
+    temp = A.udiag[i - 1] / (A.mdiag[i] - A.ldiag[i - 1] * temp);
+    x[i] =
+        (b[i] - A.ldiag[i] * x[i - 1]) / (A.mdiag[i] - A.ldiag[i - 1] * temp);
+  }
+
+  // backward iteration
+  for (int i = b.size(); i > 0; i--) {
+    b[i] = b[i] - temp * b[i + 1];
+    temp = (A.mdiag[i] - (A.ldiag[i - 1] / temp)) / A.ldiag[i - 1];
+  }
+  return x;
+}
+/*
+    cp(1) = c(1) / b(1)
+    u(1) = r(1) / b(1)
+    do i = 2, size( b )
+      denom = 1.0 / ( b(i) - a(i) * cp(i-1) )
+      cp(i) = c(i) * denom
+      u(i) = ( r(i) - a(i) * u(i-1) ) * denom
+    end do
+    do i = size( b ) - 1, 1, -1
+      u(i) = u(i) - cp(i) * u(i+1)
+    end do
+
+ */
+} // namespace linalg
+} // namespace tuvx

test/unit/linalg/test_tridiag.cpp

@@ -1,76 +1,74 @@
 #include "../../../include/tuvx/linalg/linalg.h"
-#include "../../../src/linalg/tridiag.cpp"
-
+#include "../../../src/linalg/linalg.cpp"
 #include <random>
 
-typedef tuvx::linalg::trid_mat<double> trid_matd; 
+typedef tuvx::linalg::trid_mat<double> trid_matd;
 typedef std::vector<double> vecd;
 
-template <typename T>
-void fill_rand(std::vector<T> &x, int size);
-
-template <typename T>
-void fill_mat(tuvx::linalg::trid_mat<T> &trid_mat, int n); 
-
-template <typename T>
-std::vector<T> dot(tuvx::linalg::trid_mat<T>, std::vector<T>);
-
-template <typename T>
-void print_vec(std::vector<T>);
+template <typename T> void fill_rand(std::vector<T> &x, int size);
+template <typename T> void fill_mat(tuvx::linalg::trid_mat<T> &trid_mat, int n);
+template <typename T> void print_vec(std::vector<T>);
+template <typename T> void print_trid_mat(tuvx::linalg::trid_mat<T> x);
 
 int main() {
-
-    trid_matd M;
-    vecd b; 
-    vecd x; 
-    vecd x_sol; 
-
-    fill_mat(M, 5);
-    fill_rand(b, 5);
-    fill_rand(x, 5);
-
-    x_sol = tuvx::linalg::tridiag_solve<double>(M, b);
-
-    return 0; 
+  trid_matd M;
+  vecd b;
+  vecd x = {1, 6, 2, 3, 5};
+  M.udiag = {3, 1, 5, 3};
+  M.mdiag = {3, 2, 6, 2, 4};
+  M.ldiag = {1, 3, 6, 5};
+  M.size = 5;
+  vecd x_sol;
+
+  print_trid_mat(M);
+  print_vec(x);
+  b = dot(M, x);
+  print_vec(tuvx::linalg::tridiag_solve(M, b));
+
+  /*
+  fill_mat(M, 5);
+  fill_rand(b, 5);
+  fill_rand(x, 5);
+
+  print_trid_mat(M);
+  print_vec(x);
+  x = tuvx::linalg::dot(M, x);
+  x_sol = tuvx::linalg::tridiag_solve<double>(M, b);
+
+  print_vec(x);
+  */
+
+  return 0;
 }
 
-template <typename T>
-void fill_rand(std::vector<T> &x, int size){    
-    std::random_device dev;
-    std::mt19937 rng(dev());
-    std::uniform_int_distribution<std::mt19937::result_type> dist6(1,6);
-
-    x = std::vector<T>(size);
-    for(int i = 0; i < x.size(); i++) {
-        x[i] = (T)dist6(rng);
-    }
+template <typename T> void fill_rand(std::vector<T> &x, int size) {
+  std::random_device dev;
+  std::mt19937 rng(dev());
+  std::uniform_int_distribution<std::mt19937::result_type> dist6(1, 6);
+  x = std::vector<T>(size);
+  for (int i = 0; i < x.size(); i++) {
+    x[i] = (T)dist6(rng);
+  }
 }
 
-
-template <typename T>
-void fill_mat(tuvx::linalg::trid_mat<T> &trid_mat, int n) {
-    fill_rand(trid_mat.mdiag, n);
-    fill_rand(trid_mat.ldiag, n);
-    fill_rand(trid_mat.udiag, n);
+template <typename T> void fill_mat(tuvx::linalg::trid_mat<T> &tridmat, int n) {
+  tridmat.size = n;
+  fill_rand(tridmat.mdiag, n);
+  fill_rand(tridmat.ldiag, n - 1);
+  fill_rand(tridmat.udiag, n - 1);
 }
 
-template <typename T>
-void print_vec(std::vector<T> x){
-    for (int i = 0; i < (int)x.size(); i++)
-        std::cout<< x.at(i) <<std::endl;
+template <typename T> void print_vec(std::vector<T> x) {
+  std::cout << std::endl;
+  for (int i = 0; i < (int)x.size(); i++) {
+    std::cout << x.at(i) << std::endl;
+  }
+  std::cout << std::endl;
 }
 
-template <typename T>
-std::vector<T> dot(tuvx::linalg::trid_mat<T> A, std::vector<T> x){
-    std::vector<T> v(x.size());
-    v[0] = A.mdiag[0] * x[0] + A.udiag[0]*x[1];
-
-    for (int i = 1; i < x.size(); i++) {
-        v[i] = A.mdiag[i]*x[i]   + 
-               A.udiag[i-1]*x[i+1] + 
-               A.ldiag[i-1]*x[i-1];
-    }
-
-    v[i] = A.mdiag[i]*x[i] + A.ldiag[i-1]* 
-    return v;
-}
+template <typename T> void print_trid_mat(tuvx::linalg::trid_mat<T> x) {
+  print_vec(x.udiag);
+  print_vec(x.mdiag);
+  print_vec(x.ldiag);
+  std::cout << "----" << std::endl;
+}

�LW

@@ -0,0 +1,58 @@
+#include "../../include/tuvx/linalg/linalg.h"
+
+namespace tuvx {
+namespace linalg {
+
+template <typename T>
+std::vector<T> dot(tuvx::linalg::trid_mat<T> A, std::vector<T> x) {
+  int size = x.size();
+  std::vector<T> v(size);
+  v[0] = A.mdiag[0] * x[0] + A.udiag[0] * x[1];
+
+  int i = 0;
+  for (i = 1; i < x.size() - 1; i++) {
+    v[i] =
+        A.mdiag[i] * x[i] + A.udiag[i] * x[i + 1] + A.ldiag[i - 1] * x[i - 1];
+  }
+  v[i] = A.mdiag[i] * x[i] + A.ldiag[i - 1] * x[i - 1];
+  return v;
+}
+
+/* solve tridiagonal system */
+template <typename T>
+std::vector<T> tridiag_solve(trid_mat<T> A, std::vector<T> b) {
+  T temp = A.udiag[0] / A.mdiag[0];
+  std::vector<T> x(b.size());
+
+  // solve the first lowest diag
+  x[0] = b[0] / A.mdiag[0];
+
+  // forward iteration
+  for (int i = 1; i < b.size(); i++) {
+    temp = A.udiag[i - 1] / (A.mdiag[i] - A.ldiag[i - 1] * temp);
+    x[i] =
+        (b[i] - A.ldiag[i] * x[i - 1]) / (A.mdiag[i] - A.ldiag[i - 1] * temp);
+  }
+
+  // backward iteration
+  for (int i = b.size(); i > 0; i--) {
+    temp = (A.mdiag[i] - (A.ldiag[i - 1] / temp)) / A.ldiag[i - 1];
+    b[i] = b[i] - temp * b[i + 1];
+  }
+  return x;
+}
+/*
+    cp(1) = c(1) / b(1)
+    u(1) = r(1) / b(1)
+    do i = 2, size( b )
+      denom = 1.0 / ( b(i) - a(i) * cp(i-1) )
+      cp(i) = c(i) * denom
+      u(i) = ( r(i) - a(i) * u(i-1) ) * denom
+    end do
+    do i = size( b ) - 1, 1, -1
+      u(i) = u(i) - cp(i) * u(i+1)
+    end do
+
+ */
+} // namespace linalg
+} // namespace tuvx