@@ -160,11 +160,7 @@ function eigQL2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matr
160
160
return c, s
161
161
end
162
162
163
- function eigvalsPWK! (
164
- S:: SymTridiagonal{T} ;
165
- tol = eps (T),
166
- debug:: Bool = false ,
167
- ) where {T<: Real }
163
+ function eigvalsPWK! (S:: SymTridiagonal{T} ; tol = eps (T)) where {T<: Real }
168
164
d = S. dv
169
165
e = S. ev
170
166
n = length (d)
@@ -187,15 +183,16 @@ function eigvalsPWK!(
187
183
188
184
# Deflate?
189
185
if blockstart == blockend
190
- # Yes
186
+ @debug " Deflate? Yes! "
191
187
blockstart += 1
192
188
elseif blockstart + 1 == blockend
193
- debug && println ( " Yes, but after sreolving 2x2 block" )
189
+ @ debug " Defalte? Yes, but after solving 2x2 block! "
194
190
d[blockstart], d[blockend] =
195
191
eigvals2x2 (d[blockstart], d[blockend], sqrt (e[blockstart]))
196
192
e[blockstart] = zero (T)
197
193
blockstart += 1
198
194
else
195
+ @debug " Deflate? No!"
199
196
# Calculate shift
200
197
sqrte = sqrt (e[blockstart])
201
198
μ = (d[blockstart+ 1 ] - d[blockstart]) / (2 * sqrte)
@@ -205,15 +202,9 @@ function eigvalsPWK!(
205
202
# QL bulk chase
206
203
singleShiftQLPWK! (S, μ, blockstart, blockend)
207
204
208
- debug && @printf (
209
- " QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %e, rotations: %d\n " ,
210
- blockstart,
211
- blockend,
212
- e[blockstart],
213
- e[blockend- 1 ],
214
- μ,
215
- iter += blockend - blockstart
216
- )
205
+ rotations = blockend - blockstart
206
+ iter += rotations
207
+ @debug " QL, blockstart, blockend, e[blockstart], e[blockend-1], μ, rotations" blockstart blockend e[blockstart] e[blockend- 1 ] μ rotations
217
208
end
218
209
if blockstart >= n
219
210
break
@@ -227,7 +218,6 @@ function eigQL!(
227
218
S:: SymTridiagonal{T} ;
228
219
vectors:: Matrix = zeros (T, 0 , size (S, 1 )),
229
220
tol = eps (T),
230
- debug:: Bool = false ,
231
221
) where {T<: Real }
232
222
d = S. dv
233
223
e = S. ev
@@ -262,28 +252,22 @@ function eigQL!(
262
252
end
263
253
# Deflate?
264
254
if blockstart == blockend
265
- # Yes!
255
+ @debug " Deflate? Yes!"
266
256
blockstart += 1
267
257
elseif blockstart + 1 == blockend
268
- debug && println ( " Yes, but after solving 2x2 block" )
258
+ @ debug " Deflate? Yes, but after solving 2x2 block"
269
259
eigQL2x2! (d, e, blockstart, vectors)
270
260
blockstart += 1
271
261
else
262
+ @debug " Deflate? No!"
272
263
# Calculate shift
273
264
μ = (d[blockstart+ 1 ] - d[blockstart]) / (2 * e[blockstart])
274
265
r = hypot (μ, one (T))
275
266
μ = d[blockstart] - (e[blockstart] / (μ + copysign (r, μ)))
276
267
277
268
# QL bulk chase
278
269
singleShiftQL! (S, μ, blockstart, blockend, vectors)
279
- debug && @printf (
280
- " QL, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, μ: %f\n " ,
281
- blockstart,
282
- blockend,
283
- e[blockstart],
284
- e[blockend- 1 ],
285
- μ
286
- )
270
+ @debug " QL, blockstart, blockend, e[blockstart], e[blockend-1] μ" blockstart blockend e[blockstart] e[blockend- 1 ] μ
287
271
end
288
272
if blockstart >= n
289
273
break
@@ -298,7 +282,6 @@ function eigQR!(
298
282
S:: SymTridiagonal{T} ;
299
283
vectors:: Matrix = zeros (T, 0 , size (S, 1 )),
300
284
tol = eps (T),
301
- debug:: Bool = false ,
302
285
) where {T<: Real }
303
286
d = S. dv
304
287
e = S. ev
@@ -317,13 +300,14 @@ function eigQR!(
317
300
end
318
301
# Deflate?
319
302
if blockstart == blockend
320
- # Yes!
303
+ @debug " Deflate? Yes!"
321
304
blockstart += 1
322
305
elseif blockstart + 1 == blockend
323
- debug && println ( " Yes, but after solving 2x2 block" )
306
+ @ debug " Deflate? Yes, but after solving 2x2 block! "
324
307
eigQL2x2! (d, e, blockstart, vectors)
325
308
blockstart += 1
326
309
else
310
+ @debug " Deflate? No!"
327
311
# Calculate shift
328
312
μ = (d[blockend- 1 ] - d[blockend]) / (2 * e[blockend- 1 ])
329
313
r = hypot (μ, one (T))
@@ -332,15 +316,7 @@ function eigQR!(
332
316
# QR bulk chase
333
317
singleShiftQR! (S, μ, blockstart, blockend, vectors)
334
318
335
- debug && @printf (
336
- " QR, blockstart: %d, blockend: %d, e[blockstart]: %e, e[blockend-1]:%e, d[blockend]: %f, μ: %f\n " ,
337
- blockstart,
338
- blockend,
339
- e[blockstart],
340
- e[blockend- 1 ],
341
- d[blockend],
342
- μ
343
- )
319
+ @debug " QR, blockstart, blockend, e[blockstart], e[blockend-1], d[blockend], μ" blockstart blockend e[blockstart] e[blockend- 1 ] d[blockend] μ
344
320
end
345
321
if blockstart >= n
346
322
break
@@ -601,87 +577,66 @@ function symtriUpper!(
601
577
end
602
578
603
579
604
- _eigvals! (A:: SymmetricTridiagonalFactorization ; tol = eps (real (eltype (A))), debug = false ) =
605
- eigvalsPWK! (A. diagonals, tol = tol, debug = debug )
580
+ _eigvals! (A:: SymmetricTridiagonalFactorization ; tol = eps (real (eltype (A)))) =
581
+ eigvalsPWK! (A. diagonals, tol = tol)
606
582
607
- _eigvals! (A:: SymTridiagonal ; tol = eps (real (eltype (A))), debug = false ) =
608
- eigvalsPWK! (A, tol = tol, debug = debug)
583
+ _eigvals! (A:: SymTridiagonal ; tol = eps (real (eltype (A)))) = eigvalsPWK! (A, tol = tol)
609
584
610
- _eigvals! (A:: Hermitian ; tol = eps (real (eltype (A))), debug = false ) =
611
- eigvals! (symtri! (A), tol = tol, debug = debug)
585
+ _eigvals! (A:: Hermitian ; tol = eps (real (eltype (A)))) = eigvals! (symtri! (A), tol = tol)
612
586
613
587
614
- LinearAlgebra. eigvals! (
615
- A:: SymmetricTridiagonalFactorization ;
616
- tol = eps (real (eltype (A))),
617
- debug = false ,
618
- ) = _eigvals! (A, tol = tol, debug = debug)
588
+ LinearAlgebra. eigvals! (A:: SymmetricTridiagonalFactorization ; tol = eps (real (eltype (A)))) =
589
+ _eigvals! (A, tol = tol)
619
590
620
- LinearAlgebra. eigvals! (A:: SymTridiagonal ; tol = eps (real (eltype (A))), debug = false ) =
621
- _eigvals! (A, tol = tol, debug = debug )
591
+ LinearAlgebra. eigvals! (A:: SymTridiagonal ; tol = eps (real (eltype (A)))) =
592
+ _eigvals! (A, tol = tol)
622
593
623
- LinearAlgebra. eigvals! (A:: Hermitian ; tol = eps (real (eltype (A))), debug = false ) =
624
- _eigvals! (A, tol = tol, debug = debug)
594
+ LinearAlgebra. eigvals! (A:: Hermitian ; tol = eps (real (eltype (A)))) = _eigvals! (A, tol = tol)
625
595
626
596
627
- _eigen! (A:: SymmetricTridiagonalFactorization ; tol = eps (real (eltype (A))), debug = false ) =
628
- LinearAlgebra. Eigen (
629
- eigQL! (A. diagonals, vectors = Array (A. Q), tol = tol, debug = debug)... ,
630
- )
597
+ _eigen! (A:: SymmetricTridiagonalFactorization ; tol = eps (real (eltype (A)))) =
598
+ LinearAlgebra. Eigen (eigQL! (A. diagonals, vectors = Array (A. Q), tol = tol)... )
631
599
632
- _eigen! (A:: SymTridiagonal ; tol = eps (real (eltype (A))), debug = false ) = LinearAlgebra. Eigen (
633
- eigQL! (
634
- A,
635
- vectors = Matrix {eltype(A)} (I, size (A, 1 ), size (A, 1 )),
636
- tol = tol,
637
- debug = debug,
638
- )... ,
600
+ _eigen! (A:: SymTridiagonal ; tol = eps (real (eltype (A)))) = LinearAlgebra. Eigen (
601
+ eigQL! (A, vectors = Matrix {eltype(A)} (I, size (A, 1 ), size (A, 1 )), tol = tol)... ,
639
602
)
640
603
641
- _eigen! (A:: Hermitian ; tol = eps (real (eltype (A))), debug = false ) =
642
- _eigen! (symtri! (A), tol = tol, debug = debug)
604
+ _eigen! (A:: Hermitian ; tol = eps (real (eltype (A)))) = _eigen! (symtri! (A), tol = tol)
643
605
644
606
645
- LinearAlgebra. eigen! (
646
- A:: SymmetricTridiagonalFactorization ;
647
- tol = eps (real (eltype (A))),
648
- debug = false ,
649
- ) = _eigen! (A, tol = tol, debug = debug)
607
+ LinearAlgebra. eigen! (A:: SymmetricTridiagonalFactorization ; tol = eps (real (eltype (A)))) =
608
+ _eigen! (A, tol = tol)
650
609
651
- LinearAlgebra. eigen! (A:: SymTridiagonal ; tol = eps (real (eltype (A))), debug = false ) =
652
- _eigen! (A, tol = tol, debug = debug)
610
+ LinearAlgebra. eigen! (A:: SymTridiagonal ; tol = eps (real (eltype (A)))) = _eigen! (A, tol = tol)
653
611
654
- LinearAlgebra. eigen! (A:: Hermitian ; tol = eps (real (eltype (A))), debug = false ) =
655
- _eigen! (A, tol = tol, debug = debug)
612
+ LinearAlgebra. eigen! (A:: Hermitian ; tol = eps (real (eltype (A)))) = _eigen! (A, tol = tol)
656
613
657
614
658
615
function eigen2! (
659
616
A:: SymmetricTridiagonalFactorization ;
660
617
tol = eps (real (float (one (eltype (A))))),
661
- debug = false ,
662
618
)
663
619
V = zeros (eltype (A), 2 , size (A, 1 ))
664
620
V[1 ] = 1
665
621
V[end ] = 1
666
- eigQL! (A. diagonals, vectors = rmul! (V, A. Q), tol = tol, debug = debug )
622
+ eigQL! (A. diagonals, vectors = rmul! (V, A. Q), tol = tol)
667
623
end
668
624
669
- function eigen2! (A:: SymTridiagonal ; tol = eps (real (float (one (eltype (A))))), debug = false )
625
+ function eigen2! (A:: SymTridiagonal ; tol = eps (real (float (one (eltype (A))))))
670
626
V = zeros (eltype (A), 2 , size (A, 1 ))
671
627
V[1 ] = 1
672
628
V[end ] = 1
673
- eigQL! (A, vectors = V, tol = tol, debug = debug )
629
+ eigQL! (A, vectors = V, tol = tol)
674
630
end
675
631
676
- eigen2! (A:: Hermitian ; tol = eps (float (real (one (eltype (A))))), debug = false ) =
677
- eigen2! (symtri! (A), tol = tol, debug = debug )
632
+ eigen2! (A:: Hermitian ; tol = eps (float (real (one (eltype (A)))))) =
633
+ eigen2! (symtri! (A), tol = tol)
678
634
679
635
680
- eigen2 (A:: SymTridiagonal ; tol = eps (float (real (one (eltype (A))))), debug = false ) =
681
- eigen2! (copy (A), tol = tol, debug = debug )
636
+ eigen2 (A:: SymTridiagonal ; tol = eps (float (real (one (eltype (A)))))) =
637
+ eigen2! (copy (A), tol = tol)
682
638
683
- eigen2 (A:: Hermitian , tol = eps (float (real (one (eltype (A))))), debug = false ) =
684
- eigen2! (copy (A), tol = tol, debug = debug)
639
+ eigen2 (A:: Hermitian , tol = eps (float (real (one (eltype (A)))))) = eigen2! (copy (A), tol = tol)
685
640
686
641
# First method of each type here is identical to the method defined in
687
642
# LinearAlgebra but is needed for disambiguation
0 commit comments