@@ -2,29 +2,13 @@ using Printf
22using LinearAlgebra
33using LinearAlgebra: givensAlgorithm
44
5- struct SymmetricTridiagonalFactorization{T} <: Factorization{T}
6- uplo:: Char
7- factors:: Matrix{T}
8- τ:: Vector{T}
9- diagonals:: SymTridiagonal
10- end
11-
12- Base. size (S:: SymmetricTridiagonalFactorization , i:: Integer ) = size (S. factors, i)
13-
5+ # # EigenQ
146struct EigenQ{T} <: AbstractMatrix{T}
157 uplo:: Char
168 factors:: Matrix{T}
179 τ:: Vector{T}
1810end
1911
20- function Base. getproperty (S:: SymmetricTridiagonalFactorization , s:: Symbol )
21- if s == :Q
22- return EigenQ (S. uplo, S. factors, S. τ)
23- else
24- return getfield (S, s)
25- end
26- end
27-
2812Base. size (Q:: EigenQ ) = (size (Q. factors, 1 ), size (Q. factors, 1 ))
2913function Base. size (Q:: EigenQ , i:: Integer )
3014 if i < 1
@@ -42,9 +26,8 @@ function LinearAlgebra.lmul!(Q::EigenQ, B::StridedVecOrMat)
4226 throw (DimensionMismatch (" " ))
4327 end
4428 if Q. uplo == ' L'
45- for k = length (Q. τ): - 1 : 1
29+ @inbounds for k = length (Q. τ): - 1 : 1
4630 for j = 1 : size (B, 2 )
47- b = view (B, :, j)
4831 vb = B[k+ 1 , j]
4932 for i = (k+ 2 ): m
5033 vb += Q. factors[i, k]' B[i, j]
@@ -57,9 +40,8 @@ function LinearAlgebra.lmul!(Q::EigenQ, B::StridedVecOrMat)
5740 end
5841 end
5942 elseif Q. uplo == ' U'
60- for k = length (Q. τ): - 1 : 1
43+ @inbounds for k = length (Q. τ): - 1 : 1
6144 for j = 1 : size (B, 2 )
62- b = view (B, :, j)
6345 vb = B[k+ 1 , j]
6446 for i = (k+ 2 ): m
6547 vb += Q. factors[k, i] * B[i, j]
120102
121103Base. Array (Q:: EigenQ ) = lmul! (Q, Matrix {eltype(Q)} (I, size (Q, 1 ), size (Q, 1 )))
122104
105+
106+ # # SymmetricTridiagonalFactorization
107+ struct SymmetricTridiagonalFactorization{T,Treal,S} <: Factorization{T}
108+ reflectors:: EigenQ{T}
109+ diagonals:: SymTridiagonal{Treal,S}
110+ end
111+
112+ Base. size (S:: SymmetricTridiagonalFactorization , i:: Integer ) = size (S. reflectors. factors, i)
113+
114+ function Base. getproperty (S:: SymmetricTridiagonalFactorization , s:: Symbol )
115+ if s == :Q
116+ return S. reflectors
117+ else
118+ return getfield (S, s)
119+ end
120+ end
121+
122+ # # Eigen solvers
123123function _updateVectors! (c, s, j, vectors)
124124 @inbounds for i = 1 : size (vectors, 1 )
125125 v1 = vectors[i, j+ 1 ]
@@ -510,9 +510,11 @@ function symtriLower!(
510510 end
511511 end
512512 SymmetricTridiagonalFactorization (
513- ' L' ,
514- AS,
515- τ,
513+ EigenQ (
514+ ' L' ,
515+ AS,
516+ τ,
517+ ),
516518 SymTridiagonal (real (diag (AS)), real (diag (AS, - 1 ))),
517519 )
518520end
@@ -564,9 +566,11 @@ function symtriUpper!(
564566 end
565567 end
566568 SymmetricTridiagonalFactorization (
567- ' U' ,
568- AS,
569- τ,
569+ EigenQ (
570+ ' U' ,
571+ AS,
572+ τ,
573+ ),
570574 SymTridiagonal (real (diag (AS)), real (diag (AS, 1 ))),
571575 )
572576end
0 commit comments