Writing an eigenvalue solver in Rust for WebAssembly (opens in new tab)

A while ago, I also implemented a dense eigenvalue solver in Python following a similar approach, but found that it did not converge in O(n^3) as sometimes claimed in literature. I then read about the Divide-and-conquer eigenvalue algorithm, which did the trick. It seems to have a reasonable Wikipedia page these days: https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_...

> I found the hypot method on f64 interesting. It computes the complex norm in a numerically stable way. Supposedly. The docstring gives some very hand-wavy caveats around the numerical precision and platform-dependence.

What it does is call the libm hypot function. The reason why everything is hand-wavy is because everything about the results of libm functions is incredibly hand-wavy, and there's nothing Rust can do about that except maybe provide its own implementation, which is generally inadvisable for any numerical function. (It's not even a case of "generates the LLVM intrinsic", because LLVM doesn't have an intrinsic for libm, although LLVM intrinsics are even more hand-wavy than external libm calls because oh boy is that a rabbit hole).

> But methods like scalar_mul are just screaming for a functional refactor. I assume there’s a trait similar to Haskell’s Functor typeclass that would allow me to fmap over the entries of the Matrix and apply an arbitrary function.

The way I've solved this in my own attempts to write my own sparse BLAS routines is to boil everything down to, at the very bottom, an Entry-style class where the hard part is implementing set:

    pub fn set(&mut self, value: Scalar) {
        self.entry = match (self.entry, is_zero(value)) {
            // Change zero entry to zero entry -> do nothing.
            (None, true) => { None },
            // Nonzero to nonzero entry -> set the value.
            (Some(e), false) => {
                self.vector.values[e] = value;
                Some(e)
            },
            // Nonzero to zero entry -> delete the entry.
            (Some(e), true) => {
                self.vector.soa_vecs_mut().swap_remove(e);
                None
            },
            // Zero to nonzero entry -> add the entry.
            (None, false) => {
                self.vector.soa_vecs_mut().push((self.index, value));
                self.vector.indexes.max_index()
            }
        };
    }

All of the helper methods on vectors and matrices boil down to calling that kind of matrix, and the end result is that in my sparse LU factorization routine, the core update for pivot is actually dead simple (eliding updating all of the ancillary data structures for figuring out which element to pivot on next, which is actually most of the code):

        // When you pivot a step of LU, the matrix looks like this:
        // [1 0] * [p ǔ] = [p u]
        // [ľ I]   [0 ǎ]   [l a]
        // Solving for the unknowns in L and U, we have:
        // ľ = l / p, ǔ = u, ǎ = a - outer(ľ, ǔ)
        let pivot_val = self.values.get(row, column);
        let pivot_row = self.values.get_row(row)
            .filter(|idx, _| self.is_column_unsolved(idx))
            .to_sparse_vector();
        let pivot_column = self.values.get_column(column)
            .filter(|idx, _| self.is_row_unsolved(idx))
            .scale(pivot_val.recip())
            .to_sparse_vector();

        // Compute the outer product of l and u and update A.
        for (row, li) in pivot_column.view() {
            // We didn't update A yet, so update the requisite value in A for
            // the column.
            self.values.set(row, column, li);
            for (column, uj) in pivot_row.view() {
                self.values.entry(row, column).map(|a| li.mul_add(-uj, a));
            }
        }

> I’m sure there are good technical reasons why it’s not the case, but I was surprised there was no #[derive(AddAssign)] macro that could be added to an implementation of std::ops::Add to automatically derive the AddAssign trait.

You'd want the reverse, implementing std::ops::Add on top of std::ops::AddAssign + Clone (or maybe Copy?). The advantage of op-assignment is that you don't have to any extra allocation, so it's usually the form you want if you can do it. In my sparse BLAS lib, I only implemented AddAssign and MulAssign for sparse vectors (although I did do an impl Mul<SparseVectorView> for SparseVectorView to implement dot product). Although truth be told, I'm usually trying to a += b * c in my code anyways, and I want to generate FMAs, so I can't use the math expressions anyways.