Struct Sparse 

pub struct Sparse<T> {
    pub rows: usize,
    pub cols: usize,
    pub nonzero: usize,
    pub val: Vec<T>,
    pub row_index: Vec<usize>,
    pub col_start: Vec<usize>,
}

Fields

rows: usize

cols: usize

nonzero: usize

val: Vec<T>

row_index: Vec<usize>

col_start: Vec<usize>

Implementations

impl<T: Copy + Number + Debug> Sparse<T>

pub fn from_vecs( rows: usize, cols: usize, val: Vec<T>, row_index: Vec<usize>, col_start: Vec<usize>, ) -> Self

Create a new sparse matrix of specified size using a vector of values

pub fn from_triplets( rows: usize, cols: usize, triplets: &mut Vec<(usize, usize, T)>, ) -> Self

Create a new sparse matrix of specified size using a vector of triplets

Examples found in repository

examples/sparse_linear_system.rs (line 34)

fn main() {
    println!( "------------------ Sparse linear system ------------------" );
    println!( "  * Solving the linear system Ax = b, for x, where" );
    let n: usize = 4;
    let mut a = Mat64::new( n, n, 0.0 );
    let mut triplets = Vec::new();
    for i in 0..n {
        a[(i,i)] = 2.0;
        triplets.push( ( i, i, 2.0 ) );
        if i > 0 {
            a[( i - 1 , i )] = 1.0;
            triplets.push( ( i - 1, i, 1.0 ) );
        }
        if i < n - 1 {
            a[( i + 1 , i )] = 1.0;
            triplets.push( ( i + 1, i, 1.0 ) );
        }
    }
    println!( "  * A ={}", a );
    let b = Vec64::random( n );
    println!( "  * b^T ={}", b );
    println!( "  * as both dense and sparse linear systems.");
    let dense_x = a.solve_basic( &b );
    println!( "  * The dense system gives the solution vector");
    println!( "  * x^T ={}", dense_x );
    let sparse = Sparse::<f64>::from_triplets( n, n, &mut triplets );
    //let sparse_x = sparse.solve( b ).unwrap();
    let mut sparse_x = Vector::<f64>::new( n, 0.0 );
    let max_iter = 1000;
    let tol = 1e-8;
    let result = sparse.solve_bicgstab( &b, &mut sparse_x, max_iter, tol );
    println!( "  * The sparse system gives the solution vector");
    println!( "  * x^T ={}", sparse_x );
    match result {
        Ok( iter ) => println!( "  * The sparse system converged in {} iterations.", iter ),
        Err( error ) => println!( "  * The sparse system failed to converge, error = {}", error )
    }
    let diff = dense_x - sparse_x;
    println!( "  * The maximum difference between the two is {:+.2e}", diff.norm_inf() );
    println!( "-------------------- FINISHED ---------------------" );
}

examples/sparse_laplace.rs (line 79)

fn main() {
    println!( "------------------ Laplace's equation ------------------" );
    println!( "  * Solve Laplace's equation on the unit square subject");
    println!( "  * to the given boundary conditions. The exact solution");
    println!( "  * is u(x,y) = y/[(1+x)^2 + y^2].");
    let start = Instant::now();
    let n: usize = 64;                          // Number of intervals
    let dx: f64 = 1.0 / ( n as f64 );           // Step size
    let size: usize = ( n + 1 ) * ( n + 1 );    // Size of the linear system
    let mut rhs = Vec64::zeros( size );         // Right hand side vector 

    let mut triplets = Vec::new();
    let mut row: usize = 0;
    
    // x = 0 boundary ( u = y / ( 1 + y^2 ) )
    let i = 0;
    for j in 0..n+1 {
        let y = (j as f64) * dx;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = y / ( 1.0 + y * y );
        row += 1;
    }

    for i in 1..n {
        let x = (i as f64) * dx;
        // y = 0 boundary ( u = 0 )
        let j = 0;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = 0.0;
        row += 1;
        // Interior nodes
        for j in 1..n {
            triplets.push( ( row, ( i - 1 ) * ( n + 1 ) + j, 1.0 ) );
            triplets.push( ( row, ( i + 1 ) * ( n + 1 ) + j, 1.0 ) );
            triplets.push( ( row, i * ( n + 1 ) + j - 1, 1.0 ) );
            triplets.push( ( row, i * ( n + 1 ) + j + 1, 1.0 ) );
            triplets.push( ( row, i * ( n + 1 ) + j, - 4.0 ) );
            rhs[ row ] = 0.0;
            row += 1;
        }
        // y = 1 boundary ( u = 1 / ( ( 1 + x )^2 + 1 ) )
        let j = n;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = 1. / ( ( 1. + x ) * ( 1. + x ) + 1. );
        row += 1;
    }

    // x = 1 boundary ( u = y / ( 4 + y^2 ) )
    let i = n;
    for j in 0..n+1 {
        let y = (j as f64) * dx;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = y / ( 4. + y * y );
        row += 1;
    }

    // Exact solution u = y / ( ( 1 + x )^2 + y^2 )
    let mut u_exact = Vec64::zeros( size );    
    row = 0;
    for i in 0..n+1 {
        let x = (i as f64) * dx;
        for j in 0..n+1 {
            let y = (j as f64) * dx;
            u_exact[ row ] = y / ( ( 1. + x ) * ( 1. + x ) + y * y );
            row += 1;
        }
    }
    // Create the sparse matrix from the triplets
    let sparse = Sparse::<f64>::from_triplets( size, size, &mut triplets );
    let duration = start.elapsed();
    println!("  * Time to create the matrix: {:?}", duration);
    // Solve the sparse system
    let mut u_guess = Vector::<f64>::random( size ); // A better guess would improve convergence
    let max_iter = 1000;
    let tol = 1e-8;
    let result = sparse.solve_qmr( &rhs, &mut u_guess, max_iter, tol );
    // Output time and error
    let duration = start.elapsed();
    match result {
        Ok( iter ) => println!( "  * The sparse system converged in {} iterations.", iter ),
        Err( error ) => println!( "  * The sparse system failed to converge, error = {}", error )
    }
    println!("  * Time to solve the system: {:?}", duration);
    let u_diff = u_guess - u_exact;
    println!("  * Solution error = {}", u_diff.norm_2() );
    println!( "-------------------- FINISHED ---------------------" );
}

