I'm briefly looking at rulinalg’s SliceIter. You can use the same kind of optimization used in ndarray's element iterator.

The brief idea of the "problem" is that this is efficent:

for i in rows
   for j in row(i)
       // access element at i, j

SliceIter::next() (just like ndarray's Iter::next) emulates this segmented traversal with a bit of a state machine. Unfortunately the compiler doesn't know how to turn this state machine back into two nested for loops(!)

Two main ideas I have used:

1. Use the regular [T] iterator when the matrix is contiguous. A branch on this even inside the .next() method seems to work well, the conditional is lifted out of the loop.
2. When you can't improve .next() for the non-contiguous case, you can at least special case Iterator::fold(). It is essentially the same idea as the .chain() change in rust-lang/rust#37315 ; you recurse and fold a bunch of component iterators (since each row is contiguous, you just call the slice iterator's fold on each row in turn).

In the docs we write:

Returns the row of a matrix at the given index without doing unbounds checking

Instead should read:

Returns the row of a matrix at the given index without doing any bounds checking

Current iterator functions look like this: iter_rows, iter_diag etc. We should rename these to be row_iter and diag_iter.

This change will make these functions play more nicely with IDE auto completion and fall in line with existing rust conventions.

Depends on #3 if you think in column major terms (common vector algebra habit).

From Octave, MATLAB, FORTRAN:

>> help reshape
'reshape' is a built-in function from the file libinterp/corefcn/data.cc

 -- Built-in Function: reshape (A, M, N, ...)
 -- Built-in Function: reshape (A, [M N ...])
 -- Built-in Function: reshape (A, ..., [], ...)
 -- Built-in Function: reshape (A, SIZE)
     Return a matrix with the specified dimensions (M, N, ...) whose
     elements are taken from the matrix A.

     The elements of the matrix are accessed in column-major order (like
     Fortran arrays are stored).

     The following code demonstrates reshaping a 1x4 row vector into a
     2x2 square matrix.

          reshape ([1, 2, 3, 4], 2, 2)
                =>  1  3
                    2  4

     Note that the total number of elements in the original matrix
     ('prod (size (A))') must match the total number of elements in the
     new matrix ('prod ([M N ...])').

     A single dimension of the return matrix may be left unspecified and
     Octave will determine its size automatically.  An empty matrix ([])
     is used to flag the unspecified dimension.

     See also: resize, vec, postpad, cat, squeeze.

See #25 .

We should make it clear in the documentation how the sum is computed for each case.

In this issue, I intend to maintain a list of suggestions for improvements to the assert_matrix_eq! and assert_vector_eq! macros. As more contributors get a chance to use it, I would appreciate input on missing features that would make working with it easier. Actual bugs in these macros should still be reported in their own issue.

Potential usability improvements:

• We may want to consider implementing an assert_scalar_eq! for scalar comparisons, leveraging the same comparators. Currently we have no facility for this in rulinalg, and it's a situation that comes up a lot in tests (such as testing that the determinant is almost exactly equal to an expected value).
• Print both matrices when the dimension is reasonably small. Currently it's still a little difficult to get a feel for exactly what is wrong in some cases involving small matrices.

Other improvements:

• Consider letting users provide their own comparators (low-priority unless we get requests for this).
• Avoid Copy trait bounds and avoid storing copied elements instead of references.

Would be nice to discuss the possibility of using a tuple or a struct for rows and columns in matrices constructors/structs.

For example, the following code:

pub struct Matrix<T> {
    rows: usize,
    cols: usize,
    data: Vec<T>,
}

fn from_fn<F>(rows: usize, cols: usize, mut f: F) -> Matrix<T>;
fn rows(&self) -> usize { self.rows }

Would be something like this using a tuple:

pub struct Matrix<T> {
    data: Vec<T>,
    shape: (usize, usize),
}

fn from_fn<F>(shape: (usize, usize), mut f: F) -> Matrix<T>;
fn rows(&self) -> usize { self.shape.0 }

Or something like this using a struct:

struct Shape {
    rows: usize,
    cols: usize,
}

