Allow reducing temporary allocations by making the ipvec function public.

[Voidheart88]

Hi,

First of all, I'm an electronics engineer first and a Rust enthusiast second, so maybe I'm just using your library wrong.

I'm writing a SPICE simulator and chose rsparse as my solver due to its simplicity. In each simulation step, I need to evaluate the circuit elements to get their conductance values and current values. The conductance values populate the system matrix a (a sparse matrix), and the current values go to the known values vector.

After populating the matrix, I use rsparse to solve the system. Here's a minimal example:

use rsparse::{data::{Sprs, Trpl}, lusol};

fn main() {
    // Initial matrix triples
    let mut triples = Trpl::new();
    triples.append(0, 0, 2.0);
    triples.append(1, 1, 2.0);
    triples.append(2, 2, 2.0);

    // Get System Matrix
    let mut system_matrix = Sprs::new_from_trpl(&triples);
    // Known vector
    let mut known_vector = vec![1.0,2.0,3.0];

    // Solve system
    lusol(&system_matrix, &mut known_vector, 1, 1e-6);
    assert_eq!(known_vector,vec![0.5,1.0,1.5]);

    // matrix triples change since the corresponding circuit elements change.
    // Note, that the sparsity pattern does not change
    let mut new_triples = Trpl::new();
    new_triples.append(0, 0, 1.0);
    new_triples.append(1, 1, 1.0);
    new_triples.append(2, 2, 1.0);

    // from tripl, since the memory is already allocated
    system_matrix.from_trpl(&new_triples);

    // Known vector changes too after evaluating circuit elements.
    let mut known_vector = vec![2.0,3.0,4.0];

    // solve again
    lusol(&system_matrix, &mut known_vector, 1, 1e-6);

    assert_eq!(known_vector,vec![2.0,3.0,4.0]);
    // ...
    // Carry on until the Simulation is finished
}

I've benchmarked my simulator and compared it with some publicly available alternatives (notably ngspice). As expected, mine is slower. After analyzing memory consumption, particularly the allocations, I noticed some temporary allocations and avoidable calculations occur inside lusol.

As far as I've noticed, each time I call lusol, the following happens:

An allocation of a vector filled with zeros: let mut x = vec![0.; a.n];
A symbolic analysis: s = sqr(a, order, false);
Since the sparsity pattern doesn't change, I've tried to avoid repeated allocations of x and s by calling lsolve and usolve manually and reusing s:

use rsparse::data::{Sprs, Trpl};

fn main() {
    // Initial matrix triples
    let mut triples = Trpl::new();
    triples.append(0, 0, 2.0);
    triples.append(1, 1, 2.0);
    triples.append(2, 2, 2.0);

    // Get System Matrix
    let mut system_matrix = Sprs::new_from_trpl(&triples);
    // Known vector
    let mut known_vector = vec![1.0,2.0,3.0];

    //allocate Workspace
    let mut x = vec![0.; system_matrix.n]; 
    // calculate the symbolic analysis
    let mut s= rsparse::sqr(&system_matrix, 1, false); //1 = LU factorisation

    // Solve system :
    let n = rsparse::lu(&system_matrix, &mut s, 1e-6);
    rsparse::ipvec(system_matrix.n, &n.pinv, &mut known_vector, &mut x[..]);
    rsparse::lsolve(&n.l, &mut x); // x = L\x
    rsparse::usolve(&n.u, &mut x[..]); // x = U\x
    rsparse::ipvec(system_matrix.n, &s.q, &x[..], &mut known_vector[..]);


    assert_eq!(known_vector,vec![0.5,1.0,1.5]);

    // matrix triples change since the corresponding circuit elements change.
    // Note, that the sparsity pattern does not change
    let mut new_triples = Trpl::new();
    new_triples.append(0, 0, 1.0);
    new_triples.append(1, 1, 1.0);
    new_triples.append(2, 2, 1.0);

    // from tripl, since the memory is already allocated
    system_matrix.from_trpl(&new_triples);

    // Known vector changes too after evaluating circuit elements.
    let mut known_vector = vec![2.0,3.0,4.0];

    // solve again
    let n = rsparse::lu(&system_matrix, &mut s, 1e-6); //reuse s
    rsparse::ipvec(system_matrix.n, &n.pinv, &mut known_vector, &mut x[..]); //maybe reuse x?
    rsparse::lsolve(&n.l, &mut x); // x = L\x
    rsparse::usolve(&n.u, &mut x[..]); // x = U\x
    rsparse::ipvec(system_matrix.n, &s.q, &x[..], &mut known_vector[..]);

    assert_eq!(known_vector,vec![2.0,3.0,4.0]);
    // ...
    // Carry on until the Simulation is finished
}

The problem is that I can't use ipvec because it is a private function. Is it possible to make ipvec public? Or is there a better solution to my problem?

[RLado]

Thank you for such a detailed issue. I think you are using the library correctly and your solution sounds very sensible. Perhaps I should consider making pvec and ipvec public for use cases such as yours.

Meanwhile, you could copy ipvec into your code and use it directly. It is a really simple permutation function:

fn ipvec(n: usize, p: &Option<Vec<isize>>, b: &[f64], x: &mut [f64]) {
    for k in 0..n {
        if p.is_some() {
            x[p.as_ref().unwrap()[k] as usize] = b[k];
        } else {
            x[k] = b[k];
        }
    }
}

[Voidheart88]

Thanks for the quick response and the easy fix. It seems i couldn't see the forest for the trees.