2
2
//! &
3
3
//! Methods for tridiagonal matrices
4
4
5
- use ndarray:: * ;
6
5
use cauchy:: Scalar ;
6
+ use ndarray:: * ;
7
7
use num_traits:: One ;
8
8
9
9
use crate :: opnorm:: OperationNorm ;
@@ -34,28 +34,30 @@ pub struct TriDiagonal<A: Scalar> {
34
34
pub trait ToTriDiagonal < A : Scalar > {
35
35
/// Extract tridiagonal elements and layout of the raw matrix.
36
36
/// And also calculate 1-norm.
37
- ///
37
+ ///
38
38
/// If the raw matrix has some non-tridiagonal elements,
39
39
/// they will be ignored.
40
- ///
40
+ ///
41
41
/// The shape of raw matrix should be equal to or larger than (2, 2).
42
42
fn to_tridiagonal ( & self ) -> Result < TriDiagonal < A > > ;
43
43
}
44
44
45
45
impl < A , S > ToTriDiagonal < A > for ArrayBase < S , Ix2 >
46
46
where
47
47
A : Scalar + Lapack ,
48
- S : Data < Elem = A >
48
+ S : Data < Elem = A > ,
49
49
{
50
50
fn to_tridiagonal ( & self ) -> Result < TriDiagonal < A > > {
51
51
let l = self . square_layout ( ) ?;
52
52
let ( n, _) = l. size ( ) ;
53
- if n < 2 { panic ! ( "Cannot make a tridiagonal matrix of shape=(1, 1)!" ) ; }
53
+ if n < 2 {
54
+ panic ! ( "Cannot make a tridiagonal matrix of shape=(1, 1)!" ) ;
55
+ }
54
56
let n1 = self . opnorm_one ( ) ?;
55
57
56
- let dl = self . slice ( s ! [ 1 ..n, 0 ..n- 1 ] ) . diag ( ) . to_owned ( ) ;
57
- let d = self . diag ( ) . to_owned ( ) ;
58
- let du = self . slice ( s ! [ 0 ..n- 1 , 1 ..n] ) . diag ( ) . to_owned ( ) ;
58
+ let dl = self . slice ( s ! [ 1 ..n, 0 ..n - 1 ] ) . diag ( ) . to_owned ( ) ;
59
+ let d = self . diag ( ) . to_owned ( ) ;
60
+ let du = self . slice ( s ! [ 0 ..n - 1 , 1 ..n] ) . diag ( ) . to_owned ( ) ;
59
61
Ok ( TriDiagonal { l, n1, dl, d, du } )
60
62
}
61
63
}
@@ -73,22 +75,22 @@ pub trait SolveTriDiagonal<A: Scalar, D: Dimension> {
73
75
b : ArrayBase < S , D > ,
74
76
) -> Result < ArrayBase < S , D > > ;
75
77
/// Solves a system of linear equations `A^T * x = b` with tridiagonal
76
- /// matrix `A`, where `A` is `self`, `b` is the argument, and
78
+ /// matrix `A`, where `A` is `self`, `b` is the argument, and
77
79
/// `x` is the successful result.
78
80
fn solve_t_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , D > ) -> Result < Array < A , D > > ;
79
81
/// Solves a system of linear equations `A^T * x = b` with tridiagonal
80
- /// matrix `A`, where `A` is `self`, `b` is the argument, and
82
+ /// matrix `A`, where `A` is `self`, `b` is the argument, and
81
83
/// `x` is the successful result.
82
84
fn solve_t_tridiagonal_into < S : DataMut < Elem = A > > (
83
85
& self ,
84
86
b : ArrayBase < S , D > ,
85
87
) -> Result < ArrayBase < S , D > > ;
86
88
/// Solves a system of linear equations `A^H * x = b` with tridiagonal
87
- /// matrix `A`, where `A` is `self`, `b` is the argument, and
89
+ /// matrix `A`, where `A` is `self`, `b` is the argument, and
88
90
/// `x` is the successful result.
89
91
fn solve_h_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , D > ) -> Result < Array < A , D > > ;
90
92
/// Solves a system of linear equations `A^H * x = b` with tridiagonal
91
- /// matrix `A`, where `A` is `self`, `b` is the argument, and
93
+ /// matrix `A`, where `A` is `self`, `b` is the argument, and
92
94
/// `x` is the successful result.
93
95
fn solve_h_tridiagonal_into < S : DataMut < Elem = A > > (
94
96
& self ,
@@ -155,7 +157,10 @@ where
155
157
self . solve_tridiagonal_inplace ( & mut b) ?;
156
158
Ok ( b)
157
159
}
158
- fn solve_t_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix2 > ) -> Result < Array < A , Ix2 > > {
160
+ fn solve_t_tridiagonal < S : Data < Elem = A > > (
161
+ & self ,
162
+ b : & ArrayBase < S , Ix2 > ,
163
+ ) -> Result < Array < A , Ix2 > > {
159
164
let mut b = replicate ( b) ;
160
165
self . solve_t_tridiagonal_inplace ( & mut b) ?;
161
166
Ok ( b)
@@ -167,7 +172,10 @@ where
167
172
self . solve_t_tridiagonal_inplace ( & mut b) ?;
168
173
Ok ( b)
169
174
}
170
- fn solve_h_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix2 > ) -> Result < Array < A , Ix2 > > {
175
+ fn solve_h_tridiagonal < S : Data < Elem = A > > (
176
+ & self ,
177
+ b : & ArrayBase < S , Ix2 > ,
178
+ ) -> Result < Array < A , Ix2 > > {
171
179
let mut b = replicate ( b) ;
172
180
self . solve_h_tridiagonal_inplace ( & mut b) ?;
173
181
Ok ( b)
@@ -186,7 +194,10 @@ where
186
194
A : Scalar + Lapack ,
187
195
S : Data < Elem = A > ,
188
196
{
189
- fn solve_tridiagonal < Sb : Data < Elem = A > > ( & self , b : & ArrayBase < Sb , Ix2 > ) -> Result < Array < A , Ix2 > > {
197
+ fn solve_tridiagonal < Sb : Data < Elem = A > > (
198
+ & self ,
199
+ b : & ArrayBase < Sb , Ix2 > ,
200
+ ) -> Result < Array < A , Ix2 > > {
190
201
let mut b = replicate ( b) ;
191
202
self . solve_tridiagonal_inplace ( & mut b) ?;
192
203
Ok ( b)
@@ -198,7 +209,10 @@ where
198
209
self . solve_tridiagonal_inplace ( & mut b) ?;
199
210
Ok ( b)
200
211
}
201
- fn solve_t_tridiagonal < Sb : Data < Elem = A > > ( & self , b : & ArrayBase < Sb , Ix2 > ) -> Result < Array < A , Ix2 > > {
212
+ fn solve_t_tridiagonal < Sb : Data < Elem = A > > (
213
+ & self ,
214
+ b : & ArrayBase < Sb , Ix2 > ,
215
+ ) -> Result < Array < A , Ix2 > > {
202
216
let mut b = replicate ( b) ;
203
217
self . solve_t_tridiagonal_inplace ( & mut b) ?;
204
218
Ok ( b)
@@ -210,7 +224,10 @@ where
210
224
self . solve_t_tridiagonal_inplace ( & mut b) ?;
211
225
Ok ( b)
212
226
}
213
- fn solve_h_tridiagonal < Sb : Data < Elem = A > > ( & self , b : & ArrayBase < Sb , Ix2 > ) -> Result < Array < A , Ix2 > > {
227
+ fn solve_h_tridiagonal < Sb : Data < Elem = A > > (
228
+ & self ,
229
+ b : & ArrayBase < Sb , Ix2 > ,
230
+ ) -> Result < Array < A , Ix2 > > {
214
231
let mut b = replicate ( b) ;
215
232
self . solve_h_tridiagonal_inplace ( & mut b) ?;
216
233
Ok ( b)
@@ -334,7 +351,10 @@ where
334
351
let b = self . solve_tridiagonal_into ( b) ?;
335
352
Ok ( flatten ( b) )
336
353
}
337
- fn solve_t_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix1 > ) -> Result < Array < A , Ix1 > > {
354
+ fn solve_t_tridiagonal < S : Data < Elem = A > > (
355
+ & self ,
356
+ b : & ArrayBase < S , Ix1 > ,
357
+ ) -> Result < Array < A , Ix1 > > {
338
358
let b = b. to_owned ( ) ;
339
359
self . solve_t_tridiagonal_into ( b)
340
360
}
@@ -346,7 +366,10 @@ where
346
366
let b = self . solve_t_tridiagonal_into ( b) ?;
347
367
Ok ( flatten ( b) )
348
368
}
349
- fn solve_h_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , Ix1 > ) -> Result < Array < A , Ix1 > > {
369
+ fn solve_h_tridiagonal < S : Data < Elem = A > > (
370
+ & self ,
371
+ b : & ArrayBase < S , Ix1 > ,
372
+ ) -> Result < Array < A , Ix1 > > {
350
373
let b = b. to_owned ( ) ;
351
374
self . solve_h_tridiagonal_into ( b)
352
375
}
@@ -365,7 +388,10 @@ where
365
388
A : Scalar + Lapack ,
366
389
S : Data < Elem = A > ,
367
390
{
368
- fn solve_tridiagonal < Sb : Data < Elem = A > > ( & self , b : & ArrayBase < Sb , Ix1 > ) -> Result < Array < A , Ix1 > > {
391
+ fn solve_tridiagonal < Sb : Data < Elem = A > > (
392
+ & self ,
393
+ b : & ArrayBase < Sb , Ix1 > ,
394
+ ) -> Result < Array < A , Ix1 > > {
369
395
let b = b. to_owned ( ) ;
370
396
self . solve_tridiagonal_into ( b)
371
397
}
@@ -378,7 +404,10 @@ where
378
404
let b = f. solve_tridiagonal_into ( b) ?;
379
405
Ok ( flatten ( b) )
380
406
}
381
- fn solve_t_tridiagonal < Sb : Data < Elem = A > > ( & self , b : & ArrayBase < Sb , Ix1 > ) -> Result < Array < A , Ix1 > > {
407
+ fn solve_t_tridiagonal < Sb : Data < Elem = A > > (
408
+ & self ,
409
+ b : & ArrayBase < Sb , Ix1 > ,
410
+ ) -> Result < Array < A , Ix1 > > {
382
411
let b = b. to_owned ( ) ;
383
412
self . solve_t_tridiagonal_into ( b)
384
413
}
@@ -391,7 +420,10 @@ where
391
420
let b = f. solve_t_tridiagonal_into ( b) ?;
392
421
Ok ( flatten ( b) )
393
422
}
394
- fn solve_h_tridiagonal < Sb : Data < Elem = A > > ( & self , b : & ArrayBase < Sb , Ix1 > ) -> Result < Array < A , Ix1 > > {
423
+ fn solve_h_tridiagonal < Sb : Data < Elem = A > > (
424
+ & self ,
425
+ b : & ArrayBase < Sb , Ix1 > ,
426
+ ) -> Result < Array < A , Ix1 > > {
395
427
let b = b. to_owned ( ) ;
396
428
self . solve_h_tridiagonal_into ( b)
397
429
}
@@ -429,7 +461,7 @@ where
429
461
Ok ( LUFactorizedTriDiagonal {
430
462
a : self ,
431
463
du2 : du2,
432
- ipiv : ipiv
464
+ ipiv : ipiv,
433
465
} )
434
466
}
435
467
}
@@ -448,7 +480,7 @@ where
448
480
impl < A , S > FactorizeTriDiagonal < A > for ArrayBase < S , Ix2 >
449
481
where
450
482
A : Scalar + Lapack ,
451
- S : Data < Elem = A >
483
+ S : Data < Elem = A > ,
452
484
{
453
485
fn factorize_tridiagonal ( & self ) -> Result < LUFactorizedTriDiagonal < A > > {
454
486
let mut a = self . to_tridiagonal ( ) ?;
@@ -462,19 +494,19 @@ where
462
494
/// where {a_1, a_2, ..., a_n} are diagonal elements,
463
495
/// {b_1, b_2, ..., b_{n-1}} are super-diagonal elements, and
464
496
/// {c_1, c_2, ..., c_{n-1}} are sub-diagonal elements of matrix.
465
- ///
497
+ ///
466
498
/// f[n] is used to calculate the determinant.
467
499
/// (https://en.wikipedia.org/wiki/Tridiagonal_matrix#Determinant)
468
- ///
500
+ ///
469
501
/// In the future, the vector `f` can be used to calculate the inverce matrix.
470
502
/// (https://en.wikipedia.org/wiki/Tridiagonal_matrix#Inversion)
471
503
fn rec_rel < A : Scalar > ( tridiag : & TriDiagonal < A > ) -> Vec < A > {
472
504
let n = tridiag. d . shape ( ) [ 0 ] ;
473
- let mut f = Vec :: with_capacity ( n+ 1 ) ;
505
+ let mut f = Vec :: with_capacity ( n + 1 ) ;
474
506
f. push ( One :: one ( ) ) ;
475
507
f. push ( tridiag. d [ 0 ] ) ;
476
508
for i in 1 ..n {
477
- f. push ( tridiag. d [ i] * f[ i] - tridiag. dl [ i- 1 ] * tridiag. du [ i- 1 ] * f[ i- 1 ] ) ;
509
+ f. push ( tridiag. d [ i] * f[ i] - tridiag. dl [ i - 1 ] * tridiag. du [ i - 1 ] * f[ i - 1 ] ) ;
478
510
}
479
511
f
480
512
}
@@ -483,7 +515,7 @@ fn rec_rel<A: Scalar>(tridiag: &TriDiagonal<A>) -> Vec<A> {
483
515
pub trait DeterminantTriDiagonal < A : Scalar > {
484
516
/// Computes the determinant of the matrix.
485
517
/// Unlike `.det()` of Determinant trait, this method
486
- /// doesn't returns the natural logarithm of the determinant
518
+ /// doesn't returns the natural logarithm of the determinant
487
519
/// but the determinant itself.
488
520
fn det_tridiagonal ( & self ) -> Result < A > ;
489
521
}
0 commit comments