Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "chebyshev.hpp"
5 :
6 : #include <mfem/general/forall.hpp>
7 :
8 : namespace palace
9 : {
10 :
11 : namespace
12 : {
13 :
14 0 : double GetLambdaMax(MPI_Comm comm, const Operator &A, const Vector &dinv)
15 : {
16 : // Assumes A SPD (diag(A) > 0) to use Hermitian eigenvalue solver.
17 : DiagonalOperator Dinv(dinv);
18 0 : ProductOperator DinvA(Dinv, A);
19 0 : return linalg::SpectralNorm(comm, DinvA, true);
20 : }
21 :
22 0 : double GetLambdaMax(MPI_Comm comm, const ComplexOperator &A, const ComplexVector &dinv)
23 : {
24 : // Assumes A SPD (diag(A) > 0) to use Hermitian eigenvalue solver.
25 : ComplexDiagonalOperator Dinv(dinv);
26 0 : ComplexProductOperator DinvA(Dinv, A);
27 0 : return linalg::SpectralNorm(comm, DinvA, A.IsReal());
28 : }
29 :
30 : template <bool Transpose = false>
31 : inline void ApplyOp(const Operator &A, const Vector &x, Vector &y)
32 : {
33 0 : A.Mult(x, y);
34 : }
35 :
36 : template <bool Transpose = false>
37 : inline void ApplyOp(const ComplexOperator &A, const ComplexVector &x, ComplexVector &y)
38 : {
39 : if constexpr (!Transpose)
40 : {
41 0 : A.Mult(x, y);
42 : }
43 : else
44 : {
45 : A.MultHermitianTranspose(x, y);
46 : }
47 : }
48 :
49 : template <bool Transpose = false>
50 : inline void ApplyOp(const Operator &A, const Vector &x, Vector &y, const double a)
51 : {
52 0 : A.AddMult(x, y, a);
53 : }
54 :
55 : template <bool Transpose = false>
56 : inline void ApplyOp(const ComplexOperator &A, const ComplexVector &x, ComplexVector &y,
57 : const double a)
58 : {
59 : if constexpr (!Transpose)
60 : {
61 0 : A.AddMult(x, y, a);
62 : }
63 : else
64 : {
65 : A.AddMultHermitianTranspose(x, y, a);
66 : }
67 : }
68 :
69 : template <bool Transpose = false>
70 0 : inline void ApplyOrder0(double sr, const Vector &dinv, const Vector &r, Vector &d)
71 : {
72 0 : const bool use_dev = dinv.UseDevice() || r.UseDevice() || d.UseDevice();
73 : const int N = d.Size();
74 0 : const auto *DI = dinv.Read(use_dev);
75 0 : const auto *R = r.Read(use_dev);
76 0 : auto *D = d.Write(use_dev);
77 : mfem::forall_switch(use_dev, N,
78 0 : [=] MFEM_HOST_DEVICE(int i) { D[i] = sr * DI[i] * R[i]; });
79 0 : }
80 :
81 : template <bool Transpose = false>
82 0 : inline void ApplyOrder0(const double sr, const ComplexVector &dinv, const ComplexVector &r,
83 : ComplexVector &d)
84 : {
85 0 : const bool use_dev = dinv.UseDevice() || r.UseDevice() || d.UseDevice();
86 : const int N = dinv.Size();
87 0 : const auto *DIR = dinv.Real().Read(use_dev);
88 0 : const auto *DII = dinv.Imag().Read(use_dev);
89 0 : const auto *RR = r.Real().Read(use_dev);
90 0 : const auto *RI = r.Imag().Read(use_dev);
91 0 : auto *DR = d.Real().Write(use_dev);
92 0 : auto *DI = d.Imag().Write(use_dev);
93 : if constexpr (!Transpose)
94 : {
95 : mfem::forall_switch(use_dev, N,
96 0 : [=] MFEM_HOST_DEVICE(int i)
97 : {
98 0 : DR[i] = sr * (DIR[i] * RR[i] - DII[i] * RI[i]);
99 0 : DI[i] = sr * (DII[i] * RR[i] + DIR[i] * RI[i]);
100 : });
101 : }
102 : else
103 : {
104 : mfem::forall_switch(use_dev, N,
105 : [=] MFEM_HOST_DEVICE(int i)
106 : {
107 : DR[i] = sr * (DIR[i] * RR[i] + DII[i] * RI[i]);
108 : DI[i] = sr * (-DII[i] * RR[i] + DIR[i] * RI[i]);
109 : });
110 : }
111 0 : }
112 :
113 : template <bool Transpose = false>
114 0 : inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv,
115 : const Vector &r, Vector &d)
116 : {
117 0 : const bool use_dev = dinv.UseDevice() || r.UseDevice() || d.UseDevice();
118 : const int N = dinv.Size();
119 0 : const auto *DI = dinv.Read(use_dev);
120 0 : const auto *R = r.Read(use_dev);
121 0 : auto *D = d.ReadWrite(use_dev);
122 0 : mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE(int i)
123 0 : { D[i] = sd * D[i] + sr * DI[i] * R[i]; });
124 0 : }
125 :
126 : template <bool Transpose = false>
127 0 : inline void ApplyOrderK(const double sd, const double sr, const ComplexVector &dinv,
128 : const ComplexVector &r, ComplexVector &d)
129 : {
130 0 : const bool use_dev = dinv.UseDevice() || r.UseDevice() || d.UseDevice();
131 : const int N = dinv.Size();
132 0 : const auto *DIR = dinv.Real().Read(use_dev);
133 0 : const auto *DII = dinv.Imag().Read(use_dev);
134 0 : const auto *RR = r.Real().Read(use_dev);
135 0 : const auto *RI = r.Imag().Read(use_dev);
136 0 : auto *DR = d.Real().ReadWrite(use_dev);
137 0 : auto *DI = d.Imag().ReadWrite(use_dev);
138 : if constexpr (!Transpose)
139 : {
140 : mfem::forall_switch(use_dev, N,
141 0 : [=] MFEM_HOST_DEVICE(int i)
142 : {
143 0 : DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] - DII[i] * RI[i]);
144 0 : DI[i] = sd * DI[i] + sr * (DII[i] * RR[i] + DIR[i] * RI[i]);
145 : });
146 : }
147 : else
148 : {
149 : mfem::forall_switch(use_dev, N,
150 : [=] MFEM_HOST_DEVICE(int i)
151 : {
152 : DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] + DII[i] * RI[i]);
153 : DI[i] = sd * DI[i] + sr * (-DII[i] * RR[i] + DIR[i] * RI[i]);
154 : });
155 : }
156 0 : }
157 :
158 : } // namespace
159 :
160 : template <typename OperType>
161 0 : ChebyshevSmoother<OperType>::ChebyshevSmoother(MPI_Comm comm, int smooth_it, int poly_order,
162 : double sf_max)
163 0 : : Solver<OperType>(), comm(comm), pc_it(smooth_it), order(poly_order), A(nullptr),
164 0 : lambda_max(0.0), sf_max(sf_max)
165 : {
166 0 : MFEM_VERIFY(order > 0, "Polynomial order for Chebyshev smoothing must be positive!");
167 0 : }
168 :
169 : template <typename OperType>
170 0 : void ChebyshevSmoother<OperType>::SetOperator(const OperType &op)
171 : {
172 0 : A = &op;
173 0 : d.SetSize(op.Height());
174 0 : dinv.SetSize(op.Height());
175 0 : d.UseDevice(true);
176 0 : dinv.UseDevice(true);
177 0 : op.AssembleDiagonal(dinv);
178 0 : dinv.Reciprocal();
179 :
180 : // Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. See
181 : // mfem::OperatorChebyshevSmoother or Adams et al. (2003).
182 0 : lambda_max = sf_max * GetLambdaMax(comm, *A, dinv);
183 0 : MFEM_VERIFY(lambda_max > 0.0,
184 : "Encountered zero maximum eigenvalue in Chebyshev smoother!");
185 :
186 0 : this->height = op.Height();
187 0 : this->width = op.Width();
188 0 : }
189 :
190 : template <typename OperType>
191 0 : void ChebyshevSmoother<OperType>::Mult2(const VecType &x, VecType &y, VecType &r) const
192 : {
193 : // Apply smoother: y = y + p(A) (x - A y) .
194 0 : for (int it = 0; it < pc_it; it++)
195 : {
196 0 : if (this->initial_guess || it > 0)
197 : {
198 0 : ApplyOp(*A, y, r);
199 0 : linalg::AXPBY(1.0, x, -1.0, r);
200 : }
201 : else
202 : {
203 0 : r = x;
204 0 : y = 0.0;
205 : }
206 :
207 : // 4th-kind Chebyshev smoother, from Phillips and Fischer or Lottes (with k -> k + 1
208 : // shift due to 1-based indexing).
209 0 : ApplyOrder0(4.0 / (3.0 * lambda_max), dinv, r, d);
210 0 : for (int k = 1; k < order; k++)
211 : {
212 0 : y += d;
213 0 : ApplyOp(*A, d, r, -1.0);
214 0 : const double sd = (2.0 * k - 1.0) / (2.0 * k + 3.0);
215 0 : const double sr = (8.0 * k + 4.0) / ((2.0 * k + 3.0) * lambda_max);
216 0 : ApplyOrderK(sd, sr, dinv, r, d);
217 : }
218 0 : y += d;
219 : }
220 0 : }
221 :
222 : template <typename OperType>
223 0 : ChebyshevSmoother1stKind<OperType>::ChebyshevSmoother1stKind(MPI_Comm comm, int smooth_it,
224 : int poly_order, double sf_max,
225 : double sf_min)
226 0 : : Solver<OperType>(), comm(comm), pc_it(smooth_it), order(poly_order), A(nullptr),
227 0 : theta(0.0), sf_max(sf_max), sf_min(sf_min)
228 : {
229 0 : MFEM_VERIFY(order > 0, "Polynomial order for Chebyshev smoothing must be positive!");
230 0 : }
231 :
232 : template <typename OperType>
233 0 : void ChebyshevSmoother1stKind<OperType>::SetOperator(const OperType &op)
234 : {
235 0 : A = &op;
236 0 : d.SetSize(op.Height());
237 0 : dinv.SetSize(op.Height());
238 0 : d.UseDevice(true);
239 0 : dinv.UseDevice(true);
240 0 : op.AssembleDiagonal(dinv);
241 0 : dinv.Reciprocal();
242 :
243 : // Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. The
244 : // optimized estimate of lambda_min comes from (2.24) of Phillips and Fischer (2022).
245 0 : if (sf_min <= 0.0)
246 : {
247 0 : sf_min = 1.69 / (std::pow(order, 1.68) + 2.11 * order + 1.98);
248 : }
249 0 : const double lambda_max = sf_max * GetLambdaMax(comm, *A, dinv);
250 0 : MFEM_VERIFY(lambda_max > 0.0,
251 : "Encountered zero maximum eigenvalue in Chebyshev smoother!");
252 0 : const double lambda_min = sf_min * lambda_max;
253 0 : theta = 0.5 * (lambda_max + lambda_min);
254 0 : delta = 0.5 * (lambda_max - lambda_min);
255 :
256 0 : this->height = op.Height();
257 0 : this->width = op.Width();
258 0 : }
259 :
260 : template <typename OperType>
261 0 : void ChebyshevSmoother1stKind<OperType>::Mult2(const VecType &x, VecType &y,
262 : VecType &r) const
263 : {
264 : // Apply smoother: y = y + p(A) (x - A y) .
265 0 : for (int it = 0; it < pc_it; it++)
266 : {
267 0 : if (this->initial_guess || it > 0)
268 : {
269 0 : ApplyOp(*A, y, r);
270 0 : linalg::AXPBY(1.0, x, -1.0, r);
271 : }
272 : else
273 : {
274 0 : r = x;
275 0 : y = 0.0;
276 : }
277 :
278 : // 1th-kind Chebyshev smoother, from Phillips and Fischer or Adams.
279 0 : ApplyOrder0(1.0 / theta, dinv, r, d);
280 0 : double rhop = delta / theta;
281 0 : for (int k = 1; k < order; k++)
282 : {
283 0 : y += d;
284 0 : ApplyOp(*A, d, r, -1.0);
285 0 : const double rho = 1.0 / (2.0 * theta / delta - rhop);
286 0 : const double sd = rho * rhop;
287 0 : const double sr = 2.0 * rho / delta;
288 0 : ApplyOrderK(sd, sr, dinv, r, d);
289 : rhop = rho;
290 : }
291 0 : y += d;
292 : }
293 0 : }
294 :
295 : template class ChebyshevSmoother<Operator>;
296 : template class ChebyshevSmoother<ComplexOperator>;
297 :
298 : template class ChebyshevSmoother1stKind<Operator>;
299 : template class ChebyshevSmoother1stKind<ComplexOperator>;
300 :
301 : } // namespace palace
|
