@@ -2,8 +2,16 @@ using IterativeSolvers
22using LinearMaps
33using Base. Test
44
5+ import Base. A_ldiv_B!
6+
57include (" laplace_matrix.jl" )
68
9+ struct JacobiPrec{TD}
10+ diagonal:: TD
11+ end
12+
13+ A_ldiv_B! (y, P:: JacobiPrec , x) = y .= x ./ P. diagonal
14+
715function max_err (R)
816 r = zeros (real (eltype (R)), size (R, 2 ))
917 for j in 1 : length (r)
1826@testset " Locally Optimal Block Preconditioned Conjugate Gradient" begin
1927 srand (1234321 )
2028 @testset " Single eigenvalue" begin
29+ n = 10
2130 @testset " Small full system" begin
22- n = 10
2331 @testset " Simple eigenvalue problem" begin
2432 @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
2533 @testset " largest = $largest " for largest in (true , false )
5967 end
6068 end
6169 end
62-
6370 @testset " Sparse Laplacian" begin
6471 A = laplace_matrix (Float64, 20 , 2 )
6572 rhs = randn (size (A, 2 ), 1 )
7481 end
7582 end
7683 end
84+ @testset " Zero initial solution" begin
85+ @testset " Simple eigenvalue problem" begin
86+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
87+ @testset " largest = $largest " for largest in (true , false )
88+ A = rand (T, n, n)
89+ A = A' * A + I
90+ b = zeros (T, n, 1 )
91+ tol = √ eps (real (T))
92+
93+ r = lobpcg (A, largest, b; tol= tol, maxiter= Inf , log= false )
94+ λ, X = r. λ, r. X
95+ @test norm (A* X - X* λ) ≤ tol
96+ end
97+ end
98+ end
99+ @testset " Generalized eigenvalue problem" begin
100+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
101+ @testset " largest = $largest " for largest in (true , false )
102+ A = rand (T, n, n)
103+ A = A' * A + I
104+ B = rand (T, n, n)
105+ B = B' * B + I
106+ b = zeros (T, n, 1 )
107+ tol = √ eps (real (T))
108+
109+ r = lobpcg (A, B, largest, b; tol= tol, maxiter= Inf , log= true )
110+ λ, X = r. λ, r. X
111+ @test max_err (A* X - B* X* λ) ≤ tol
112+ end
113+ end
114+ end
115+ end
116+ @testset " No initial solution" begin
117+ @testset " Simple eigenvalue problem" begin
118+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
119+ @testset " largest = $largest " for largest in (true , false )
120+ A = rand (T, n, n)
121+ A = A' * A + I
122+ tol = √ eps (real (T))
123+
124+ r = lobpcg (A, largest, 1 ; tol= tol, maxiter= Inf , log= false )
125+ λ, X = r. λ, r. X
126+ @test norm (A* X - X* λ) ≤ tol
127+ end
128+ end
129+ end
130+ @testset " Generalized eigenvalue problem" begin
131+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
132+ @testset " largest = $largest " for largest in (true , false )
133+ A = rand (T, n, n)
134+ A = A' * A + I
135+ B = rand (T, n, n)
136+ B = B' * B + I
137+ tol = √ eps (real (T))
138+
139+ r = lobpcg (A, B, largest, 1 ; tol= tol, maxiter= Inf , log= true )
140+ λ, X = r. λ, r. X
141+ @test max_err (A* X - B* X* λ) ≤ tol
142+ end
143+ end
144+ end
145+ end
146+ @testset " Inplace" begin
147+ @testset " Simple eigenvalue problem" begin
148+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
149+ @testset " largest = $largest " for largest in (true , false )
150+ A = rand (T, n, n)
151+ A = A' * A + I
152+ tol = √ eps (real (T))
153+ b = rand (T, n, 1 )
154+ itr = LOBPCGIterator (A, b, largest)
155+
156+ r = lobpcg! (itr; tol= tol, maxiter= Inf , log= false )
157+ λ, X = r. λ, r. X
158+ @test norm (A* X - X* λ) ≤ tol
159+ end
160+ end
161+ end
162+ @testset " Generalized eigenvalue problem" begin
163+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
164+ @testset " largest = $largest " for largest in (true , false )
165+ A = rand (T, n, n)
166+ A = A' * A + I
167+ B = rand (T, n, n)
168+ B = B' * B + I
169+ b = rand (T, n, 1 )
170+ tol = √ eps (real (T))
171+ itr = LOBPCGIterator (A, B, b, largest)
172+
173+ r = lobpcg! (itr; tol= tol, maxiter= Inf , log= true )
174+ λ, X = r. λ, r. X
175+ @test max_err (A* X - B* X* λ) ≤ tol
176+ end
177+ end
178+ end
179+ end
180+ @testset " Jacobi preconditioner" begin
181+ @testset " Simple eigenvalue problem" begin
182+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
183+ @testset " largest = $largest " for largest in (true , false )
184+ A = rand (T, n, n)
185+ A = A' * A + I
186+ tol = √ eps (real (T))
187+ P = JacobiPrec (diag (A))
188+ r = lobpcg (A, largest, 1 ; P= P, tol= tol, maxiter= Inf , log= false )
189+ λ, X = r. λ, r. X
190+ @test norm (A* X - X* λ) ≤ tol
191+ end
192+ end
193+ end
194+ @testset " Generalized eigenvalue problem" begin
195+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
196+ @testset " largest = $largest " for largest in (true , false )
197+ A = rand (T, n, n)
198+ A = A' * A + I
199+ P = JacobiPrec (diag (A))
200+ B = rand (T, n, n)
201+ B = B' * B + I
202+ tol = √ eps (real (T))
203+
204+ r = lobpcg (A, B, largest, 1 ; P= P, tol= tol, maxiter= Inf , log= true )
205+ λ, X = r. λ, r. X
206+ @test max_err (A* X - B* X* λ) ≤ tol
207+ end
208+ end
209+ end
210+ end
211+
212+ @testset " Constraint" begin
213+ @testset " Simple eigenvalue problem" begin
214+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
215+ @testset " largest = $largest " for largest in (true , false )
216+ A = rand (T, n, n)
217+ A = A' * A + I
218+ tol = √ eps (real (T))
219+ r = lobpcg (A, largest, 1 ; tol= tol, maxiter= Inf , log= false )
220+ λ1, X1 = r. λ, r. X
221+ r = lobpcg (A, largest, 1 ; C= copy (r. X), tol= tol, maxiter= Inf , log= false )
222+ λ2, X2 = r. λ, r. X
223+ @test norm (A* X2 - X2* λ2) ≤ tol
224+ @test isapprox (real (Ac_mul_B (X1, X2)[1 ,1 ]), 0 , atol= n* tol)
225+ end
226+ end
227+ end
228+ @testset " Generalized eigenvalue problem" begin
229+ @testset " Matrix{$T }" for T in (Float32, Float64, Complex64, Complex128)
230+ @testset " largest = $largest " for largest in (true , false )
231+ A = rand (T, n, n)
232+ A = A' * A + I
233+ B = rand (T, n, n)
234+ B = B' * B + I
235+ tol = √ eps (real (T))
236+ r = lobpcg (A, B, largest, 1 ; tol= tol, maxiter= Inf , log= false )
237+ λ1, X1 = r. λ, r. X
238+ r = lobpcg (A, B, largest, 1 ; C= copy (r. X), tol= tol, maxiter= Inf , log= false )
239+ λ2, X2 = r. λ, r. X
240+ @test norm (A* X2 - B* X2* λ2) ≤ tol
241+ @test isapprox (real (Ac_mul_B (X1, B* X2)[1 ,1 ]), 0 , atol= n* tol)
242+ end
243+ end
244+ end
245+ end
77246 end
78247 @testset " Two eigenvalues" begin
79248 @testset " Small full system" begin
104273 B = rand (T, n, n)
105274 B = B' * B + I
106275 b = rand (T, n, 2 )
107- tol = √ eps (real (T))
108-
276+ tol = eps (real (T))^ (real (T)(4 / 10 ))
109277 r = lobpcg (A, B, largest, b; tol= tol, maxiter= Inf , log= true )
110278 λ, X = r. λ, r. X
111279 @test max_err (A* X - B* X* diagm (λ)) ≤ tol
0 commit comments