pub struct Matrix<T> {
    data: Vec<T>,
    shape: Shape,
}

fn from_fn<F>(shape: Shape, mut f: F) -> Matrix<T>;
fn rows(&self) -> usize { self.shape.rows }

The iter_mut method for BaseMatrixMut produces a SliceIterMut. For Matrix we can do better than this as we know that the matrix data is contiguous. This has notable effects on the apply function:

current:
test linalg::matrix::apply_1000_100                ... bench:     346,624 ns/iter (+/- 29,463)

improved:
test linalg::matrix::apply_1000_100                ... bench:     274,132 ns/iter (+/- 9,047)

These numbers are from #117.

We will probably need to modify the signature of iter_mut to return an Iterator and then override it for Matrix.

Similarly to Matrix we should implement iter and iter_mut methods directly on Vector. This should simply call self.data.iter(). We can then implement IntoIterator using these methods.

As discussed in #21 we should have a Display implementation for Vector.

We have this currently for Matrix but not Vector. The implementation should be similar to this one (allowing precision to be specified) but will list the entries in an array with natural line wrapping.

As with most project names I come up with I've never been too happy with rulinalg. It doesn't really roll off the tongue.

I wanted to open this issue to give people a chance to comment on whether we should stick with this name or not. Changing the name now will be less painful than doing so later.

I think that unless there is a strong backing for a particular new name I'll be inclined to keep the name as . I plan to keep this issue open for a month or so and if there's no consensus on a new name we'll stick with rulinalg.

In places we are checking if entries in the matrix are close to zero by using: entry < n * T::min_positive_value(), where n is a small integer. It is better for us to use a larger value than the min_positive_value (which is somewhere around 1^-300). Instead we should use the new MachineEpsilon trait - or alternatively the new epsilon function in the num crate once it is released.

Places where this occurs (maybe non-exhaustive):

• matrix/mod.rs in back_substitution
• matrix/mod.rs in forward_substitution

We should simply replace the min_positive_value function with the epsilon value.

We could reduce the number of places in the code with get_unchecked calls if we were able to return a MatrixSlice and a MatrixSliceMut that consists of a diagonal.

This would only require setting row_stride equal to cols() + 1, so that the column is incremented each time. We should have assert! checks to ensure that the inputs are good, then use get_unchecked_mut if all is well.

With #47 merged we now have a macro for small matrices! This macro helps to enforce at compile time that our dimensions are correct for our input data and is stricter than the existing new constructor. Because of this we should use this for future tests and also convert existing tests to use this.

This proposed change gives an easier access point to more efficient actions on the main diagonal. This function will be equivalent to the new mat.iter_diag(DiagOffset::Main) function.

The current diag function copies the diagonal into a new Vec. This functionality can be easily retained using: mat.diag().collect::<Vec<_>>().

Currently we do self.clone() at the start. We should instead consume self to pass this choice to the user.

See #90

sum_rows seems to return a vector equal to the sum of the matrices columns instead of the matrices rows.
Similarly, sum_cols calculates the sum of the matrices rows.

Here a small code example:

let m = Matrix::new(2, 3, vec![1.0, 2.0, 3.0,
                               4.0, 5.0, 6.0]);
// => Matrix { rows: 2, cols: 3, data: [1, 2, 3, 4, 5, 6] }

m.sum_cols();
// => Vector { size: 2, data: [6, 15] }
// I expected: Vector { size: 3, data: [5, 7, 9]}

m.sum_rows();
// => Vector { size: 3, data: [5, 7, 9] }
// I expected: Vector { size: 2, data: [6, 15] }

The current Metric trait is pretty lackluster. It has a norm function which computes the euclidean norm. We could do a lot more to flesh it out. This issue is to open discussion about what this trait should support.

Some suggestions:

• Adding L-norms.
  • How should we define L, usize, float, enum?
  • How would the L^infinity norm be defined?
• Should we directly allow distances (i.e. metrics...). Perhaps we should split out the norms into a new trait and use this trait in the definition of Metric: pub trait Metric : Norm.
• We could cover a lot here. Where should we define our limit?