pub fn col_index(&self) -> Vector<usize>

Return a vector of column indices relating to each value ( triplet form )

pub fn get(&self, row: usize, col: usize) -> Option<T>

Get a specific element from the sparse matrix if it exists

pub fn col_start_from_index(&self, col_index: &Vector<usize>) -> Vec<usize>

Calculate column start vector for a given column index vector

pub fn scale(&mut self, value: &T)

Insert a non-zero element into the sparse matrix at the specified row and column Scale the non-zero elements of sparse matrix by a given value

pub fn multiply(&self, x: &Vector<T>) -> Vector<T>

Multiply the sparse matrix A by a Vector x to the right and return the result vector A * x

pub fn transpose_multiply(&self, x: &Vector<T>) -> Vector<T>

Multiply the transpose matrix A^T by a Vector x and return the result vector A^T * x

pub fn transpose(&self) -> Sparse<T>

Return the transpose of the sparse matrix

pub fn to_triplets(&self) -> Vec<(usize, usize, T)>

Return a vector of triplets representing the sparse matrix

pub fn insert(&mut self, row: usize, col: usize, value: T)

Insert a non-zero element into the sparse matrix at the specified row and column If the element already exists, the value is updated. This is costly and should be avoided.

pub fn to_dense(&self) -> Matrix<T>

Return a dense matrix representation of the sparse matrix

impl Sparse<f64>

pub fn solve_bicg( &self, b: &Vector<f64>, x: &mut Vector<f64>, max_iter: usize, tol: f64, itol: usize, ) -> Result<usize, f64>

Solve the system of equations Ax=b using the biconjugate gradient method with a specified maximum number of iterations and tolerance. itol = 1: relative residual norm itol = 2: relative error norm Returns the number of iterations if successful or the error if the method fails

pub fn solve_bicgstab( &self, b: &Vector<f64>, x: &mut Vector<f64>, max_iter: usize, tol: f64, ) -> Result<usize, f64>

Solve the system of equations Ax=b using the stabilised biconjugate gradient method BiCGSTAB with a specified maximum number of iterations and tolerance. Returns the number of iterations if successful or the error if the method fails

Examples found in repository

examples/sparse_linear_system.rs (line 39)

fn main() {
    println!( "------------------ Sparse linear system ------------------" );
    println!( "  * Solving the linear system Ax = b, for x, where" );
    let n: usize = 4;
    let mut a = Mat64::new( n, n, 0.0 );
    let mut triplets = Vec::new();
    for i in 0..n {
        a[(i,i)] = 2.0;
        triplets.push( ( i, i, 2.0 ) );
        if i > 0 {
            a[( i - 1 , i )] = 1.0;
            triplets.push( ( i - 1, i, 1.0 ) );
        }
        if i < n - 1 {
            a[( i + 1 , i )] = 1.0;
            triplets.push( ( i + 1, i, 1.0 ) );
        }
    }
    println!( "  * A ={}", a );
    let b = Vec64::random( n );
    println!( "  * b^T ={}", b );
    println!( "  * as both dense and sparse linear systems.");
    let dense_x = a.solve_basic( &b );
    println!( "  * The dense system gives the solution vector");
    println!( "  * x^T ={}", dense_x );
    let sparse = Sparse::<f64>::from_triplets( n, n, &mut triplets );
    //let sparse_x = sparse.solve( b ).unwrap();
    let mut sparse_x = Vector::<f64>::new( n, 0.0 );
    let max_iter = 1000;
    let tol = 1e-8;
    let result = sparse.solve_bicgstab( &b, &mut sparse_x, max_iter, tol );
    println!( "  * The sparse system gives the solution vector");
    println!( "  * x^T ={}", sparse_x );
    match result {
        Ok( iter ) => println!( "  * The sparse system converged in {} iterations.", iter ),
        Err( error ) => println!( "  * The sparse system failed to converge, error = {}", error )
    }
    let diff = dense_x - sparse_x;
    println!( "  * The maximum difference between the two is {:+.2e}", diff.norm_inf() );
    println!( "-------------------- FINISHED ---------------------" );
}

