Skip to content

Commit 261800f

Browse files
cpmechcpmech
committed
[Important] Rewrite Real Sparse Solvers such that they take a CooMatrix instead of SparseMatrix (to be removed).
1 parent 4c5c104 commit 261800f

26 files changed

+236
-510
lines changed

‎russell_ode/examples/pde_laplace_equation.rs

+3-4
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
use plotpy::{Contour, Plot};
22
use russell_lab::{StrError, Vector};
33
use russell_ode::{PdeDiscreteLaplacian2d, Side};
4-
use russell_sparse::{Genie, LinSolver, SparseMatrix};
4+
use russell_sparse::{Genie, LinSolver};
55

66
fn main() -> Result<(), StrError> {
77
// Approximate (with the Finite Differences Method, FDM) the solution of
@@ -53,10 +53,9 @@ fn main() -> Result<(), StrError> {
5353
});
5454

5555
// solve the linear system
56-
let mut mat = SparseMatrix::from_coo(aa);
5756
let mut solver = LinSolver::new(Genie::Umfpack)?;
58-
solver.actual.factorize(&mut mat, None)?;
59-
solver.actual.solve(&mut phi, &mut mat, &rhs, false)?;
57+
solver.actual.factorize(&aa, None)?;
58+
solver.actual.solve(&mut phi, &rhs, false)?;
6059

6160
// plot results
6261
let mut contour = Contour::new();

‎russell_ode/src/euler_backward.rs

+6-8
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
use crate::StrError;
22
use crate::{OdeSolverTrait, Params, System, Workspace};
33
use russell_lab::{vec_copy, vec_rms_scaled, vec_update, Vector};
4-
use russell_sparse::{numerical_jacobian, LinSolver, SparseMatrix};
4+
use russell_sparse::{numerical_jacobian, CooMatrix, LinSolver};
55

66
/// Implements the backward Euler (implicit) solver (implicit, order 1, unconditionally stable)
77
pub(crate) struct EulerBackward<'a, A> {
@@ -26,7 +26,7 @@ pub(crate) struct EulerBackward<'a, A> {
2626
dy: Vector,
2727

2828
/// Coefficient matrix K = h J - I
29-
kk: SparseMatrix,
29+
kk: CooMatrix,
3030

3131
/// Linear solver
3232
solver: LinSolver<'a>,
@@ -50,7 +50,7 @@ impl<'a, A> EulerBackward<'a, A> {
5050
w: Vector::new(ndim),
5151
r: Vector::new(ndim),
5252
dy: Vector::new(ndim),
53-
kk: SparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(),
53+
kk: CooMatrix::new(ndim, ndim, nnz, sym).unwrap(),
5454
solver: LinSolver::new(params.newton.genie).unwrap(),
5555
}
5656
}
@@ -103,7 +103,7 @@ impl<'a, A> OdeSolverTrait<A> for EulerBackward<'a, A> {
103103
work.stats.n_jacobian += 1;
104104

105105
// calculate J_new := h J
106-
let kk = self.kk.get_coo_mut().unwrap();
106+
let kk = &mut self.kk;
107107
if self.params.newton.use_numerical_jacobian || self.system.jacobian.is_none() {
108108
work.stats.n_function += ndim;
109109
let w1 = &mut self.k; // workspace
@@ -124,16 +124,14 @@ impl<'a, A> OdeSolverTrait<A> for EulerBackward<'a, A> {
124124
// perform factorization
125125
work.stats.sw_factor.reset();
126126
work.stats.n_factor += 1;
127-
self.solver
128-
.actual
129-
.factorize(&mut self.kk, self.params.newton.lin_sol_params)?;
127+
self.solver.actual.factorize(kk, self.params.newton.lin_sol_params)?;
130128
work.stats.stop_sw_factor();
131129
}
132130

133131
// solve the linear system
134132
work.stats.sw_lin_sol.reset();
135133
work.stats.n_lin_sol += 1;
136-
self.solver.actual.solve(&mut self.dy, &self.kk, &self.r, false)?;
134+
self.solver.actual.solve(&mut self.dy, &self.r, false)?;
137135
work.stats.stop_sw_lin_sol();
138136

139137
// update y

‎russell_ode/src/radau5.rs

+13-18
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ use crate::{OdeSolverTrait, Params, System, Workspace};
33
use russell_lab::math::SQRT_6;
44
use russell_lab::{complex_vec_zip, cpx, format_fortran, vec_copy, Complex64, ComplexVector, Vector};
55
use russell_sparse::{numerical_jacobian, ComplexCscMatrix, CooMatrix, CscMatrix};
6-
use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, Genie, LinSolver, SparseMatrix};
6+
use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, Genie, LinSolver};
77
use std::thread;
88

