From f264b75ce6a22b505bb0c1b2d017c4402612ce9d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= <developer@crozet.re>
Date: Fri, 31 May 2013 11:28:42 +0200
Subject: [PATCH] Add dynamically sized matrix.

---
 src/nalgebra.rc  |   2 +-
 src/ndim/dmat.rs | 244 +++++++++++++++++++++++++++++++++++++++++++++++
 src/ndim/dvec.rs |   8 +-
 src/ndim/nmat.rs | 136 +++-----------------------
 src/ndim/nvec.rs |   6 +-
 5 files changed, 265 insertions(+), 131 deletions(-)
 create mode 100644 src/ndim/dmat.rs

diff --git a/src/nalgebra.rc b/src/nalgebra.rc
index 42111bed..348546e0 100644
--- a/src/nalgebra.rc
+++ b/src/nalgebra.rc
@@ -39,7 +39,7 @@ pub mod dim3
 /// n-dimensional linear algebra (slower).
 pub mod ndim
 {
-  // pub mod dmat;
+  pub mod dmat;
   pub mod dvec;
   pub mod nvec;
   pub mod nmat;
diff --git a/src/ndim/dmat.rs b/src/ndim/dmat.rs
new file mode 100644
index 00000000..aee81c94
--- /dev/null
+++ b/src/ndim/dmat.rs
@@ -0,0 +1,244 @@
+use core::num::{One, Zero};
+use core::vec::{from_elem, swap, all, all2, len};
+use core::cmp::ApproxEq;
+use traits::inv::Inv;
+use traits::transpose::Transpose;
+use traits::workarounds::rlmul::{RMul, LMul};
+use ndim::dvec::{DVec, zero_vec_with_dim};
+
+#[deriving(Eq, ToStr, Clone)]
+pub struct DMat<T>
+{
+  dim: uint, // FIXME: handle more than just square matrices
+  mij: ~[T]
+}
+
+pub fn zero_mat_with_dim<T: Zero + Copy>(dim: uint) -> DMat<T>
+{ DMat { dim: dim, mij: from_elem(dim * dim, Zero::zero()) } }
+
+pub fn is_zero_mat<T: Zero>(mat: &DMat<T>) -> bool
+{ all(mat.mij, |e| e.is_zero()) }
+
+pub fn one_mat_with_dim<T: Copy + One + Zero>(dim: uint) -> DMat<T>
+{
+  let mut res = zero_mat_with_dim(dim);
+  let     _1  = One::one::<T>();
+
+  for uint::range(0u, dim) |i|
+  { res.set(i, i, &_1); }
+
+  res
+}
+
+impl<T: Copy> DMat<T>
+{
+  pub fn offset(&self, i: uint, j: uint) -> uint
+  { i * self.dim + j }
+
+  pub fn set(&mut self, i: uint, j: uint, t: &T)
+  {
+    assert!(i < self.dim);
+    assert!(j < self.dim);
+    self.mij[self.offset(i, j)] = *t
+  }
+}
+
+impl<T: Copy> Index<(uint, uint), T> for DMat<T>
+{
+  fn index(&self, &(i, j): &(uint, uint)) -> T
+  { self.mij[self.offset(i, j)] }
+}
+
+impl<T: Copy + Mul<T, T> + Add<T, T> + Zero>
+Mul<DMat<T>, DMat<T>> for DMat<T>
+{
+  fn mul(&self, other: &DMat<T>) -> DMat<T>
+  {
+    assert!(self.dim == other.dim);
+
+    let     dim = self.dim;
+    let mut res = zero_mat_with_dim(dim);
+
+    for uint::range(0u, dim) |i|
+    {
+      for uint::range(0u, dim) |j|
+      {
+        let mut acc: T = Zero::zero();
+
+        for uint::range(0u, dim) |k|
+        { acc += self[(i, k)] * other[(k, j)]; }
+
+        res.set(i, j, &acc);
+      }
+    }
+
+    res
+  }
+}
+
+impl<T: Copy + Add<T, T> + Mul<T, T> + Zero>
+RMul<DVec<T>> for DMat<T>
+{
+  fn rmul(&self, other: &DVec<T>) -> DVec<T>
+  {
+    assert!(self.dim == len(other.at));
+
+    let     dim           = self.dim;
+    let mut res : DVec<T> = zero_vec_with_dim(dim);
+
+    for uint::range(0u, dim) |i|
+    {
+      for uint::range(0u, dim) |j|
+      { res.at[i] = res.at[i] + other.at[j] * self[(i, j)]; }
+    }
+
+    res
+  }
+}
+
+impl<T: Copy + Add<T, T> + Mul<T, T> + Zero>
+LMul<DVec<T>> for DMat<T>
+{
+  fn lmul(&self, other: &DVec<T>) -> DVec<T>
+  {
+    assert!(self.dim == len(other.at));
+
+    let     dim           = self.dim;
+    let mut res : DVec<T> = zero_vec_with_dim(dim);
+
+    for uint::range(0u, dim) |i|
+    {
+      for uint::range(0u, dim) |j|
+      { res.at[i] = res.at[i] + other.at[j] * self[(j, i)]; }
+    }
+
+    res
+  }
+}
+
+impl<T: Clone + Copy + Eq + One + Zero +
+        Mul<T, T> + Div<T, T> + Sub<T, T> + Neg<T>>
+Inv for DMat<T>
+{
+  fn inverse(&self) -> DMat<T>
+  {
+    let mut res : DMat<T> = self.clone();
+
+    res.invert();
+
+    res
+  }
+
+  fn invert(&mut self)
+  {
+    let     dim = self.dim;
+    let mut res = one_mat_with_dim::<T>(dim);
+    let     _0T = Zero::zero::<T>();
+
+    // inversion using Gauss-Jordan elimination
+    for uint::range(0u, dim) |k|
+    {
+      // search a non-zero value on the k-th column
+      // FIXME: would it be worth it to spend some more time searching for the
+      // max instead?
+
+      // FIXME: this is kind of uggly…
+      // … but we cannot use position_between since we are iterating on one
+      // columns
+      let mut n0 = 0u; // index of a non-zero entry
+
+      while (n0 != dim)
+      {
+        if (self[(n0, k)] != _0T)
+        { break; }
+
+        n0 += 1;
+      }
+
+      assert!(n0 != dim); // non inversible matrix
+
+      // swap pivot line
+      if (n0 != k)
+      {
+        for uint::range(0u, dim) |j|
+        {
+          let off_n0_j = self.offset(n0, j);
+          let off_k_j  = self.offset(k, j);
+
+          swap(self.mij, off_n0_j, off_k_j);
+          swap(res.mij,  off_n0_j, off_k_j);
+        }
+      }
+
+      let pivot = self[(k, k)];
+
+      for uint::range(k, dim) |j|
+      {
+        // FIXME: not to putting selfal exression directly on the nuction call
+        // is uggly but does not seem to compile any more…
+        let selfval = &(self[(k, j)] / pivot);
+        let resval  = &(res[(k, j)] / pivot);
+
+        self.set(k, j, selfval);
+        res.set(k, j, resval);
+      }
+
+      for uint::range(0u, dim) |l|
+      {
+        if (l != k)
+        {
+          let normalizer = self[(l, k)] / pivot;
+
+          for uint::range(k, dim) |j|
+          {
+            let selfval = &(self[(l, j)] - self[(k, j)] * normalizer);
+            let resval  = &(res[(l, j)] - res[(k, j)] * normalizer);
+
+            self.set(k, j, selfval);
+            res.set(k, j, resval);
+          }
+        }
+      }
+    }
+  }
+}
+
+impl<T:Copy> Transpose for DMat<T>
+{
+  fn transposed(&self) -> DMat<T>
+  {
+    let mut res = copy *self;
+
+    res.transpose();
+
+    res
+  }
+
+  fn transpose(&mut self)
+  {
+    let dim = self.dim;
+
+    for uint::range(1u, dim) |i|
+    {
+      for uint::range(0u, dim - 1) |j|
+      {
+        let off_i_j = self.offset(i, j);
+        let off_j_i = self.offset(j, i);
+
+        swap(self.mij, off_i_j, off_j_i);
+      }
+    }
+  }
+}
+
+impl<T: ApproxEq<T>> ApproxEq<T> for DMat<T>
+{
+  fn approx_epsilon() -> T
+  { ApproxEq::approx_epsilon::<T, T>() }
+
+  fn approx_eq(&self, other: &DMat<T>) -> bool
+  { all2(self.mij, other.mij, |a, b| a.approx_eq(b)) }
+
+  fn approx_eq_eps(&self, other: &DMat<T>, epsilon: &T) -> bool
+  { all2(self.mij, other.mij, |a, b| a.approx_eq_eps(b, epsilon)) }
+}
diff --git a/src/ndim/dvec.rs b/src/ndim/dvec.rs
index 1d1f1a38..9fc8e236 100644
--- a/src/ndim/dvec.rs
+++ b/src/ndim/dvec.rs
@@ -15,10 +15,10 @@ pub struct DVec<T>
   at: ~[T]
 }
 