pub fn solve_cg( &self, b: &Vector<f64>, x: &mut Vector<f64>, max_iter: usize, tol: f64, ) -> Result<usize, f64>

Solve the system of equations Ax=b using the conjugate gradient method with a specified maximum number of iterations and tolerance. Returns the number of iterations if successful or the error if the method fails

pub fn solve_qmr( &self, b: &Vector<f64>, x: &mut Vector<f64>, max_iter: usize, tol: f64, ) -> Result<usize, f64>

Solve the system of equations Ax=b using the Quasi-minimal residual method with a specified maximum number of iterations and tolerance. Returns the number of iterations if successful or the error if the method fails

Examples found in repository

examples/sparse_laplace.rs (line 86)

fn main() {
    println!( "------------------ Laplace's equation ------------------" );
    println!( "  * Solve Laplace's equation on the unit square subject");
    println!( "  * to the given boundary conditions. The exact solution");
    println!( "  * is u(x,y) = y/[(1+x)^2 + y^2].");
    let start = Instant::now();
    let n: usize = 64;                          // Number of intervals
    let dx: f64 = 1.0 / ( n as f64 );           // Step size
    let size: usize = ( n + 1 ) * ( n + 1 );    // Size of the linear system
    let mut rhs = Vec64::zeros( size );         // Right hand side vector 

    let mut triplets = Vec::new();
    let mut row: usize = 0;
    
    // x = 0 boundary ( u = y / ( 1 + y^2 ) )
    let i = 0;
    for j in 0..n+1 {
        let y = (j as f64) * dx;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = y / ( 1.0 + y * y );
        row += 1;
    }

    for i in 1..n {
        let x = (i as f64) * dx;
        // y = 0 boundary ( u = 0 )
        let j = 0;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = 0.0;
        row += 1;
        // Interior nodes
        for j in 1..n {
            triplets.push( ( row, ( i - 1 ) * ( n + 1 ) + j, 1.0 ) );
            triplets.push( ( row, ( i + 1 ) * ( n + 1 ) + j, 1.0 ) );
            triplets.push( ( row, i * ( n + 1 ) + j - 1, 1.0 ) );
            triplets.push( ( row, i * ( n + 1 ) + j + 1, 1.0 ) );
            triplets.push( ( row, i * ( n + 1 ) + j, - 4.0 ) );
            rhs[ row ] = 0.0;
            row += 1;
        }
        // y = 1 boundary ( u = 1 / ( ( 1 + x )^2 + 1 ) )
        let j = n;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = 1. / ( ( 1. + x ) * ( 1. + x ) + 1. );
        row += 1;
    }

    // x = 1 boundary ( u = y / ( 4 + y^2 ) )
    let i = n;
    for j in 0..n+1 {
        let y = (j as f64) * dx;
        triplets.push( ( row, i * ( n + 1 ) + j, 1.0 ) );
        rhs[ row ] = y / ( 4. + y * y );
        row += 1;
    }

    // Exact solution u = y / ( ( 1 + x )^2 + y^2 )
    let mut u_exact = Vec64::zeros( size );    
    row = 0;
    for i in 0..n+1 {
        let x = (i as f64) * dx;
        for j in 0..n+1 {
            let y = (j as f64) * dx;
            u_exact[ row ] = y / ( ( 1. + x ) * ( 1. + x ) + y * y );
            row += 1;
        }
    }
    // Create the sparse matrix from the triplets
    let sparse = Sparse::<f64>::from_triplets( size, size, &mut triplets );
    let duration = start.elapsed();
    println!("  * Time to create the matrix: {:?}", duration);
    // Solve the sparse system
    let mut u_guess = Vector::<f64>::random( size ); // A better guess would improve convergence
    let max_iter = 1000;
    let tol = 1e-8;
    let result = sparse.solve_qmr( &rhs, &mut u_guess, max_iter, tol );
    // Output time and error
    let duration = start.elapsed();
    match result {
        Ok( iter ) => println!( "  * The sparse system converged in {} iterations.", iter ),
        Err( error ) => println!( "  * The sparse system failed to converge, error = {}", error )
    }
    println!("  * Time to solve the system: {:?}", duration);
    let u_diff = u_guess - u_exact;
    println!("  * Solution error = {}", u_diff.norm_2() );
    println!( "-------------------- FINISHED ---------------------" );
}