M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Make DMat able to represent rectangular matrices. 

The code is largely untested.
This commit is contained in:
Sébastien Crozet 2013-09-05 00:01:10 +02:00
parent 8973e0d67c
commit 628066cdc8
2 changed files with 113 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

182
src/dmat.rs
View File

@ -1,106 +1,110 @@
use std::num::{One, Zero}; use std::num::{One, Zero};
use std::vec::from_elem; use std::vec::from_elem;
use std::cmp::ApproxEq; use std::cmp::ApproxEq;
use std::util;
use traits::inv::Inv; use traits::inv::Inv;
use traits::transpose::Transpose; use traits::transpose::Transpose;
use traits::rlmul::{RMul, LMul}; 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)] #[deriving(Eq, ToStr, Clone)]
pub struct DMat<N> { pub struct DMat<N> {
priv dim: uint, // FIXME: handle more than just square matrices priv nrows: uint,
priv ncols: uint,
priv mij: ~[N] priv mij: ~[N]
} }
/// Builds a matrix filled with zeros. impl<N: Zero + Clone> DMat<N> {
/// /// Builds a matrix filled with zeros.
/// # Arguments ///
/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim` /// # Arguments
/// components. /// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim`
#[inline] /// components.
pub fn zero_mat_with_dim<N: Zero + Clone>(dim: uint) -> DMat<N> { #[inline]
DMat { dim: dim, mij: from_elem(dim * dim, Zero::zero()) } pub fn new_zeros(nrows: uint, ncols: uint) -> DMat<N> {
} DMat {
nrows: nrows,
/// Tests if all components of the matrix are zeroes. ncols: ncols,
#[inline] mij: from_elem(nrows * ncols, Zero::zero())
pub fn is_zero_mat<N: Zero>(mat: &DMat<N>) -> 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<N: Clone + One + Zero>(dim: uint) -> DMat<N> {
let mut res = zero_mat_with_dim(dim);
let _1: N = One::one();
for i in range(0u, dim) {
res.set(i, i, &_1);
} }
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<N: One + Zero + Clone> DMat<N> {
/// 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<N> {
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<N: Clone> DMat<N> { impl<N: Clone> DMat<N> {
#[inline] #[inline]
fn offset(&self, i: uint, j: uint) -> uint { 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. /// Changes the value of a component of the matrix.
/// ///
/// # Arguments /// # Arguments
/// * `i` - 0-based index of the line to be changed /// * `row` - 0-based index of the line to be changed
/// * `j` - 0-based index of the column to be changed /// * `col` - 0-based index of the column to be changed
#[inline] #[inline]
pub fn set(&mut self, i: uint, j: uint, t: &N) { pub fn set(&mut self, row: uint, col: uint, val: N) {
assert!(i < self.dim); assert!(row < self.nrows);
assert!(j < self.dim); assert!(col < self.ncols);
self.mij[self.offset(i, j)] = t.clone() self.mij[self.offset(row, col)] = val
} }
/// Reads the value of a component of the matrix. /// Reads the value of a component of the matrix.
/// ///
/// # Arguments /// # Arguments
/// * `i` - 0-based index of the line to be read /// * `row` - 0-based index of the line to be read
/// * `j` - 0-based index of the column to be read /// * `col` - 0-based index of the column to be read
#[inline] #[inline]
pub fn at(&self, i: uint, j: uint) -> N { pub fn at(&self, row: uint, col: uint) -> N {
assert!(i < self.dim); assert!(row < self.nrows);
assert!(j < self.dim); assert!(col < self.ncols);
self.mij[self.offset(i, j)].clone() self.mij[self.offset(row, col)].clone()
} }
} }
impl<N: Clone> Index<(uint, uint), N> for DMat<N> { impl<N: Clone + Mul<N, N> + Add<N, N> + Zero> Mul<DMat<N>, DMat<N>> for DMat<N> {
#[inline]
fn index(&self, &(i, j): &(uint, uint)) -> N {
self.at(i, j)
}
}
impl<N: Clone + Mul<N, N> + Add<N, N> + Zero>
Mul<DMat<N>, DMat<N>> for DMat<N> {
fn mul(&self, other: &DMat<N>) -> DMat<N> { fn mul(&self, other: &DMat<N>) -> DMat<N> {
assert!(self.dim == other.dim); assert!(self.ncols == other.nrows);
let dim = self.dim; let mut res = DMat::new_zeros(self.nrows, other.ncols);
let mut res = zero_mat_with_dim(dim);
for i in range(0u, dim) { for i in range(0u, self.nrows) {
for j in range(0u, dim) { for j in range(0u, other.ncols) {
let mut acc: N = Zero::zero(); 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); 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<N>, DMat<N>> for DMat<N> {
impl<N: Clone + Add<N, N> + Mul<N, N> + Zero> impl<N: Clone + Add<N, N> + Mul<N, N> + Zero>
RMul<DVec<N>> for DMat<N> { RMul<DVec<N>> for DMat<N> {
fn rmul(&self, other: &DVec<N>) -> DVec<N> { fn rmul(&self, other: &DVec<N>) -> DVec<N> {
assert!(self.dim == other.at.len()); assert!(self.ncols == other.at.len());
let dim = self.dim; let mut res : DVec<N> = DVec::new_zeros(self.nrows);
let mut res : DVec<N> = zero_vec_with_dim(dim);
for i in range(0u, dim) { for i in range(0u, self.nrows) {
for j in range(0u, dim) { let mut acc: N = Zero::zero();
res.at[i] = res.at[i] + other.at[j] * self.at(i, j);
for j in range(0u, self.ncols) {
acc = acc + other.at[j] * self.at(i, j);
} }
res.at[i] = acc;
} }
res res
@ -129,15 +136,18 @@ RMul<DVec<N>> for DMat<N> {
impl<N: Clone + Add<N, N> + Mul<N, N> + Zero> impl<N: Clone + Add<N, N> + Mul<N, N> + Zero>
LMul<DVec<N>> for DMat<N> { LMul<DVec<N>> for DMat<N> {
fn lmul(&self, other: &DVec<N>) -> DVec<N> { fn lmul(&self, other: &DVec<N>) -> DVec<N> {
assert!(self.dim == other.at.len()); assert!(self.nrows == other.at.len());
let dim = self.dim; let mut res : DVec<N> = DVec::new_zeros(self.ncols);
let mut res : DVec<N> = zero_vec_with_dim(dim);
for i in range(0u, dim) { for i in range(0u, self.ncols) {
for j in range(0u, dim) { let mut acc: N = Zero::zero();
res.at[i] = res.at[i] + other.at[j] * self.at(j, i);
for j in range(0u, self.nrows) {
acc = acc + other.at[j] * self.at(j, i);
} }
res.at[i] = acc;
} }
res res
@ -159,9 +169,11 @@ Inv for DMat<N> {
} }
fn inplace_inverse(&mut self) -> bool { fn inplace_inverse(&mut self) -> bool {
let dim = self.dim; assert!(self.nrows == self.ncols);
let mut res = one_mat_with_dim::<N>(dim);
let _0T: N = Zero::zero(); let dim = self.nrows;
let mut res: DMat<N> = DMat::new_identity(dim);
let _0T: N = Zero::zero();
// inversion using Gauss-Jordan elimination // inversion using Gauss-Jordan elimination
for k in range(0u, dim) { for k in range(0u, dim) {
@ -197,12 +209,12 @@ Inv for DMat<N> {
let pivot = self.at(k, k); let pivot = self.at(k, k);
for j in range(k, dim) { 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); self.set(k, j, selfval);
} }
for j in range(0u, dim) { 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); res.set(k, j, resval);
} }
@ -211,12 +223,12 @@ Inv for DMat<N> {
let normalizer = self.at(l, k); let normalizer = self.at(l, k);
for j in range(k, dim) { 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); self.set(l, j, selfval);
} }
for j in range(0u, dim) { 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); res.set(l, j, resval);
} }
} }
@ -240,23 +252,23 @@ impl<N: Clone> Transpose for DMat<N> {
} }
fn transpose(&mut self) { fn transpose(&mut self) {
let dim = self.dim; for i in range(1u, self.nrows) {
for j in range(0u, self.ncols - 1) {
for i in range(1u, dim) {
for j in range(0u, dim - 1) {
let off_i_j = self.offset(i, j); let off_i_j = self.offset(i, j);
let off_j_i = self.offset(j, i); let off_j_i = self.offset(j, i);
self.mij.swap(off_i_j, off_j_i); self.mij.swap(off_i_j, off_j_i);
} }
} }
util::swap(&mut self.nrows, &mut self.ncols);
} }
} }
impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> { impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
#[inline] #[inline]
fn approx_epsilon() -> N { fn approx_epsilon() -> N {
fail!("Fix this.") fail!("This function cannot work due to a compiler bug.")
// let res: N = ApproxEq::<N>::approx_epsilon(); // let res: N = ApproxEq::<N>::approx_epsilon();
// res // res

30
src/dvec.rs
View File

@ -15,19 +15,21 @@ pub struct DVec<N> {
at: ~[N] at: ~[N]
} }
/// Builds a vector filled with zeros. impl<N: Zero + Clone> DVec<N> {
/// /// Builds a vector filled with zeros.
/// # Arguments ///
/// * `dim` - The dimension of the vector. /// # Arguments
#[inline] /// * `dim` - The dimension of the vector.
pub fn zero_vec_with_dim<N: Zero + Clone>(dim: uint) -> DVec<N> { #[inline]
DVec { at: from_elem(dim, Zero::zero()) } pub fn new_zeros(dim: uint) -> DVec<N> {
} DVec { at: from_elem(dim, Zero::zero()) }
}
/// Tests if all components of the vector are zeroes. /// Tests if all components of the vector are zeroes.
#[inline] #[inline]
pub fn is_zero_vec<N: Zero>(vec: &DVec<N>) -> bool { pub fn is_zero(&self) -> bool {
vec.at.iter().all(|e| e.is_zero()) self.at.iter().all(|e| e.is_zero())
}
} }
impl<N> Iterable<N> for DVec<N> { impl<N> Iterable<N> for DVec<N> {
@ -62,7 +64,7 @@ impl<N: Clone + Num + Algebraic + ApproxEq<N>> DVec<N> {
let mut res : ~[DVec<N>] = ~[]; let mut res : ~[DVec<N>] = ~[];
for i in range(0u, dim) { for i in range(0u, dim) {
let mut basis_element : DVec<N> = zero_vec_with_dim(dim); let mut basis_element : DVec<N> = DVec::new_zeros(dim);
basis_element.at[i] = One::one(); basis_element.at[i] = One::one();
@ -81,7 +83,7 @@ impl<N: Clone + Num + Algebraic + ApproxEq<N>> DVec<N> {
let mut res : ~[DVec<N>] = ~[]; let mut res : ~[DVec<N>] = ~[];
for i in range(0u, dim) { for i in range(0u, dim) {
let mut basis_element : DVec<N> = zero_vec_with_dim(self.at.len()); let mut basis_element : DVec<N> = DVec::new_zeros(self.at.len());
basis_element.at[i] = One::one(); basis_element.at[i] = One::one();