99
/// Implements the Radau5 method (Radau IIA) (implicit, order 5, embedded) for ODEs and DAEs
@@ -35,10 +35,10 @@ pub(crate) struct Radau5<'a, A> {
3535
mass: Option<CooMatrix>,
3636

3737
/// Holds the Jacobian matrix. J = df/dy
38-
jj: SparseMatrix,
38+
jj: CooMatrix,
3939

4040
/// Coefficient matrix (for real system). K_real = γ M - J
41-
kk_real: SparseMatrix,
41+
kk_real: CooMatrix,
4242

4343
/// Coefficient matrix (for real system). K_comp = (α + βi) M - J
4444
kk_comp: ComplexSparseMatrix,
@@ -146,8 +146,8 @@ impl<'a, A> Radau5<'a, A> {
146146
params,
147147
system,
148148
mass,
149-
jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, sym).unwrap(),
150-
kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(),
149+
jj: CooMatrix::new(ndim, ndim, jac_nnz, sym).unwrap(),
150+
kk_real: CooMatrix::new(ndim, ndim, nnz, sym).unwrap(),
151151
kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(),
152152
solver_real: LinSolver::new(params.newton.genie).unwrap(),
153153
solver_comp: ComplexLinSolver::new(params.newton.genie).unwrap(),
@@ -193,8 +193,8 @@ impl<'a, A> Radau5<'a, A> {
193193
/// Assembles the K_real and K_comp matrices
194194
fn assemble(&mut self, work: &mut Workspace, x: f64, y: &Vector, h: f64, args: &mut A) -> Result<(), StrError> {
195195
// auxiliary
196-
let jj = self.jj.get_coo_mut().unwrap(); // J = df/dy
197-
let kk_real = self.kk_real.get_coo_mut().unwrap(); // K_real = γ M - J
196+
let jj = &mut self.jj; // J = df/dy
197+
let kk_real = &mut self.kk_real; // K_real = γ M - J
198198
let kk_comp = self.kk_comp.get_coo_mut().unwrap(); // K_comp = (α + βi) M - J
199199

200200
// Jacobian matrix
@@ -258,7 +258,7 @@ impl<'a, A> Radau5<'a, A> {
258258
fn factorize(&mut self) -> Result<(), StrError> {
259259
self.solver_real
260260
.actual
261-
.factorize(&mut self.kk_real, self.params.newton.lin_sol_params)?;
261+
.factorize(&self.kk_real, self.params.newton.lin_sol_params)?;
262262
self.solver_comp
263263
.actual
264264
.factorize(&mut self.kk_comp, self.params.newton.lin_sol_params)
@@ -270,7 +270,7 @@ impl<'a, A> Radau5<'a, A> {
270270
let handle_real = scope.spawn(|| {
271271
self.solver_real
272272
.actual
273-
.factorize(&mut self.kk_real, self.params.newton.lin_sol_params)
273+
.factorize(&self.kk_real, self.params.newton.lin_sol_params)
274274
.unwrap();
275275
});
276276
let handle_comp = scope.spawn(|| {
@@ -295,9 +295,7 @@ impl<'a, A> Radau5<'a, A> {
295295

296296
/// Solves the real and complex linear systems
297297
fn solve_lin_sys(&mut self) -> Result<(), StrError> {
298-
self.solver_real
299-
.actual
300-
.solve(&mut self.dw0, &self.kk_real, &self.v0, false)?;
298+
self.solver_real.actual.solve(&mut self.dw0, &self.v0, false)?;
301299
self.solver_comp
302300
.actual
303301
.solve(&mut self.dw12, &self.kk_comp, &self.v12, false)?;
@@ -308,10 +306,7 @@ impl<'a, A> Radau5<'a, A> {
308306
fn solve_lin_sys_concurrently(&mut self) -> Result<(), StrError> {
309307
thread::scope(|scope| {
310308
let handle_real = scope.spawn(|| {
311-
self.solver_real
312-
.actual
313-
.solve(&mut self.dw0, &self.kk_real, &self.v0, false)
314-
.unwrap();
309+
self.solver_real.actual.solve(&mut self.dw0, &self.v0, false).unwrap();
315310
});
316311
let handle_comp = scope.spawn(|| {
317312
self.solver_comp
@@ -567,7 +562,7 @@ impl<'a, A> OdeSolverTrait<A> for Radau5<'a, A> {
567562
}
568563

569564
// err := K_real⁻¹ rhs = (γ M - J)⁻¹ rhs (HW-VII p123 Eq.(8.20))
570-
self.solver_real.actual.solve(err, &self.kk_real, rhs, false)?;
565+
self.solver_real.actual.solve(err, rhs, false)?;
571566
work.rel_error = rms_norm(err, &self.scaling);
572567

573568
// done with the error estimate
@@ -587,7 +582,7 @@ impl<'a, A> OdeSolverTrait<A> for Radau5<'a, A> {
587582
for m in 0..ndim {
588583
rhs[m] = mez[m] + fpe[m];
589584
}
590-
self.solver_real.actual.solve(err, &self.kk_real, rhs, false)?;
585+
self.solver_real.actual.solve(err, rhs, false)?;
591586
work.rel_error = rms_norm(err, &self.scaling);
592587
}
593588
Ok(())

‎russell_ode/tests/test_pde_laplace_1.rs

+3-4
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
use russell_lab::{vec_approx_eq, Vector};
22
use russell_ode::{PdeDiscreteLaplacian2d, Side};
3-
use russell_sparse::{Genie, LinSolver, SparseMatrix};
3+
use russell_sparse::{Genie, LinSolver};
44

55
#[test]
66
fn test_pde_laplace_1() {
@@ -66,10 +66,9 @@ fn test_pde_laplace_1() {
6666
});
6767

6868
// solve the linear system
69-
let mut mat = SparseMatrix::from_coo(aa);
7069
let mut solver = LinSolver::new(Genie::Umfpack).unwrap();
71-
solver.actual.factorize(&mut mat, None).unwrap();
72-
solver.actual.solve(&mut x, &mut mat, &b, false).unwrap();
70+
solver.actual.factorize(&aa, None).unwrap();
71+
solver.actual.solve(&mut x, &b, false).unwrap();
7372

7473
// check
7574
let x_correct = [

‎russell_ode/tests/test_pde_poisson_1.rs

+3-4
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
use plotpy::{Contour, Plot};
22
use russell_lab::{vec_approx_eq, Vector};
33
use russell_ode::PdeDiscreteLaplacian2d;
4-
use russell_sparse::{Genie, LinSolver, SparseMatrix};
4+
use russell_sparse::{Genie, LinSolver};
55

66
const SAVE_FIGURE: bool = false;
77

@@ -51,10 +51,9 @@ fn test_pde_poisson_1() {
5151
});
5252

5353
// solve the linear system
54-
let mut mat = SparseMatrix::from_coo(aa);
5554
let mut solver = LinSolver::new(Genie::Umfpack).unwrap();
56-
solver.actual.factorize(&mut mat, None).unwrap();
57-
solver.actual.solve(&mut phi, &mut mat, &rhs, false).unwrap();
55+
solver.actual.factorize(&aa, None).unwrap();
56+
solver.actual.solve(&mut phi, &rhs, false).unwrap();
5857

5958
// check
6059
let mut phi_correct = Vector::new(dim);

‎russell_ode/tests/test_pde_poisson_2.rs

+3-4
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
use plotpy::{Contour, Plot};
22
use russell_lab::{math::PI, vec_approx_eq, Vector};
33
use russell_ode::{PdeDiscreteLaplacian2d, Side};
4-
use russell_sparse::{Genie, LinSolver, SparseMatrix};
4+
use russell_sparse::{Genie, LinSolver};
55

66
const SAVE_FIGURE: bool = false;
77

@@ -77,10 +77,9 @@ fn test_pde_poisson_2() {
7777
});
7878

7979
// solve the linear system
80-
let mut mat = SparseMatrix::from_coo(aa);
8180
let mut solver = LinSolver::new(Genie::Umfpack).unwrap();
82-
solver.actual.factorize(&mut mat, None).unwrap();
83-
solver.actual.solve(&mut phi, &mut mat, &rhs, false).unwrap();
81+
solver.actual.factorize(&aa, None).unwrap();
82+
solver.actual.solve(&mut phi, &rhs, false).unwrap();
8483

8584
// check
8685
let mut phi_correct = Vector::new(dim);

‎russell_ode/tests/test_pde_poisson_3.rs

+3-4
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
use plotpy::{Contour, Plot};
22
use russell_lab::{vec_approx_eq, Vector};
33
use russell_ode::PdeDiscreteLaplacian2d;
4-
use russell_sparse::{Genie, LinSolver, SparseMatrix};
4+
use russell_sparse::{Genie, LinSolver};
55

66
const SAVE_FIGURE: bool = false;
77

@@ -59,10 +59,9 @@ fn test_pde_poisson_3() {
5959
});
6060

6161
// solve the linear system
62-
let mut mat = SparseMatrix::from_coo(aa);
6362
let mut solver = LinSolver::new(Genie::Umfpack).unwrap();
64-
solver.actual.factorize(&mut mat, None).unwrap();
65-
solver.actual.solve(&mut phi, &mut mat, &rhs, false).unwrap();
63+
solver.actual.factorize(&aa, None).unwrap();
64+
solver.actual.solve(&mut phi, &rhs, false).unwrap();
6665

6766
// check
6867
let mut phi_correct = Vector::new(dim);

‎russell_sparse/README.md

+3-3
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,7 @@ fn main() -> Result<(), StrError> {
9696
let mut umfpack = SolverUMFPACK::new()?;
9797

9898
// allocate the coefficient matrix
99-
let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?;
99+
let mut coo = CooMatrix::new(ndim, ndim, nnz, Sym::No)?;
100100
coo.put(0, 0, 0.2)?;
101101
coo.put(0, 1, 0.2)?;
102102
coo.put(1, 0, 0.5)?;
@@ -113,14 +113,14 @@ fn main() -> Result<(), StrError> {
113113
assert_eq!(format!("{}", a), correct);
114114

115115
// call factorize
116-
umfpack.factorize(&mut coo, None)?;
116+
umfpack.factorize(&coo, None)?;
117117

118118
// allocate two right-hand side vectors
119119
let b = Vector::from(&[1.0, 1.0, 1.0]);
120120

121121
// calculate the solution
122122
let mut x = Vector::new(ndim);
123-
umfpack.solve(&mut x, &coo, &b, false)?;
123+
umfpack.solve(&mut x, &b, false)?;
124124
let correct = vec![3.0, 2.0, 4.0];
125125
vec_approx_eq(&x, &correct, 1e-14);
126126
Ok(())

‎russell_sparse/examples/doc_lin_solver_compute.rs

+2-2
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ fn main() -> Result<(), StrError> {
88
let nnz = 5; // number of non-zero values
99

1010
// allocate the coefficient matrix
11-
let mut mat = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?;
11+
let mut mat = CooMatrix::new(ndim, ndim, nnz, Sym::No)?;
1212
mat.put(0, 0, 0.2)?;
1313
mat.put(0, 1, 0.2)?;
1414
mat.put(1, 0, 0.5)?;
@@ -29,7 +29,7 @@ fn main() -> Result<(), StrError> {
2929

3030
// calculate the solution
3131
let mut x = Vector::new(ndim);
32-
LinSolver::compute(Genie::Umfpack, &mut x, &mut mat, &rhs, None)?;
32+
LinSolver::compute(Genie::Umfpack, &mut x, &mat, &rhs, None)?;
3333
let correct = vec![3.0, 2.0, 4.0];
3434
vec_approx_eq(&x, &correct, 1e-14);
3535
Ok(())

‎russell_sparse/examples/doc_lin_solver_umfpack_tiny.rs

+4-4
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ fn main() -> Result<(), StrError> {
1111
let mut solver = LinSolver::new(Genie::Umfpack)?;
1212

1313
// allocate the coefficient matrix
14-
let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?;
14+
let mut coo = CooMatrix::new(ndim, ndim, nnz, Sym::No)?;
1515
coo.put(0, 0, 0.2)?;
1616
coo.put(0, 1, 0.2)?;
1717
coo.put(1, 0, 0.5)?;
@@ -28,21 +28,21 @@ fn main() -> Result<(), StrError> {
2828
assert_eq!(format!("{}", a), correct);
2929

3030
// call factorize
31-
solver.actual.factorize(&mut coo, None)?;
31+
solver.actual.factorize(&coo, None)?;
3232

3333
// allocate two right-hand side vectors
3434
let rhs1 = Vector::from(&[1.0, 1.0, 1.0]);
3535
let rhs2 = Vector::from(&[2.0, 2.0, 2.0]);
3636

3737
// calculate the solution
3838
let mut x1 = Vector::new(ndim);
39-
solver.actual.solve(&mut x1, &coo, &rhs1, false)?;
39+
solver.actual.solve(&mut x1, &rhs1, false)?;
4040
let correct = vec![3.0, 2.0, 4.0];
4141
vec_approx_eq(&x1, &correct, 1e-14);
4242

4343
// calculate the solution again
4444
let mut x2 = Vector::new(ndim);
45-
solver.actual.solve(&mut x2, &coo, &rhs2, false)?;
45+
solver.actual.solve(&mut x2, &rhs2, false)?;
4646
let correct = vec![6.0, 4.0, 8.0];
4747
vec_approx_eq(&x2, &correct, 1e-14);
4848
Ok(())

‎russell_sparse/examples/doc_umfpack_quickstart_coo.rs

+3-3
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ fn main() -> Result<(), StrError> {
1616
// . -1 -3 2 .
1717
// . . 1 . .
1818
// . 4 2 . 1
19-
let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?;
19+
let mut coo = CooMatrix::new(ndim, ndim, nnz, Sym::No)?;
2020
coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
2121
coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
2222
coo.put(1, 0, 3.0)?;
@@ -37,14 +37,14 @@ fn main() -> Result<(), StrError> {
3737
params.compute_determinant = true;
3838

3939
// call factorize
40-
umfpack.factorize(&mut coo, Some(params))?;
40+
umfpack.factorize(&coo, Some(params))?;
4141

4242
// allocate x and rhs
4343
let mut x = Vector::new(ndim);
4444
let rhs = Vector::from(&[8.0, 45.0, -3.0, 3.0, 19.0]);
4545

4646
// calculate the solution
47-
umfpack.solve(&mut x, &coo, &rhs, false)?;
47+
umfpack.solve(&mut x, &rhs, false)?;
4848
println!("x =\n{}", x);
4949

5050
// check the results

0 commit comments

Comments
 (0)