-pub fn zero_with_dim<T: Zero + Copy>(dim: uint) -> DVec<T>
+pub fn zero_vec_with_dim<T: Zero + Copy>(dim: uint) -> DVec<T>
 { DVec { at: from_elem(dim, Zero::zero::<T>()) } }
 
-pub fn is_zero<T: Zero>(vec: &DVec<T>) -> bool
+pub fn is_zero_vec<T: Zero>(vec: &DVec<T>) -> bool
 { all(vec.at, |e| e.is_zero()) }
 
 // FIXME: is Clone needed?
@@ -30,7 +30,7 @@ impl<T: Copy + DivisionRing + Algebraic + Clone + ApproxEq<T>> DVec<T>
 
     for uint::range(0u, dim) |i|
     {
-      let mut basis_element : DVec<T> = zero_with_dim(dim);
+      let mut basis_element : DVec<T> = zero_vec_with_dim(dim);
 
       basis_element.at[i] = One::one();
 
@@ -49,7 +49,7 @@ impl<T: Copy + DivisionRing + Algebraic + Clone + ApproxEq<T>> DVec<T>
 
     for uint::range(0u, dim) |i|
     {
-      let mut basis_element : DVec<T> = zero_with_dim(len(self.at));
+      let mut basis_element : DVec<T> = zero_vec_with_dim(len(self.at));
 
       basis_element.at[i] = One::one();
 
diff --git a/src/ndim/nmat.rs b/src/ndim/nmat.rs
index a93c25a9..08c475eb 100644
--- a/src/ndim/nmat.rs
+++ b/src/ndim/nmat.rs
@@ -1,11 +1,11 @@
 use core::num::{One, Zero};
 use core::rand::{Rand, Rng, RngUtil};
-use core::vec::{from_elem, swap, all, all2};
 use core::cmp::ApproxEq;
 use traits::dim::Dim;
 use traits::inv::Inv;
 use traits::transpose::Transpose;
 use traits::workarounds::rlmul::{RMul, LMul};
+use ndim::dmat::{DMat, one_mat_with_dim, zero_mat_with_dim, is_zero_mat};
 use ndim::nvec::NVec;
 
 // D is a phantom type parameter, used only as a dimensional token.
@@ -14,15 +14,7 @@ use ndim::nvec::NVec;
 // using d0, d1, d2, d3 and d4 tokens are prefered.
 #[deriving(Eq, ToStr)]
 pub struct NMat<D, T>
-{
-  mij: ~[T]
-}
-
-impl<D, T: Clone> Clone for NMat<D, T>
-{
-  fn clone(&self) -> NMat<D, T>
-  { NMat{ mij: self.mij.clone() } }
-}
+{ mij: DMat<T> }
 
 impl<D: Dim, T: Copy> NMat<D, T>
 {
@@ -30,7 +22,7 @@ impl<D: Dim, T: Copy> NMat<D, T>
   { i * Dim::dim::<D>() + j }
 
   fn set(&mut self, i: uint, j: uint, t: &T)
-  { self.mij[NMat::offset::<D, T>(i, j)] = *t }
+  { self.mij.set(i, j, t) }
 }
 
 impl<D: Dim, T> Dim for NMat<D, T>
@@ -41,36 +33,23 @@ impl<D: Dim, T> Dim for NMat<D, T>
 
 impl<D: Dim, T: Copy> Index<(uint, uint), T> for NMat<D, T>
 {
-  fn index(&self, &(i, j): &(uint, uint)) -> T
-  { self.mij[NMat::offset::<D, T>(i, j)] }
+  fn index(&self, &idx: &(uint, uint)) -> T
+  { self.mij[idx] }
 }
 
 impl<D: Dim, T: Copy + One + Zero> One for NMat<D, T>
 {
   fn one() -> NMat<D, T>
-  {
-    let     dim = Dim::dim::<D>();
-    let mut res = NMat{ mij: from_elem(dim * dim, Zero::zero()) };
-    let     _1  = One::one::<T>();
-
-    for uint::range(0u, dim) |i|
-    { res.set(i, i, &_1); }
-
-    res
-  }
+  { NMat { mij: one_mat_with_dim(Dim::dim::<D>()) } }
 }
 
 impl<D: Dim, T: Copy + Zero> Zero for NMat<D, T>
 {
   fn zero() -> NMat<D, T>
-  {
-    let dim = Dim::dim::<D>();
-
-    NMat{ mij: from_elem(dim * dim, Zero::zero()) }
-  }
+  { NMat { mij: zero_mat_with_dim(Dim::dim::<D>()) } }
 
   fn is_zero(&self) -> bool
-  { all(self.mij, |e| e.is_zero()) }
+  { is_zero_mat(&self.mij) }
 }
 
 impl<D: Dim, T: Copy + Mul<T, T> + Add<T, T> + Zero>
@@ -140,87 +119,10 @@ impl<D: Dim,
 Inv for NMat<D, T>
 {
   fn inverse(&self) -> NMat<D, T>
-  {
-    let mut res : NMat<D, T> = self.clone();
-
-    res.invert();
-
-    res
-  }
+  { NMat { mij: self.mij.inverse() } }
 
   fn invert(&mut self)
-  {
-    let     dim = Dim::dim::<D>();
-    let mut res = One::one::<NMat<D, T>>();
-    let     _0T = Zero::zero::<T>();
-
-    // inversion using Gauss-Jordan elimination
-    for uint::range(0u, dim) |k|
-    {
-      // search a non-zero value on the k-th column
-      // FIXME: would it be worth it to spend some more time searching for the
-      // max instead?
-
-      // FIXME: this is kind of uggly…
-      // … but we cannot use position_between since we are iterating on one
-      // columns
-      let mut n0 = 0u; // index of a non-zero entry
-
-      while (n0 != dim)
-      {
-        if (self[(n0, k)] != _0T)
-        { break; }
-
-        n0 += 1;
-      }
-
-      assert!(n0 != dim); // non inversible matrix
-
-      // swap pivot line
-      if (n0 != k)
-      {
-        for uint::range(0u, dim) |j|
-        {
-          swap(self.mij,
-               NMat::offset::<D, T>(n0, j),
-               NMat::offset::<D, T>(k, j));
-          swap(res.mij,
-               NMat::offset::<D, T>(n0, j),
-               NMat::offset::<D, T>(k, j));
-        }
-      }
-
-      let pivot = self[(k, k)];
-
-      for uint::range(k, dim) |j|
-      {
-        // FIXME: not to putting selfal exression directly on the nuction call
-        // is uggly but does not seem to compile any more…
-        let selfval = &(self[(k, j)] / pivot);
-        let resval  = &(res[(k, j)] / pivot);
-
-        self.set(k, j, selfval);
-        res.set(k, j, resval);
-      }
-
-      for uint::range(0u, dim) |l|
-      {
-        if (l != k)
-        {
-          let normalizer = self[(l, k)] / pivot;
-
-          for uint::range(k, dim) |j|
-          {
-            let selfval = &(self[(l, j)] - self[(k, j)] * normalizer);
-            let resval  = &(res[(l, j)] - res[(k, j)] * normalizer);
-
-            self.set(k, j, selfval);
-            res.set(k, j, resval);
-          }
-        }
-      }
-    }
-  }
+  { self.mij.invert() }
 }
 
 impl<D: Dim, T:Copy> Transpose for NMat<D, T>
@@ -235,19 +137,7 @@ impl<D: Dim, T:Copy> Transpose for NMat<D, T>
   }
 
   fn transpose(&mut self)
-  {
-    let dim = Dim::dim::<D>();
-
-    for uint::range(1u, dim) |i|
-    {
-      for uint::range(0u, dim - 1) |j|
-      {
-        swap(self.mij,
-             NMat::offset::<D, T>(i, j),
-             NMat::offset::<D, T>(j, i));
-      }
-    }
-  }
+  { self.mij.transpose() }
 }
 
 impl<D, T: ApproxEq<T>> ApproxEq<T> for NMat<D, T>
