68
68
69
69
impl < A , S , F , Ortho > Iterator for Arnoldi < A , S , F , Ortho >
70
70
where
71
- A : Scalar ,
71
+ A : Scalar + Lapack ,
72
72
S : DataMut < Elem = A > ,
73
73
F : Fn ( & mut ArrayBase < S , Ix1 > ) ,
74
74
Ortho : Orthogonalizer < Elem = A > ,
@@ -78,15 +78,14 @@ where
78
78
fn next ( & mut self ) -> Option < Self :: Item > {
79
79
( self . a ) ( & mut self . v ) ;
80
80
let result = self . ortho . div_append ( & mut self . v ) ;
81
- azip ! ( mut v( & mut self . v) in { * v = v. div_real( result. residual_norm( ) ) } ) ;
81
+ let norm = self . v . norm_l2 ( ) ;
82
+ azip ! ( mut v( & mut self . v) in { * v = v. div_real( norm) } ) ;
82
83
match result {
83
84
AppendResult :: Added ( coef) => {
84
- dbg ! ( & coef) ;
85
85
self . h . push ( coef. clone ( ) ) ;
86
86
Some ( coef)
87
87
}
88
88
AppendResult :: Dependent ( coef) => {
89
- dbg ! ( & coef) ;
90
89
self . h . push ( coef) ;
91
90
None
92
91
}
@@ -135,7 +134,7 @@ where
135
134
#[ cfg( test) ]
136
135
mod tests {
137
136
use super :: * ;
138
- use crate :: generate:: * ;
137
+ use crate :: { assert :: * , generate:: * } ;
139
138
140
139
#[ test]
141
140
fn aq_qh ( ) {
@@ -146,23 +145,26 @@ mod tests {
146
145
println ! ( "A = \n {:?}" , & a) ;
147
146
println ! ( "Q = \n {:?}" , & q) ;
148
147
println ! ( "H = \n {:?}" , & h) ;
148
+ let aq = a. dot ( & q) ;
149
+ let qh = q. dot ( & h) ;
149
150
println ! ( "AQ = \n {:?}" , a. dot( & q) ) ;
150
151
println ! ( "QH = \n {:?}" , q. dot( & h) ) ;
151
- panic ! ( )
152
+ close_l2 ( & aq , & qh , 1e-9 ) ;
152
153
}
153
154
154
155
#[ test]
155
156
fn aq_qh_random ( ) {
156
157
let a: Array2 < f64 > = random ( ( 5 , 5 ) ) ;
157
- let mut v = Array :: zeros ( 5 ) ;
158
- v[ 0 ] = 1.0 ;
158
+ let v: Array1 < f64 > = random ( 5 ) ;
159
159
let ( q, h) = arnoldi_mgs ( a. clone ( ) , v, 1e-9 ) ;
160
160
println ! ( "A = \n {:?}" , & a) ;
161
161
println ! ( "Q = \n {:?}" , & q) ;
162
162
println ! ( "H = \n {:?}" , & h) ;
163
- println ! ( "AQ = \n {:?}" , a. dot( & q) ) ;
164
- println ! ( "QH = \n {:?}" , q. dot( & h) ) ;
165
- panic ! ( )
163
+ let aq = a. dot ( & q) ;
164
+ let qh = q. dot ( & h) ;
165
+ println ! ( "AQ = \n {:?}" , & aq) ;
166
+ println ! ( "QH = \n {:?}" , & qh) ;
167
+ close_l2 ( & aq, & qh, 1e-9 ) ;
166
168
}
167
169
168
170
}
0 commit comments