sparse_linear_system/
sparse_linear_system.rs

1// cargo run --example sparse_linear_system
2
3extern crate ohsl;
4
5pub use ohsl::vector::{Vector, Vec64};
6pub use ohsl::matrix::{Matrix, Mat64};
7use ohsl::sparse::Sparse;
8
9fn main() {
10    println!( "------------------ Sparse linear system ------------------" );
11    println!( "  * Solving the linear system Ax = b, for x, where" );
12    let n: usize = 4;
13    let mut a = Mat64::new( n, n, 0.0 );
14    let mut triplets = Vec::new();
15    for i in 0..n {
16        a[(i,i)] = 2.0;
17        triplets.push( ( i, i, 2.0 ) );
18        if i > 0 {
19            a[( i - 1 , i )] = 1.0;
20            triplets.push( ( i - 1, i, 1.0 ) );
21        }
22        if i < n - 1 {
23            a[( i + 1 , i )] = 1.0;
24            triplets.push( ( i + 1, i, 1.0 ) );
25        }
26    }
27    println!( "  * A ={}", a );
28    let b = Vec64::random( n );
29    println!( "  * b^T ={}", b );
30    println!( "  * as both dense and sparse linear systems.");
31    let dense_x = a.solve_basic( &b );
32    println!( "  * The dense system gives the solution vector");
33    println!( "  * x^T ={}", dense_x );
34    let sparse = Sparse::<f64>::from_triplets( n, n, &mut triplets );
35    //let sparse_x = sparse.solve( b ).unwrap();
36    let mut sparse_x = Vector::<f64>::new( n, 0.0 );
37    let max_iter = 1000;
38    let tol = 1e-8;
39    let result = sparse.solve_bicgstab( &b, &mut sparse_x, max_iter, tol );
40    println!( "  * The sparse system gives the solution vector");
41    println!( "  * x^T ={}", sparse_x );
42    match result {
43        Ok( iter ) => println!( "  * The sparse system converged in {} iterations.", iter ),
44        Err( error ) => println!( "  * The sparse system failed to converge, error = {}", error )
45    }
46    let diff = dense_x - sparse_x;
47    println!( "  * The maximum difference between the two is {:+.2e}", diff.norm_inf() );
48    println!( "-------------------- FINISHED ---------------------" );
49}
sparse_linear_system/sparse_linear_system.rs

sparse_linear_system/
sparse_linear_system.rs