In rusty-machine we often want to define a metric to be used as a prediction/classification error. For example we can choose the metric in k-means used to compute distance between elements in each set. It would be nice if we could provide a trait in rulinalg which we can be freely generic over in rulinalg.

Any work defined from this issue would be a breaking change - we should aim to include these in a 0.4 release.

One would expect the matrix! macro to allow the creation of empty matrices: let empty_mat = matrix![]. However currently this fails to compile as there are multiple viable interpretations of the macro.

I'm verifying the results of the FD port I have written using rulinalg, which consists of the following two files:

• https://gitlab.com/bright-star/electro/blob/master/examples/empty_rect.rs
• https://gitlab.com/bright-star/electro/blob/master/src/lib.rs

Running that example cargo run --example empty_rect with nightly gives an LUP-related error:

thread 'main' panicked at 'called `Result::unwrap()` on an `Err` value: 
Error { kind: DivByZero, error: StringError("Singular matrix found in LUP decomposition. A value in the diagonal of U == 0.0.") }', ../src/libcore/result.rs:799

It also prints out the matrix and vector it was trying to solve with, which I verified against the MATLAB source I am porting with the same params: https://github.com/bright-star/computational-em/blob/master/finite-difference/e1_uniform.m So I know that this a solvable system.

For my sanity's sake: these are identical matrices, right?

-4   1   0   1   0   0   0   0   0
 1  -4   1   0   1   0   0   0   0
 0   1  -4   0   0   1   0   0   0
 1   0   0  -4   1   0   1   0   0
 0   1   0   1  -4   1   0   1   0
 0   0   1   0   1  -4   0   0   1
 0   0   0   1   0   0  -4   1   0
 0   0   0   0   1   0   1  -4   1
 0   0   0   0   0   1   0   1  -4

-4   1   0   1   0   0   0   0   0
 1  -4   1   0   1   0   0   0   0
 0   1  -4   0   0   1   0   0   0
 1   0   0  -4   1   0   1   0   0
 0   1   0   1  -4   1   0   1   0
 0   0   1   0   1  -4   0   0   1
 0   0   0   1   0   0  -4   1   0
 0   0   0   0   1   0   1  -4   1
 0   0   0   0   0   1   0   1  -4

Vector does not implement the Hash trait, which makes it - for example - unusable as the key of a HashMap

Following requirements:

• Should be behind a feature flag (to avoid new users running into dependency hell by default).
• Hopefully will fail at compile time if bindings are not found. (Prefer this to silently continuing without failure).

It would also be good to include benchmarks for matrix multiplication across different sizes.

From @AtheMathmo on July 6, 2016 22:10

We currently have a transpose(&self) -> Matrix<T> method. It would be nice if we could do this operation in place, i.e. in_place_transpose(&mut self).

We can achieve this by swapping places in the data. We could do this fairly easily by utilizing the existing slice::swap method. This is acceptable but will include unnecessary bound checking, to speed things up we could use the underlying technique from that swap method:

unsafe {
    let pa: *mut T = self.get_unchecked_mut(i);
    let pb: *mut T = self.get_unchecked_mut(j);
    ptr::swap(pa, pb);
}

It would be nice to have this method in place for both Matrix and MatrixSliceMut. The Matrix implementation could simply call self.as_mut_slice().in_place_transpose().

To resolve this issue

1. In place transposition for square Matrix.
2. In place transposition for square MatrixSliceMut.
3. Tests for both.

Bonus points

1. Using ptr::swap to avoid bounds checking.
2. Explore vectorization and cache-line utilization.
3. Explore more advanced techniques.
4. Use permutations to transpose non-square matrices.

Copied from original issue: AtheMathmo/rusty-machine#97

This example would not compile:

let mut x = Vector::new(vec![0; 10]);

for i in 0..10 {
     x[i] = i;
}

The compiler complains that "cannot assign to immutable indexed content". This is because only Index is implemented for Vector<T>. Looking in the docs, IndexMut is clearly missing. By now, the workaround I have found was:

let mut x = Vector::new(vec![0; 10]);

for i in 0..10 {
     x.mut_data()[i] = i;
}

which does compile.

Similar to #15, would be nice to add clBLAS and clSPARSE bindings, which are OpenCL libraries freely available for general public developed by AMD.
Adding bindings to those libs would require much research and effort but I think it is something to consider in the long term.

More info:
Benchmarking clBLAS SGEMM on the AMD FirePro™ Professional Graphics Card
Open Source clSPARSE Beta Released, Library Accelerates Performance
AutoGemm: A flexible and high-performance solution to multiplying matrices

From @zackmdavis on April 14, 2016 15:52

There exist more efficient algorithms, but just using an eigenbasis of AᵀA is probably okay to start with? I understand this should also get us principal components analysis pretty easily.

Copied from original issue: AtheMathmo/rusty-machine#40

For any vector v, v is equal to itself.
This is then full equality and the vector needs to implement the trait Eq.
Right now vectors only implement PartialEq

Apart from logical reason to do this, some traits require Eq (for example the trait Ord).

We will introduce a new generic function to enable slicing in #18 . This means reslice can be phased out.

We can get a slight performance gain in the Matrix::variance function.

let v = Vector::new(t);
variance = variance + &(&v - &mean).elemul(&(&v - &mean));

Here we compute v - mean twice. Instead we should compute it once and then modify.

// Computes diff using the space allocated for `t`
let diff = Vector::new(t) - &mean;
variance = variance + diff.elemul(&diff);

This may well be optimized away already (?) - but we should make it explicit.

From @AtheMathmo on April 22, 2016 2:0

We have added row iterators with #48, #49 . These were easy to do because our matrices are row-major. as mentioned in #48 it would be nice to have column iterators too.

This ticket is to track the implementation of column iterators. The first step should be to write up a specification for how exactly these should look.

Copied from original issue: AtheMathmo/rusty-machine#52

My program has been panicking when I invoke solve_l_triangular on this:

⎡                    1.68933665291405 0.0000000000000000027626478732696488⎤
⎣                  0.0210783862932372                   1.7853800754297995⎦

Ok, I know! This is technically not lower triangular. But could not rulinalg be a bit more forgiving? Besides, there is no way to tell the code "I know what I am doing. Carry on!".

Alternatively, there is no easy way to force a matrix to be lower triangular, is there?

Somewhat related to #17 and #46.

There are two known bugs in the SVD implementation.

• When computing the bidiagonal blocks we use the wrong criteria.

// This is wrong:
diag_abs_sum = T::min_positive_value() *
                                   (b_ii.abs() + *b.get_unchecked([i + 1, i + 1]));

// Instead:
diag_abs_sum = T::min_positive_value() *
                                   (b_ii.abs() + b.get_unchecked([i + 1, i + 1]).abs());

• We zero the off-diagonals incorrectly

We also zero the off-diagonals incorrectly before running the Golub-Kahan step. Currently we use a givens rotation to zero the off diagonal element in b. This is essentially equivalent to setting the value to zero directly... Instead we should use a (implicit) Given's rotation on the whole matrix and update the orthogonal matrices in turn.

• We do not need to transpose b if it was flipped before

Even if we transposed at the start, the matrix b is diagonal by definition and so we gain nothing by transposing it.

• Our stopping criterion is far too harsh

We frequently check that we are below T::min_positive_value(). This makes the algorithm take a very long time to converge. Instead we should use a larger (but still small) value.

In addition to these known issues we should explore the following:

• Should the algorithm order the singular values automatically?
• Why does the algorithm return negative singular values?

These last two may not be bugs but rather things we clean up after the algorithm.

This is a mirror to rusty-machine/109.

These methods can fail and it would be good for this to be recognized in the API.

The mean method may fail if the input Matrix is empty - though it may make sense to simply return an empty mean vector in response.

The variance method may fail if the input Matrix has only one row/column of data in the working axis. In this case we will end up dividing by zero (as this method computes the unbiased sample variance). For now we should fail in this case and perhaps later add in a new method for the biased sample variance which will simply return 0 in this case.

