1
+ extern crate ndarray;
1
2
///! This file demonstrates basic usage of the sprs library,
2
3
///! where a heat diffusion problem with Dirichlet boundary condition
3
4
///! is solved for equilibrium. For simplicity we omit any relevant
14
15
///! This shows how a laplacian matrix can be constructed by directly
15
16
///! constructing the compressed structure, and how the resulting linear
16
17
///! system can be solved using an iterative method.
17
-
18
18
extern crate sprs;
19
- extern crate ndarray;
20
19
21
20
type VecView < ' a , T > = ndarray:: ArrayView < ' a , T , ndarray:: Ix1 > ;
22
21
type VecViewMut < ' a , T > = ndarray:: ArrayViewMut < ' a , T , ndarray:: Ix1 > ;
@@ -49,7 +48,7 @@ fn grid_laplacian(shape: (usize, usize)) -> sprs::CsMat<f64> {
49
48
let ( rows, cols) = shape;
50
49
let nb_vert = rows * cols;
51
50
let mut indptr = Vec :: with_capacity ( nb_vert + 1 ) ;
52
- let nnz = 5 * nb_vert + 5 ;
51
+ let nnz = 5 * nb_vert + 5 ;
53
52
let mut indices = Vec :: with_capacity ( nnz) ;
54
53
let mut data = Vec :: with_capacity ( nnz) ;
55
54
let mut cumsum = 0 ;
@@ -67,8 +66,7 @@ fn grid_laplacian(shape: (usize, usize)) -> sprs::CsMat<f64> {
67
66
if is_border ( i, j, shape) {
68
67
// establish Dirichlet boundary conditions
69
68
add_elt ( i, j, 1. ) ;
70
- }
71
- else {
69
+ } else {
72
70
add_elt ( i - 1 , j, 1. ) ;
73
71
add_elt ( i, j - 1 , 1. ) ;
74
72
add_elt ( i, j, -4. ) ;
@@ -84,17 +82,18 @@ fn grid_laplacian(shape: (usize, usize)) -> sprs::CsMat<f64> {
84
82
}
85
83
86
84
/// Set a dirichlet boundary condition
87
- fn set_boundary_condition < F > ( mut rhs : VecViewMut < f64 > ,
88
- grid_shape : ( usize , usize ) ,
89
- f : F
90
- )
91
- where F : Fn ( usize , usize ) -> f64
85
+ fn set_boundary_condition < F > (
86
+ mut rhs : VecViewMut < f64 > ,
87
+ grid_shape : ( usize , usize ) ,
88
+ f : F ,
89
+ ) where
90
+ F : Fn ( usize , usize ) -> f64 ,
92
91
{
93
92
let ( rows, cols) = grid_shape;
94
93
for i in 0 ..rows {
95
94
for j in 0 ..cols {
96
95
if is_border ( i, j, grid_shape) {
97
- let index = i* rows + j;
96
+ let index = i * rows + j;
98
97
rhs[ [ index] ] = f ( i, j) ;
99
98
}
100
99
}
@@ -103,12 +102,13 @@ where F: Fn(usize, usize) -> f64
103
102
104
103
/// Gauss-Seidel method to solve the system
105
104
/// see https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method#Algorithm
106
- fn gauss_seidel ( mat : sprs:: CsMatView < f64 > ,
107
- mut x : VecViewMut < f64 > ,
108
- rhs : VecView < f64 > ,
109
- max_iter : usize ,
110
- eps : f64
111
- ) -> Result < ( usize , f64 ) , f64 > {
105
+ fn gauss_seidel (
106
+ mat : sprs:: CsMatView < f64 > ,
107
+ mut x : VecViewMut < f64 > ,
108
+ rhs : VecView < f64 > ,
109
+ max_iter : usize ,
110
+ eps : f64 ,
111
+ ) -> Result < ( usize , f64 ) , f64 > {
112
112
assert ! ( mat. rows( ) == mat. cols( ) ) ;
113
113
assert ! ( mat. rows( ) == x. shape( ) [ 0 ] ) ;
114
114
let mut error = ( & mat * & x - rhs) . scalar_sum ( ) . sqrt ( ) ;
@@ -119,16 +119,15 @@ fn gauss_seidel(mat: sprs::CsMatView<f64>,
119
119
for ( col_ind, & val) in vec. iter ( ) {
120
120
if row_ind != col_ind {
121
121
sigma += val * x[ [ col_ind] ] ;
122
- }
123
- else {
122
+ } else {
124
123
diag = Some ( val) ;
125
124
}
126
125
}
127
126
// Gauss-Seidel requires a non-zero diagonal, which
128
127
// is satisfied for a laplacian matrix
129
128
let diag = diag. unwrap ( ) ;
130
129
let cur_rhs = rhs[ [ row_ind] ] ;
131
- x[ [ row_ind] ] = ( cur_rhs - sigma) / diag;
130
+ x[ [ row_ind] ] = ( cur_rhs - sigma) / diag;
132
131
}
133
132
134
133
error = ( & mat * & x - rhs) . scalar_sum ( ) . sqrt ( ) ;
@@ -153,9 +152,10 @@ fn main() {
153
152
154
153
match gauss_seidel ( lap. view ( ) , x. view_mut ( ) , rhs. view ( ) , 300 , 1e-8 ) {
155
154
Ok ( ( iters, error) ) => {
156
- println ! ( "Solved system in {} iterations with residual error {}" ,
157
- iters,
158
- error) ;
155
+ println ! (
156
+ "Solved system in {} iterations with residual error {}" ,
157
+ iters, error
158
+ ) ;
159
159
let grid = x. view ( ) . into_shape ( ( rows, cols) ) . unwrap ( ) ;
160
160
for i in 0 ..rows {
161
161
for j in 0 ..cols {
0 commit comments