@@ -91,3 +91,131 @@ function pruneeigs(λ,V,ds,tolerance)
91
91
end
92
92
retλ,retV
93
93
end
94
+
95
+ abstract type EigenSystem end
96
+
97
+ """
98
+ SymmetricEigensystem(L::Operator, B::Operator, QuotientSpaceType::Type = QuotientSpace)
99
+
100
+ Represent the eigensystem `L v = λ v` subject to `B v = 0`, where `L` is self-adjoint with respect to the standard `L2`
101
+ inner product given the boundary conditions `B`.
102
+
103
+ !!! note
104
+ No tests are performed to assert that the operator `L` is self-adjoint, and it's the user's responsibility
105
+ to ensure that the operators are compliant.
106
+
107
+ The optional argument `QuotientSpaceType` specifies the type of space to be used to denote the quotient space in the basis
108
+ recombination process. In most cases, the default choice of `QuotientSpace` is a good one. In specific instances where `B`
109
+ is rank-deficient (e.g. it contains a column of zeros,
110
+ which typically happens if one of the basis elements already satiafies the boundary conditions),
111
+ one may need to choose this to be a `PathologicalQuotientSpace`.
112
+
113
+ !!! note
114
+ No checks on the rank of `B` are carried out, and it's up to the user to specify the correct type.
115
+ """
116
+ SymmetricEigensystem
117
+
118
+ """
119
+ SkewSymmetricEigensystem(L::Operator, B::Operator, QuotientSpaceType::Type = QuotientSpace)
120
+
121
+ Represent the eigensystem `L v = λ v` subject to `B v = 0`, where `L` is skew-symmetric with respect to the standard `L2`
122
+ inner product given the boundary conditions `B`.
123
+
124
+ !!! note
125
+ No tests are performed to assert that the operator `L` is skew-symmetric, and it's the user's responsibility
126
+ to ensure that the operators are compliant.
127
+
128
+ The optional argument `QuotientSpaceType` specifies the type of space to be used to denote the quotient space in the basis
129
+ recombination process. In most cases, the default choice of `QuotientSpace` is a good one. In specific instances where `B`
130
+ is rank-deficient (e.g. it contains a column of zeros,
131
+ which typically happens if one of the basis elements already satiafies the boundary conditions),
132
+ one may need to choose this to be a `PathologicalQuotientSpace`.
133
+
134
+ !!! note
135
+ No checks on the rank of `B` are carried out, and it's up to the user to specify the correct type.
136
+ """
137
+ SkewSymmetricEigensystem
138
+
139
+ for SET in (:SymmetricEigensystem , :SkewSymmetricEigensystem )
140
+ @eval begin
141
+ struct $ SET{LT,QST} <: EigenSystem
142
+ L :: LT
143
+ QS :: QST
144
+
145
+ function $SET (L, B, :: Type{QST} = QuotientSpace) where {QST}
146
+ L2, B2 = promotedomainspace ((L, B))
147
+ if isambiguous (domainspace (L))
148
+ throw (ArgumentError (" could not detect spaces, please specify the domain spaces for the operators" ))
149
+ end
150
+
151
+ QS = QST (B2)
152
+ new {typeof(L2),typeof(QS)} (L2, QS)
153
+ end
154
+ end
155
+ end
156
+ end
157
+
158
+ function basis_recombination (SE:: EigenSystem )
159
+ L, QS = SE. L, SE. QS
160
+ S = domainspace (L)
161
+ D1 = Conversion_normalizedspace (S)
162
+ D2 = Conversion_normalizedspace (S, Val (:backward ))
163
+ Q = Conversion (QS, S)
164
+ R = D1* Q;
165
+ C = Conversion (S, rangespace (L))
166
+ P = cache (PartialInverseOperator (C, (0 , bandwidth (L, 1 ) + bandwidth (R, 1 ) + bandwidth (C, 2 ))));
167
+ A = R' D1* P* L* D2* R
168
+ BB = R' R;
169
+
170
+ return A, BB
171
+ end
172
+
173
+ """
174
+ bandmatrices_eigen(S::Union{SymmetricEigensystem, SkewSymmetricEigensystem}, n::Integer)
175
+
176
+ Recast the symmetric/skew-symmetric eigenvalue problem `L v = λ v` subject to `B v = 0` to the generalized
177
+ eigenvalue problem `SA v = λ SB v`, where `SA` and `SB` are banded operators, and
178
+ return the `n × n` matrix representations of `SA` and `SB`.
179
+ If `S isa SymmetricEigensystem`, the returned matrices will be `Symmetric`.
180
+
181
+ !!! note
182
+ No tests are performed to assert that the system is symmetric/skew-symmetric, and it's the user's responsibility
183
+ to ensure that the operators are compliant.
184
+ """
185
+ function bandmatrices_eigen (S:: SymmetricEigensystem , n:: Integer )
186
+ A, B = _bandmatrices_eigen (S, n)
187
+ SA = Symmetric (A, :L )
188
+ SB = Symmetric (B, :L )
189
+ return SA, SB
190
+ end
191
+
192
+ function bandmatrices_eigen (S:: EigenSystem , n:: Integer )
193
+ A, B = _bandmatrices_eigen (S, n)
194
+ A2 = tril (A, bandwidth (A,1 ))
195
+ B2 = tril (B, bandwidth (B,1 ))
196
+ Matrix (A2), Matrix (B2)
197
+ end
198
+
199
+ function _bandmatrices_eigen (S:: EigenSystem , n:: Integer )
200
+ AA, BB = basis_recombination (S)
201
+ A = AA[1 : n,1 : n]
202
+ B = BB[1 : n,1 : n]
203
+ return A, B
204
+ end
205
+
206
+ function eigvals (S:: EigenSystem , n:: Integer )
207
+ SA, SB = bandmatrices_eigen (S, n)
208
+ eigvals (SA, SB)
209
+ end
210
+
211
+ function eigs (Seig:: EigenSystem , n:: Integer ; tolerance:: Float64 = 100 eps ())
212
+ SA, SB = bandmatrices_eigen (Seig, n)
213
+ λ, v = eigen (SA, SB)
214
+ vm = Matrix (v)
215
+ L, QS = Seig. L, Seig. QS
216
+ S = domainspace (L)
217
+ Q = Conversion (QS, S)
218
+ QM = Q[FiniteRange, axes (vm, 1 )]
219
+ V = QM * vm
220
+ pruneeigs (λ,V,S,tolerance)
221
+ end
0 commit comments