It would be good to add a slicing equivalent for Vector (similar to MatrixSlice for Matrix).

For Vector we can probably do this in a nicer way than Matrix. We can implement a fat pointer type to represent the underlying data in a Vector - and then implement the relevant methods on a reference to this type (similar to the Vec - &[T] relationship in std).

Using Deref/Borrow/AsRef we can cut down on code by implementing all functions taking &self on the VectorSlice type.

It would be acceptable to simply mirror the current behaviour for matrices - but we should try to get the above working.

Currently the singular values are not ordered though it is convention to do so.

We are currently using the Golub-Reinsch scheme to compute the singular value decomposition. I don't think there is a way to implicitly order the singular values so we should do a standard sort whilst reordering the columns and rows of the unitary matrices as appropriate.

Currently it is a bit of a headache to convert between types for matrices. For example say we have the following 2x2 matrix of integers:

[1, 2;
 3, 4]

And we want to convert this to a matrix of floats. To do this currently in rulinalg we need something like this:

use num::ToPrimitive;

let (n, m) = (mat.rows(), mat.cols());
let float_data = mat.into_iter().map(|x| x.to_f64()).collect::<Vec<_>>();
let float_mat = Matrix::new(n, m, float_data);

Instead of this we could leverage the NumCast trait. This would allow us to convert between the primitives. We could add a new function like the following:

impl<T: ToPrimitive> Matrix<T> {
    pub fn try_into<U: NumCast>(self) -> Result<Matrix<U>, ...> {
        let (n, m) = (self.rows(), self.cols());
        let new_data = self.into_iter().map(|x| try!(U::from(x).ok_or(...))).collect::<Vec<_>>();
        Matrix::new(n, m, float_data)
    }
}

This will allow us to cast between the primitive types, or any custom types that can go between primitives using the ToPrimitive and NumCast traits. We can get a more complete coverage if the TryInto and TryFrom traits are stabilized and implemented between primitives.

Hi.

When I used this library in my own machine learning project, I tried to implement my own method of typing the values of vector and matrix, such as :

`pub fn show_vector_f64(vec: Vector) {

print!("[");
for x in (vec.data()) {
    print!("{},", x);
}
print!("]\n");

}
pub fn show_matrix_f64(mat: Matrix) {

let row_num = mat.rows();
let col_num = mat.cols();
print!("[\n");
for i in 0..row_num {
    print!("[");
    for j in 0..col_num {
        print!("{},", mat.data()[i*j]);
    }
    print!("]\n");
}
print!("]\n");

}`

I'm not sure if it is useful to add them or some similar but more clear code in the library. However it is hard to implement a fmt::Display method to the generic type of Matrix, Vector, but Matrix and Vector is often used so... maybe it's helpful.

This issue is to move discussion from #84 and #78 in one place. I'll fill out this issue properly soon, but primarily we want to discuss how best to handle the return values on the row and column access functions.

Hello,

Does anyone mind if I send a sparse matrix pull request? I am implement it and should be done in ~2 weeks.

Thanks

We added the MachineEpsilon trait to get generic access to the f32::EPSILON and f64::EPSILON constants. But these are now available via the Float trait in the num crate.

We should remove the MachineEpsilon trait and replace uses of it with the num crate version.

This requires num version 0.1.36

From @AtheMathmo on April 17, 2016 0:55

Would be nice if we could get things running on more than once core! I've been playing around with getting this working for matrix multiplication for a while. Now that we have MatrixSlice we can get something decent working. My initial tests produced the following benchmarks:

test linalg::matrix::mat_mul_128_100 ... bench: 221,813 ns/iter (+/- 28,576)
test linalg::matrix::mat_paramul_128_100 ... bench: 213,257 ns/iter (+/- 16,667)
test linalg::matrix::mat_blasmul_128_100 ... bench: 107,305 ns/iter (+/- 14,451)

