del_ls/
linearsystem.rs

1//! Linear system solver class
2
3use num_traits::AsPrimitive;
4
5/// class of linear system solver
6/// * `sparse` - sparse square coefficient matrix
7/// * `r_vec` - residual vector (i.e., rhs vector)
8pub struct Solver<T> {
9    pub sparse: crate::sparse_square::Matrix<T>,
10    pub ilu: crate::sparse_ilu::Preconditioner<T>,
11    pub merge_buffer: Vec<usize>,
12    pub r_vec: Vec<T>,
13    pub u_vec: Vec<T>,
14    pub conv: Vec<T>,
15    pub conv_ratio_tol: T,
16    pub max_num_iteration: usize,
17    pub ap_vec: Vec<T>,
18    pub p_vec: Vec<T>,
19}
20
21impl<T> Solver<T>
22where
23    T: 'static
24        + Copy
25        + num_traits::Float
26        + std::default::Default
27        + std::ops::AddAssign
28        + std::fmt::Display
29        + std::ops::MulAssign
30        + std::ops::SubAssign,
31    f32: num_traits::AsPrimitive<T>,
32{
33    pub fn new() -> Self {
34        Solver {
35            sparse: crate::sparse_square::Matrix::<T>::new(),
36            ilu: crate::sparse_ilu::Preconditioner::new(),
37            merge_buffer: Vec::<usize>::new(),
38            ap_vec: Vec::<T>::new(),
39            p_vec: Vec::<T>::new(),
40            r_vec: Vec::<T>::new(),
41            u_vec: Vec::<T>::new(),
42            conv: Vec::<T>::new(),
43            conv_ratio_tol: 1.0e-5_f32.as_(),
44            max_num_iteration: 100,
45        }
46    }
47
48    pub fn initialize(&mut self, colind: &Vec<usize>, rowptr: &Vec<usize>) {
49        self.sparse.symbolic_initialization(&colind, &rowptr);
50        //self.ilu.initialize(&self.sparse);
51        let nblk = colind.len() - 1;
52        self.r_vec.resize(nblk, T::zero());
53        self.ilu.initialize_ilu0(&self.sparse);
54        //self.ilu.Initialize_ILUk(&self.sparse, 0);
55        //self.ilu.initialize_iluk(&self.sparse, 10000);
56        //self.ilu.initialize_full(self.sparse.num_blk);
57    }
58
59    /// set zero to the matrix
60    pub fn begin_merge(&mut self) {
61        self.sparse.set_zero();
62        let nblk = self.sparse.num_blk;
63        self.r_vec.resize(nblk, T::zero());
64        self.r_vec.iter_mut().for_each(|v| v.set_zero());
65    }
66
67    pub fn end_merge(&mut self) {
68        // apply boundary condition here
69        crate::sparse_ilu::copy_value(&mut self.ilu, &self.sparse);
70        crate::sparse_ilu::decompose(&mut self.ilu);
71    }
72
73    pub fn solve_cg(&mut self) {
74        self.conv = crate::solver_sparse::conjugate_gradient(
75            &mut self.r_vec,
76            &mut self.u_vec,
77            &mut self.ap_vec,
78            &mut self.p_vec,
79            self.conv_ratio_tol,
80            self.max_num_iteration,
81            &self.sparse,
82        );
83    }
84
85    pub fn solve_pcg(&mut self) {
86        self.conv = crate::solver_sparse::preconditioned_conjugate_gradient(
87            &mut self.r_vec,
88            &mut self.u_vec,
89            &mut self.ap_vec,
90            &mut self.p_vec,
91            self.conv_ratio_tol,
92            self.max_num_iteration,
93            &self.sparse,
94            &self.ilu,
95        );
96    }
97
98    pub fn clone(&self) -> Self {
99        Solver {
100            sparse: self.sparse.clone(),
101            ilu: self.ilu.clone(),
102            merge_buffer: self.merge_buffer.clone(),
103            ap_vec: self.ap_vec.clone(),
104            p_vec: self.p_vec.clone(),
105            r_vec: self.r_vec.clone(),
106            u_vec: self.u_vec.clone(),
107            conv: self.conv.clone(),
108            conv_ratio_tol: self.conv_ratio_tol,
109            max_num_iteration: 100,
110        }
111    }
112}
del_ls/linearsystem.rs

del_ls/
linearsystem.rs
