From 628066cdc8b4ec561622eb61c3725bc8ea23da4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Thu, 5 Sep 2013 00:01:10 +0200 Subject: [PATCH] Make DMat able to represent rectangular matrices. The code is largely untested. --- src/dmat.rs | 182 ++++++++++++++++++++++++++++------------------------ src/dvec.rs | 30 +++++---- 2 files changed, 113 insertions(+), 99 deletions(-) diff --git a/src/dmat.rs b/src/dmat.rs index 150b9326..881d315c 100644 --- a/src/dmat.rs +++ b/src/dmat.rs @@ -1,106 +1,110 @@ use std::num::{One, Zero}; use std::vec::from_elem; use std::cmp::ApproxEq; +use std::util; use traits::inv::Inv; use traits::transpose::Transpose; use traits::rlmul::{RMul, LMul}; -use dvec::{DVec, zero_vec_with_dim}; +use dvec::DVec; -/// Square matrix with a dimension unknown at compile-time. +/// Matrix with dimensions unknown at compile-time. #[deriving(Eq, ToStr, Clone)] pub struct DMat { - priv dim: uint, // FIXME: handle more than just square matrices + priv nrows: uint, + priv ncols: uint, priv mij: ~[N] } -/// Builds a matrix filled with zeros. -/// -/// # Arguments -/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim` -/// components. -#[inline] -pub fn zero_mat_with_dim(dim: uint) -> DMat { - DMat { dim: dim, mij: from_elem(dim * dim, Zero::zero()) } -} - -/// Tests if all components of the matrix are zeroes. -#[inline] -pub fn is_zero_mat(mat: &DMat) -> bool { - mat.mij.iter().all(|e| e.is_zero()) -} - -/// Builds an identity matrix. -/// -/// # Arguments -/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim` -/// components. -#[inline] -pub fn one_mat_with_dim(dim: uint) -> DMat { - let mut res = zero_mat_with_dim(dim); - let _1: N = One::one(); - - for i in range(0u, dim) { - res.set(i, i, &_1); +impl DMat { + /// Builds a matrix filled with zeros. + /// + /// # Arguments + /// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim` + /// components. + #[inline] + pub fn new_zeros(nrows: uint, ncols: uint) -> DMat { + DMat { + nrows: nrows, + ncols: ncols, + mij: from_elem(nrows * ncols, Zero::zero()) + } } - res + /// Tests if all components of the matrix are zeroes. + #[inline] + pub fn is_zero(&self) -> bool { + self.mij.iter().all(|e| e.is_zero()) + } } +// FIXME: add a function to modify the dimension (to avoid useless allocations)? + +impl DMat { + /// Builds an identity matrix. + /// + /// # Arguments + /// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim` + /// components. + #[inline] + pub fn new_identity(dim: uint) -> DMat { + let mut res = DMat::new_zeros(dim, dim); + + for i in range(0u, dim) { + let _1: N = One::one(); + res.set(i, i, _1); + } + + res + } +} + + impl DMat { #[inline] fn offset(&self, i: uint, j: uint) -> uint { - i * self.dim + j + i * self.ncols + j } /// Changes the value of a component of the matrix. /// /// # Arguments - /// * `i` - 0-based index of the line to be changed - /// * `j` - 0-based index of the column to be changed + /// * `row` - 0-based index of the line to be changed + /// * `col` - 0-based index of the column to be changed #[inline] - pub fn set(&mut self, i: uint, j: uint, t: &N) { - assert!(i < self.dim); - assert!(j < self.dim); - self.mij[self.offset(i, j)] = t.clone() + pub fn set(&mut self, row: uint, col: uint, val: N) { + assert!(row < self.nrows); + assert!(col < self.ncols); + self.mij[self.offset(row, col)] = val } /// Reads the value of a component of the matrix. /// /// # Arguments - /// * `i` - 0-based index of the line to be read - /// * `j` - 0-based index of the column to be read + /// * `row` - 0-based index of the line to be read + /// * `col` - 0-based index of the column to be read #[inline] - pub fn at(&self, i: uint, j: uint) -> N { - assert!(i < self.dim); - assert!(j < self.dim); - self.mij[self.offset(i, j)].clone() + pub fn at(&self, row: uint, col: uint) -> N { + assert!(row < self.nrows); + assert!(col < self.ncols); + self.mij[self.offset(row, col)].clone() } } -impl Index<(uint, uint), N> for DMat { - #[inline] - fn index(&self, &(i, j): &(uint, uint)) -> N { - self.at(i, j) - } -} - -impl + Add + Zero> -Mul, DMat> for DMat { +impl + Add + Zero> Mul, DMat> for DMat { fn mul(&self, other: &DMat) -> DMat { - assert!(self.dim == other.dim); + assert!(self.ncols == other.nrows); - let dim = self.dim; - let mut res = zero_mat_with_dim(dim); + let mut res = DMat::new_zeros(self.nrows, other.ncols); - for i in range(0u, dim) { - for j in range(0u, dim) { + for i in range(0u, self.nrows) { + for j in range(0u, other.ncols) { let mut acc: N = Zero::zero(); - for k in range(0u, dim) { + for k in range(0u, self.ncols) { acc = acc + self.at(i, k) * other.at(k, j); } - res.set(i, j, &acc); + res.set(i, j, acc); } } @@ -111,15 +115,18 @@ Mul, DMat> for DMat { impl + Mul + Zero> RMul> for DMat { fn rmul(&self, other: &DVec) -> DVec { - assert!(self.dim == other.at.len()); + assert!(self.ncols == other.at.len()); - let dim = self.dim; - let mut res : DVec = zero_vec_with_dim(dim); + let mut res : DVec = DVec::new_zeros(self.nrows); - for i in range(0u, dim) { - for j in range(0u, dim) { - res.at[i] = res.at[i] + other.at[j] * self.at(i, j); + for i in range(0u, self.nrows) { + let mut acc: N = Zero::zero(); + + for j in range(0u, self.ncols) { + acc = acc + other.at[j] * self.at(i, j); } + + res.at[i] = acc; } res @@ -129,15 +136,18 @@ RMul> for DMat { impl + Mul + Zero> LMul> for DMat { fn lmul(&self, other: &DVec) -> DVec { - assert!(self.dim == other.at.len()); + assert!(self.nrows == other.at.len()); - let dim = self.dim; - let mut res : DVec = zero_vec_with_dim(dim); + let mut res : DVec = DVec::new_zeros(self.ncols); - for i in range(0u, dim) { - for j in range(0u, dim) { - res.at[i] = res.at[i] + other.at[j] * self.at(j, i); + for i in range(0u, self.ncols) { + let mut acc: N = Zero::zero(); + + for j in range(0u, self.nrows) { + acc = acc + other.at[j] * self.at(j, i); } + + res.at[i] = acc; } res @@ -159,9 +169,11 @@ Inv for DMat { } fn inplace_inverse(&mut self) -> bool { - let dim = self.dim; - let mut res = one_mat_with_dim::(dim); - let _0T: N = Zero::zero(); + assert!(self.nrows == self.ncols); + + let dim = self.nrows; + let mut res: DMat = DMat::new_identity(dim); + let _0T: N = Zero::zero(); // inversion using Gauss-Jordan elimination for k in range(0u, dim) { @@ -197,12 +209,12 @@ Inv for DMat { let pivot = self.at(k, k); for j in range(k, dim) { - let selfval = &(self.at(k, j) / pivot); + let selfval = self.at(k, j) / pivot; self.set(k, j, selfval); } for j in range(0u, dim) { - let resval = &(res.at(k, j) / pivot); + let resval = res.at(k, j) / pivot; res.set(k, j, resval); } @@ -211,12 +223,12 @@ Inv for DMat { let normalizer = self.at(l, k); for j in range(k, dim) { - let selfval = &(self.at(l, j) - self.at(k, j) * normalizer); + let selfval = self.at(l, j) - self.at(k, j) * normalizer; self.set(l, j, selfval); } for j in range(0u, dim) { - let resval = &(res.at(l, j) - res.at(k, j) * normalizer); + let resval = res.at(l, j) - res.at(k, j) * normalizer; res.set(l, j, resval); } } @@ -240,23 +252,23 @@ impl Transpose for DMat { } fn transpose(&mut self) { - let dim = self.dim; - - for i in range(1u, dim) { - for j in range(0u, dim - 1) { + for i in range(1u, self.nrows) { + for j in range(0u, self.ncols - 1) { let off_i_j = self.offset(i, j); let off_j_i = self.offset(j, i); self.mij.swap(off_i_j, off_j_i); } } + + util::swap(&mut self.nrows, &mut self.ncols); } } impl> ApproxEq for DMat { #[inline] fn approx_epsilon() -> N { - fail!("Fix this.") + fail!("This function cannot work due to a compiler bug.") // let res: N = ApproxEq::::approx_epsilon(); // res diff --git a/src/dvec.rs b/src/dvec.rs index a79a700c..1dc8bd68 100644 --- a/src/dvec.rs +++ b/src/dvec.rs @@ -15,19 +15,21 @@ pub struct DVec { at: ~[N] } -/// Builds a vector filled with zeros. -/// -/// # Arguments -/// * `dim` - The dimension of the vector. -#[inline] -pub fn zero_vec_with_dim(dim: uint) -> DVec { - DVec { at: from_elem(dim, Zero::zero()) } -} +impl DVec { + /// Builds a vector filled with zeros. + /// + /// # Arguments + /// * `dim` - The dimension of the vector. + #[inline] + pub fn new_zeros(dim: uint) -> DVec { + DVec { at: from_elem(dim, Zero::zero()) } + } -/// Tests if all components of the vector are zeroes. -#[inline] -pub fn is_zero_vec(vec: &DVec) -> bool { - vec.at.iter().all(|e| e.is_zero()) + /// Tests if all components of the vector are zeroes. + #[inline] + pub fn is_zero(&self) -> bool { + self.at.iter().all(|e| e.is_zero()) + } } impl Iterable for DVec { @@ -62,7 +64,7 @@ impl> DVec { let mut res : ~[DVec] = ~[]; for i in range(0u, dim) { - let mut basis_element : DVec = zero_vec_with_dim(dim); + let mut basis_element : DVec = DVec::new_zeros(dim); basis_element.at[i] = One::one(); @@ -81,7 +83,7 @@ impl> DVec { let mut res : ~[DVec] = ~[]; for i in range(0u, dim) { - let mut basis_element : DVec = zero_vec_with_dim(self.at.len()); + let mut basis_element : DVec = DVec::new_zeros(self.at.len()); basis_element.at[i] = One::one();