test linalg::matrix::mat_mul_128_1000 ... bench: 1,994,442 ns/iter (+/- 79,774)
test linalg::matrix::mat_paramul_128_1000 ... bench: 1,147,764 ns/iter (+/- 136,592)
test linalg::matrix::mat_blasmul_128_1000 ... bench: 996,405 ns/iter (+/- 109,778)

test linalg::matrix::mat_mul_128_10000 ... bench: 21,185,583 ns/iter (+/- 794,584)
test linalg::matrix::mat_paramul_128_10000 ... bench: 11,687,473 ns/iter (+/- 638,582)
test linalg::matrix::mat_blasmul_128_10000 ... bench: 10,278,981 ns/iter (+/- 973,273)

test linalg::matrix::mat_mul_128_100000 ... bench: 210,618,866 ns/iter (+/- 4,908,516)
test linalg::matrix::mat_paramul_128_100000 ... bench: 112,120,346 ns/iter (+/- 6,052,281)
test linalg::matrix::mat_blasmul_128_100000 ... bench: 102,699,089 ns/iter (+/- 9,024,207)

We get roughly a 2x increase in performance (on my sub-par laptop) when using the parallel implementation (that is currently on the paramul branch). The above results are for f32 only. For f64 the largest benchmark produces:

test linalg::matrix::mat_mul_f64_128_100000 ... bench: 445,007,480 ns/iter (+/- 71,323,075)
test linalg::matrix::mat_paramul_f64_128_100000 ... bench: 254,693,413 ns/iter (+/- 57,254,546)

This is a promising start. This issue will track progress.

Copied from original issue: AtheMathmo/rusty-machine#44

It would be nice if we could implement a get_row method to get a slice of a single row from a Matrix. Something like the following:

impl Matrix<T> {
    pub fn get_row(&self, idx: usize) -> &[T] {}

    pub fn get_row_mut(&mut self, idx: usize) -> &mut [T] {}
}

The contents of the method will need to be unsafe as we'll have to play with pointers a little within the slices. The above methods should do safety checks to ensure the indices are within the right range.

It may also be useful to implement unsafe get_row_unchecked and get_row_unchecked_mut which can be reused in the above methods.

The code in the iter_rows methods can likely be adapted to produce the get_row methods etc.

Rulinalg has implementations of a number of decompositions. Currently these are returned as Result<tuple, Error>, where tuple is a tuple of Matrix<T>. This is in itself useful, as it gives the user easy access to the decomposed parts. However, often you are not necessarily so interested in the decomposed parts themselves, but rather what you can do with them.

For example, in the case of LU or QR and a linear system Ax = b, you may want to perform the decomposition once, and then use it many times over to solve for different instances of b. In this case, as a user, you don't want to have to perform the transposition of Q and the solution of the triangular system. Rather, it would be convenient if you could do something along the lines of (ignoring error handling for to make my point more clear):

let qr = A.qr();
for b in right_hand_sides {
    let x = qr.solve(b);
    do_something(x);
}

This is possible if qr() returns a struct QR that has a specialized solve method, as well as accessors to the decomposed parts Q and R. Similarly for LU. In my opinion, you would not lose anything over the current API, but you'd make it easier to use, and it could also open up for optimizations in how the decompositions are stored, and how they are used to perform certain actions.

For SVD, you could have convenient methods to retrieve the singular values, perhaps for least squares solution etc. In addition to accessors to the decomposed parts, of course.

Eigen (C++ library) does something like this. I'll be quick to add that I do not want rulinalg to mirror Eigen's API too closely, as it is... perhaps overly complex at times.

Do you think this style of API could be suitable for rulinalg? I'd be interested to hear your opinions!

Currently some of our tests look like this:

assert!(!mat.data()
    .iter()
    .zip(recovered.data().iter())
    .any(|(&x, &y)| (x - y).abs() > 1e-10));

This is functional, and it works, but it has some issues:

• It is not immediately obvious what it does, and as such subtle errors may easily creep in.
• It is verbose.
• If the assertion fails, it's very hard/inconvenient to determine what actually went wrong, because it is boolean.
• Most of our tests currently use arbitrary thresholds (*)

