We can leverage much more on trait to avoid repeating code and to simplify logic.

• Enhance BaseSlice (iter, ...), we should probably also have iter_rows, selectetc ...
• Have Matrix implement BaseSlice (can answer issue #14 )
• Create a BaseSliceMut trait, which add mutability
• Have Matrix and MatrixSliceMut implement BaseSliceMut

This is still a WIP, I think we should for instance be able to replace most arithmetic operations by something like (I cannot make it work yet):

impl<'a, T, S1, S2> Add<S1> for S2 
    where S1: BaseSlice<T>, S2: BaseSlice<T>,
{
// ...
}

This PR introduces new iterators (Diagonal, DiagonalMut) for Matrix/MatrixSlice/MatrixSliceMut. It is intended to resolve #53 .

Here are some benchmarks for collecting the main diagonal of a matrix:

test linalg::iter::mat_diag_iter_100_500           ... bench:         272 ns/iter (+/- 34)
test linalg::iter::mat_diag_manual_100_500         ... bench:         309 ns/iter (+/- 87)
test linalg::iter::mat_diag_iter_10_50             ... bench:          62 ns/iter (+/- 8)
test linalg::iter::mat_diag_manual_10_50           ... bench:          68 ns/iter (+/- 21)

The manual version uses unsafe get_unchecked indexing and is slower!

Following is kept in place for prosperity but this PR is no longer WIP (pending review).

This PR is currently a work in progress and the following should be addressed:

• How do we handle iteration of off-diagonals?
• How do we expose these iterators?

I'm thinking we do the following:

Add the following:

pub enum DiagOffset {
    Above(usize),
    Below(usize)
}

The iterator will then have this enum as a field - which will determine which diagonal offset we use. When creating we will have to check that Above(x) has x < mat.cols(), and Below(y) has y < mat.rows().

In terms of exposing the iterator - I think the following is a good approach:

• Add fn iter_diag(&self, k: DiagOffset) to BaseMatrix.
• Add fn iter_diag_mut(&mut self, k: DiagOffset) to BaseMatrixMut.

We will keep the diag function as returning a Vec for now - mostly to prevent breaking changes and we can reassess in future.

We will probably need to move the iterators into matrix/mod (they are currently in matrix/iter). This is so that we can construct them in matrix/slice.

This pull request is not intended as a candidate for merge (yet).

I'm working on the matrix equality macro mentioned in #67. I added an example in the examples/ subfolder, which one can run with cargo run --example mismatched_f64_elements. See also the small number of tests currently included as well. The example looks like this:

#[macro_use]
extern crate rulinalg;

fn main() {
    // The purpose of this example is to demonstrate the type of errors generated
    // by `assert_matrix_eq!` when the comparison fails (which the below call will).
    let a = matrix![1.00, 2.00;
                    3.00, 4.00];
    let b = matrix![1.01, 2.00;
                    3.40, 4.00];
    assert_matrix_eq!(a, b, comp = abs, tol = 1e-8);
}

and it gives the following output:

     Running `target/debug/examples/mismatched_f64_elements`
thread 'main' panicked at '

Matrices X and Y have 2 mismatched element pairs. The mismatched elements are listed below, in the format
(row, col): x = x[[row, col]], y = y[[row, col]].

 (0, 0): x = 1, y = 1.01. Error: 0.010000000000000009
 (1, 0): x = 3, y = 3.4. Error: 0.3999999999999999

Comparison criterion: absolute difference, defined by
    |x - y| <= 0.00000001.

', examples/mismatched_f64_elements.rs:11
note: Run with `RUST_BACKTRACE=1` for a backtrace.

Now, I'd like to give people the chance to give some input on how the macro should work. Comparing floating numbers is seriously hard. See for example this famous article. A preliminary conclusion is that there is no one-size-fits-all solution for comparing floating point numbers. Hence, it's impossible for us to have a default assert_matrix_eq!(a, b) macro that works in all cases. This is why I've opted for explicitly making the user choose which criterion to assert with (the comp = abs and tol = 1e-8 parameters).

However, I think we can provide a very conservative default assert_matrix_eq!(a, b) macro. I need to think a bit more on exactly how it's going to work, but I have the following in mind:

1. The default invocation is very strict and conservative.
2. The user simply writes assert_matrix_eq!(a, b) for his/her tests, and this might be fine if the numbers are very close to exactly equal.
3. However, in the case when the numbers are near but not sufficiently near, the macro will fail.
4. When the macro fails, the error message should explain the next steps the user should take (ideally by referring to a fleshed out part in the documentation on this) to ensure that he/she is using the correct comparison criteria.

The reasoning in step 4 again comes back to the fact that a one-size-fits-all comparison does not exist. I'm a bit inspired by how rustc guides users around compilation failures, by giving error messages that are truly helpful. Hopefully we can do the same for our users (and for our contributors).

Currently, the macro only supports absolute difference as a comparison criterion, but I want to at least support ULP-based comparison too (probably basing the default criterion on something involving ULP). I also need to write a lot more tests, as I hope this could form the cornerstone of future test assertions in rulinalg as a whole, as well as flesh out the documentation much more. Also, there should perhaps be an analogous assert_vector_eq macro?

Would really appreciate any opinions, questions, reservations or advice anyone might have!

Note: This is not a complete implementation! Intended to eventually solve #17.

My idea was to first extend the existing SVD tests to verify that the singular values were sorted. Once these were failing, I'd go ahead and implement the sort, which should restore them back to a passing state. As it happened, for the existing test data, the SVD algorithm in use already returned sorted singular values, so I constructed some new matrices for which I'd expect the algorithm not to return sorted singular values. However, it seems that for these test matrices, the SVD algorithm does not complete (i.e. it is stuck in an endless loop).

Hence, as long as the SVD algorithm does not finish, I'm also unable to complete the feature in question. I'm not at all familiar with the algorithm in use, so I can't really debug it without spending a lot of time dividing into the algorithm and how it works. Unfortunately, at the moment, that is time I don't have, so I wanted to put the PR up here, in case anybody has any idea how to fix the SVD.

See also commit messages.

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.

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?

Let me know if I should add anything, or if I need to refactor things! I didn't add a number of functions that might be useful, like is_invertible, null_space, chosen_pivots, column_space, etc.

For now, the iterator is returning Vec<&T> since it is not possible to return slices. Another possible returning type is MatrixSlice as defined in #3. The column iterator will iterate the full defined slice shape because I didn't use col_stride as the row iterator uses row_stride. There are some automatic formatting code.

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.

I still don't know for sure what the *_stride attribute does but the column iterator is working as expected without it. Didn't include raw_slice or similar for Column and ColumnMut since it is not possible to return slices for non-contiguous data. We could return Vector but that would lead to the same ambiguity discussed in #84.

Here is a much simpler but surprisingly slower version of the iterator (at least on my system). I suspect it has been slowed down by MatrxSlice initialisation, but I can't see why exactly, perhaps we're doing something wrong elsewhere. This is not intended to be merged but I wanted to show a simpler alternative without using pointers.

In this trait you use the unsafe keyword for many safe expressions. However, I found that only 4 functions are real unsafe operations (see the list below).

We need to mark unsafe operations more precisely using unsafe keyword. Keeping unsafe blocks small can bring many benefits. For example, when mistakes happen, we can locate any errors related to memory safety within an unsafe block. This is the balance between Safe and Unsafe Rust. The separation is designed to make using Safe Rust as ergonomic as possible, but requires extra effort and care when writing Unsafe Rust. Real unsafe operation list:

1. the from_raw_parts_mut()\offset()\mem::transmute() function(these are unsafe functions from the standard library)
2. the row_unchecked() function(row_unchecked() is a native unsafe function)

Hope this PR can help you. Best regards. References https://doc.rust-lang.org/nomicon/safe-unsafe-meaning.html https://doc.rust-lang.org/book/ch19-01-unsafe-rust.html

The current definition of raw_slice and raw_slice_mut creates 'a bounded reference from &self. Since the returned slice is created from a stored pointer in &self, it should be bounded by 'self lifetime instead of 'a.

With the current definitions of those methods, it is possible to cause data race with safe Rust code.

use rulinalg::matrix;
use rulinalg::matrix::BaseMatrixMut;

fn main() {
    let mut mat = matrix![0];

    let mut row = mat.row_mut(0);

    // this creates mutable aliases to the same location
    let raw_slice1 = row.raw_slice_mut();
    let raw_slice2 = row.raw_slice_mut();

    assert_eq!(raw_slice1[0], 0);
    raw_slice2[0] = 1;
    assert_eq!(raw_slice1[0], 0);
}

Use slice::iter instead of into_iter to avoid future breakage

an_array.into_iter() currently just works because of the autoref feature, which then calls <[T] as IntoIterator>::into_iter. But in the future, arrays will implement IntoIterator, too. In order to avoid problems in the future, the call is replaced by iter() which is shorter and more explicit.

A crater run showed that your crate is affected by a potential future change. See https://github.com/rust-lang/rust/pull/65819 for more information.

I want to use rulinalg on Complex Matrices, mainly the svd decomposition. I don't know much about linear algebra so I wanted to know that if I add the required traits to Complex32 numbers such as Float, would I get the right result?

Right multiplication of a matrix by a permutation permutes the columns of said matrix. The permutation that is performed however, corresponds to the transpose (inverse) of the permutation matrix, which leads to unexpected results. Consider the simple sample program:

extern crate rulinalg;

fn main()
{
    let p = rulinalg::matrix::PermutationMatrix::from_array(vec![1, 2, 0]).unwrap();
    let i = rulinalg::matrix::Matrix::<u32>::identity(3);

    println!("p.as_matrix() =\n{}", p.as_matrix());
    println!("p*i =\n{}", &p * &i);
    println!("i*p =\n{}", &i * &p);
}

The output of this program is

p.as_matrix() =
⎡0 0 1⎤
⎢1 0 0⎥
⎣0 1 0⎦
p*i =
⎡0 0 1⎤
⎢1 0 0⎥
⎣0 1 0⎦
i*p =
⎡0 1 0⎤
⎢0 0 1⎥
⎣1 0 0⎦

So left multiplying the identity matrix by a permutation gives its matrix representation (as per as_matrix()), but right multiplication results in the transpose matrix.

I can see how this happens (the same permutation is simply applied to the columns when right-multiplying), but for matrix multiplication this is not what I would expect.