@@ -256,10 +146,10 @@ impl<D, T: ApproxEq<T>> ApproxEq<T> for NMat<D, T>
   { ApproxEq::approx_epsilon::<T, T>() }
 
   fn approx_eq(&self, other: &NMat<D, T>) -> bool
-  { all2(self.mij, other.mij, |a, b| a.approx_eq(b)) }
+  { self.mij.approx_eq(&other.mij) }
 
   fn approx_eq_eps(&self, other: &NMat<D, T>, epsilon: &T) -> bool
-  { all2(self.mij, other.mij, |a, b| a.approx_eq_eps(b, epsilon)) }
+  { self.mij.approx_eq_eps(&other.mij, epsilon) }
 }
 
 impl<D: Dim, T: Rand + Zero + Copy> Rand for NMat<D, T>
diff --git a/src/ndim/nvec.rs b/src/ndim/nvec.rs
index b887908d..a18f6d3f 100644
--- a/src/ndim/nvec.rs
+++ b/src/ndim/nvec.rs
@@ -2,7 +2,7 @@ use core::num::{Zero, Algebraic};
 use core::rand::{Rand, Rng, RngUtil};
 use core::vec::{map};
 use core::cmp::ApproxEq;
-use ndim::dvec::{DVec, zero_with_dim, is_zero};
+use ndim::dvec::{DVec, zero_vec_with_dim, is_zero_vec};
 use traits::basis::Basis;
 use traits::ring::Ring;
 use traits::division_ring::DivisionRing;
@@ -164,10 +164,10 @@ Basis for NVec<D, T>
 impl<D: Dim, T: Copy + Zero> Zero for NVec<D, T>
 {
   fn zero() -> NVec<D, T>
-  { NVec { at: zero_with_dim(Dim::dim::<D>()) } }
+  { NVec { at: zero_vec_with_dim(Dim::dim::<D>()) } }
 
   fn is_zero(&self) -> bool
-  { is_zero(&self.at) }
+  { is_zero_vec(&self.at) }
 }
 
 impl<D, T: ApproxEq<T>> ApproxEq<T> for NVec<D, T>