It would be immensely useful for us, if we instead could write something along the lines of

// Elementwise absolute difference comparison
assert_matrix_eq!(mat, recovered, tol = 1e-10, kind = absolute);

// Frobenius norm comparison, perhaps?
assert_matrix_eq!(mat, recovered, tol = 1e-10, kind = frobenius);

// Default comparison (X amount of Units in Last Place (ULP) perhaps?)
assert_matrix_eq!(mat, recovered);

With such a macro, we could have much more detailed feedback when the assertion fails. For example, for a small matrix:

Expected:
1.00000, 0.00000, 0.00000,
0.00000, 1.00000, 0.00000,
0.00000, 0.00000, 1.00000

Actual:
1.00000, 0.00000, 0.00000,
0.00031, 1.00001, 0.00002,
0.00000, 0.00312, 1.00812

or for a larger matrix (bigger than 5x5 or so?)

397 out of 400 matrix elements compared equal.
The following elements did not meet comparison criteria:
(15,  3) - expected: 103.151351, actual: -99.001141
(14,  2) - expected:  13.000000, actual:  12.500000
( 3,  3) - expected:   0.000000, actual: - 3.153135

While this is immensely useful for the developers of rulinalg, keep in mind that our users also need to write test cases for their applications, and as such it would be very useful also for them.

Note that I have not considered the above very thoroughly. I rather wanted to throw out some ideas, which might inspire other good ideas.

I'd like to start experimenting with this and see what I can come up with, when I have time.

(*) On arbitrary thresholds: this obviously isn't fixed by a macro. Instead, at some point we need to try to design our tests to use theoretically proven bounds for errors. This is sort of next level though, but it has to be done if rulinalg wants to be a serious alternative in the longer run.

It would be valuable to benchmark rulinalg. It would also be great to do comparisons against other libraries, both utilizing LAPACK and standalone - to see how well rulinalg holds its own.

Things to benchmark:

• Matrix Multiplication
• Matrix Inversion
• Eigen Decomposition
• Singular Value Decomposition

It would be good to show rulinalg vs. LAPACK (or some library consuming, i.e. scipy). Even just having these benchmarks without a direct comparison would be very valuable.

It would be nice to produce a BENCHMARKS.md file in the root project directory. I see this file containing a table (or set of tables) with the details of the machine used to produce the table clearly visible. The table should have rows per function benchmarked with columns for each framework tested (rulinalg in first column for each).

In slice.rs there are function signatures such as elemul that require the same type for both sides of the operation. This leads to pretty silly code as you can't do an elemul operation with a mut slice on the one side and a normal matrix on the other.

(This issue could probably be tagged "easy," and I'm willing to help review PRs if needed.)

Gulp...

Is this worth doing?

It will allow us to integrate more nicely with BLAS/LAPACK - and we get things like free transposition as a result. We currently have a lot of the framework to do this but will need to make a lot of changes - especially to the decomposition algorithms and other areas optimized for row-major.

Note also that our iter_rows usage will become difficult...

This issue should track potential pitfalls and criteria. We'll make a new issue to track the actual implementation if it happens.

We should try to move towards a standardized formatting. We should almost definitely use rustfmt and use a project config file to enforce the standard.

As suggested by @Andlon in #57 we could also check the formatting using Travis.

I'm happy to use the rustfmt defaults but I think it would be best to have the config file in place to override contributors global configs.

If you call

matrix.swap_rows(index, index);

I believe it currently invokes undefined behavior, as it will perform mem::swap on elements that reside in the same memory location. The same is currently not true for swap_cols(), as it uses ptr::swap.

Moreover, even if this might work with the current Rust compiler, it still performs a full swap when we could let it do a no-op. If swap_rows(index, index) was a no-op, we could simplify some calling code which currently explicitly checks that we're not swapping a row with itself (see for example the lup_decomp code). This is also true for swap_cols().

I think a fix should at the very least only perform the swap if a != b in swap_rows() and swap_cols(). Even better if we can rewrite it to avoid unsafe code without a performance penalty (I don't know if this is doable).