22//! &
33//! Methods for tridiagonal matrices
44
5- use ndarray:: * ;
65use cauchy:: Scalar ;
6+ use ndarray:: * ;
77use num_traits:: One ;
88
99use crate :: opnorm:: OperationNorm ;
@@ -34,28 +34,30 @@ pub struct TriDiagonal<A: Scalar> {
3434pub trait ToTriDiagonal < A : Scalar > {
3535 /// Extract tridiagonal elements and layout of the raw matrix.
3636 /// And also calculate 1-norm.
37- ///
37+ ///
3838 /// If the raw matrix has some non-tridiagonal elements,
3939 /// they will be ignored.
40- ///
40+ ///
4141 /// The shape of raw matrix should be equal to or larger than (2, 2).
4242 fn to_tridiagonal ( & self ) -> Result < TriDiagonal < A > > ;
4343}
4444
4545impl < A , S > ToTriDiagonal < A > for ArrayBase < S , Ix2 >
4646where
4747 A : Scalar + Lapack ,
48- S : Data < Elem = A >
48+ S : Data < Elem = A > ,
4949{
5050 fn to_tridiagonal ( & self ) -> Result < TriDiagonal < A > > {
5151 let l = self . square_layout ( ) ?;
5252 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+ }
5456 let n1 = self . opnorm_one ( ) ?;
5557
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 ( ) ;
5961 Ok ( TriDiagonal { l, n1, dl, d, du } )
6062 }
6163}
@@ -73,22 +75,22 @@ pub trait SolveTriDiagonal<A: Scalar, D: Dimension> {
7375 b : ArrayBase < S , D > ,
7476 ) -> Result < ArrayBase < S , D > > ;
7577 /// 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
7779 /// `x` is the successful result.
7880 fn solve_t_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , D > ) -> Result < Array < A , D > > ;
7981 /// 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
8183 /// `x` is the successful result.
8284 fn solve_t_tridiagonal_into < S : DataMut < Elem = A > > (
8385 & self ,
8486 b : ArrayBase < S , D > ,
8587 ) -> Result < ArrayBase < S , D > > ;
8688 /// 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
8890 /// `x` is the successful result.
8991 fn solve_h_tridiagonal < S : Data < Elem = A > > ( & self , b : & ArrayBase < S , D > ) -> Result < Array < A , D > > ;
9092 /// 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
9294 /// `x` is the successful result.
9395 fn solve_h_tridiagonal_into < S : DataMut < Elem = A > > (
9496 & self ,
@@ -155,7 +157,10 @@ where
155157 self . solve_tridiagonal_inplace ( & mut b) ?;
156158 Ok ( b)
157159 }
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 > > {
159164 let mut b = replicate ( b) ;
160165 self . solve_t_tridiagonal_inplace ( & mut b) ?;
161166 Ok ( b)
@@ -167,7 +172,10 @@ where
167172 self . solve_t_tridiagonal_inplace ( & mut b) ?;
168173 Ok ( b)
169174 }
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 > > {
171179 let mut b = replicate ( b) ;
172180 self . solve_h_tridiagonal_inplace ( & mut b) ?;
173181 Ok ( b)
@@ -186,7 +194,10 @@ where
186194 A : Scalar + Lapack ,
187195 S : Data < Elem = A > ,
188196{
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 > > {
190201 let mut b = replicate ( b) ;
191202 self . solve_tridiagonal_inplace ( & mut b) ?;
192203 Ok ( b)
@@ -198,7 +209,10 @@ where
198209 self . solve_tridiagonal_inplace ( & mut b) ?;
199210 Ok ( b)
200211 }
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 > > {
202216 let mut b = replicate ( b) ;
203217 self . solve_t_tridiagonal_inplace ( & mut b) ?;
204218 Ok ( b)
@@ -210,7 +224,10 @@ where
210224 self . solve_t_tridiagonal_inplace ( & mut b) ?;
211225 Ok ( b)
212226 }
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 > > {
214231 let mut b = replicate ( b) ;
215232 self . solve_h_tridiagonal_inplace ( & mut b) ?;
216233 Ok ( b)
@@ -334,7 +351,10 @@ where
334351 let b = self . solve_tridiagonal_into ( b) ?;
335352 Ok ( flatten ( b) )
336353 }
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 > > {
338358 let b = b. to_owned ( ) ;
339359 self . solve_t_tridiagonal_into ( b)
340360 }
@@ -346,7 +366,10 @@ where
346366 let b = self . solve_t_tridiagonal_into ( b) ?;
347367 Ok ( flatten ( b) )
348368 }
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 > > {
350373 let b = b. to_owned ( ) ;
351374 self . solve_h_tridiagonal_into ( b)
352375 }
@@ -365,7 +388,10 @@ where
365388 A : Scalar + Lapack ,
366389 S : Data < Elem = A > ,
367390{
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 > > {
369395 let b = b. to_owned ( ) ;
370396 self . solve_tridiagonal_into ( b)
371397 }
@@ -378,7 +404,10 @@ where
378404 let b = f. solve_tridiagonal_into ( b) ?;
379405 Ok ( flatten ( b) )
380406 }
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 > > {
382411 let b = b. to_owned ( ) ;
383412 self . solve_t_tridiagonal_into ( b)
384413 }
@@ -391,7 +420,10 @@ where
391420 let b = f. solve_t_tridiagonal_into ( b) ?;
392421 Ok ( flatten ( b) )
393422 }
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 > > {
395427 let b = b. to_owned ( ) ;
396428 self . solve_h_tridiagonal_into ( b)
397429 }
@@ -429,7 +461,7 @@ where
429461 Ok ( LUFactorizedTriDiagonal {
430462 a : self ,
431463 du2 : du2,
432- ipiv : ipiv
464+ ipiv : ipiv,
433465 } )
434466 }
435467}
@@ -448,7 +480,7 @@ where
448480impl < A , S > FactorizeTriDiagonal < A > for ArrayBase < S , Ix2 >
449481where
450482 A : Scalar + Lapack ,
451- S : Data < Elem = A >
483+ S : Data < Elem = A > ,
452484{
453485 fn factorize_tridiagonal ( & self ) -> Result < LUFactorizedTriDiagonal < A > > {
454486 let mut a = self . to_tridiagonal ( ) ?;
@@ -462,19 +494,19 @@ where
462494/// where {a_1, a_2, ..., a_n} are diagonal elements,
463495/// {b_1, b_2, ..., b_{n-1}} are super-diagonal elements, and
464496/// {c_1, c_2, ..., c_{n-1}} are sub-diagonal elements of matrix.
465- ///
497+ ///
466498/// f[n] is used to calculate the determinant.
467499/// (https://en.wikipedia.org/wiki/Tridiagonal_matrix#Determinant)
468- ///
500+ ///
469501/// In the future, the vector `f` can be used to calculate the inverce matrix.
470502/// (https://en.wikipedia.org/wiki/Tridiagonal_matrix#Inversion)
471503fn rec_rel < A : Scalar > ( tridiag : & TriDiagonal < A > ) -> Vec < A > {
472504 let n = tridiag. d . shape ( ) [ 0 ] ;
473- let mut f = Vec :: with_capacity ( n+ 1 ) ;
505+ let mut f = Vec :: with_capacity ( n + 1 ) ;
474506 f. push ( One :: one ( ) ) ;
475507 f. push ( tridiag. d [ 0 ] ) ;
476508 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 ] ) ;
478510 }
479511 f
480512}
@@ -483,7 +515,7 @@ fn rec_rel<A: Scalar>(tridiag: &TriDiagonal<A>) -> Vec<A> {
483515pub trait DeterminantTriDiagonal < A : Scalar > {
484516 /// Computes the determinant of the matrix.
485517 /// 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
487519 /// but the determinant itself.
488520 fn det_tridiagonal ( & self ) -> Result < A > ;
489521}
0 commit comments