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 "distrelaxation.hpp"
5 :
6 : #include <mfem.hpp>
7 : #include "linalg/chebyshev.hpp"
8 : #include "linalg/rap.hpp"
9 :
10 : namespace palace
11 : {
12 :
13 : template <typename OperType>
14 0 : DistRelaxationSmoother<OperType>::DistRelaxationSmoother(
15 : MPI_Comm comm, const Operator &G, int smooth_it, int cheby_smooth_it, int cheby_order,
16 : double cheby_sf_max, double cheby_sf_min, bool cheby_4th_kind)
17 0 : : Solver<OperType>(), pc_it(smooth_it), G(&G), A(nullptr), A_G(nullptr),
18 0 : dbc_tdof_list_G(nullptr)
19 : {
20 : // Initialize smoothers.
21 0 : if (cheby_4th_kind)
22 : {
23 0 : B = std::make_unique<ChebyshevSmoother<OperType>>(comm, cheby_smooth_it, cheby_order,
24 : cheby_sf_max);
25 0 : B_G = std::make_unique<ChebyshevSmoother<OperType>>(comm, cheby_smooth_it, cheby_order,
26 : cheby_sf_max);
27 : }
28 : else
29 : {
30 0 : B = std::make_unique<ChebyshevSmoother1stKind<OperType>>(
31 : comm, cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min);
32 0 : B_G = std::make_unique<ChebyshevSmoother1stKind<OperType>>(
33 : comm, cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min);
34 : }
35 0 : B_G->SetInitialGuess(false);
36 0 : }
37 :
38 : template <typename OperType>
39 0 : void DistRelaxationSmoother<OperType>::SetOperators(const OperType &op,
40 : const OperType &op_G)
41 : {
42 : using ParOperType =
43 : typename std::conditional<std::is_same<OperType, ComplexOperator>::value,
44 : ComplexParOperator, ParOperator>::type;
45 :
46 0 : MFEM_VERIFY(op.Height() == G->Height() && op.Width() == G->Height() &&
47 : op_G.Height() == G->Width() && op_G.Width() == G->Width(),
48 : "Invalid operator sizes for DistRelaxationSmoother!");
49 0 : A = &op;
50 0 : A_G = &op_G;
51 0 : x_G.SetSize(op_G.Height());
52 0 : y_G.SetSize(op_G.Height());
53 0 : r_G.SetSize(op_G.Height());
54 0 : x_G.UseDevice(true);
55 0 : y_G.UseDevice(true);
56 0 : r_G.UseDevice(true);
57 :
58 0 : const auto *PtAP_G = dynamic_cast<const ParOperType *>(&op_G);
59 0 : MFEM_VERIFY(PtAP_G,
60 : "ChebyshevSmoother requires a ParOperator or ComplexParOperator operator!");
61 0 : dbc_tdof_list_G = PtAP_G->GetEssentialTrueDofs();
62 :
63 : // Set up smoothers for A and A_G.
64 0 : B->SetOperator(op);
65 0 : B_G->SetOperator(op_G);
66 :
67 0 : this->height = op.Height();
68 0 : this->width = op.Width();
69 0 : }
70 :
71 : namespace
72 : {
73 :
74 : inline void RealAddMult(const Operator &op, const Vector &x, Vector &y)
75 : {
76 0 : op.AddMult(x, y, 1.0);
77 : }
78 :
79 : inline void RealAddMult(const Operator &op, const ComplexVector &x, ComplexVector &y)
80 : {
81 0 : op.AddMult(x.Real(), y.Real(), 1.0);
82 0 : op.AddMult(x.Imag(), y.Imag(), 1.0);
83 : }
84 :
85 : inline void RealMultTranspose(const Operator &op, const Vector &x, Vector &y)
86 : {
87 0 : op.MultTranspose(x, y);
88 0 : }
89 :
90 : inline void RealMultTranspose(const Operator &op, const ComplexVector &x, ComplexVector &y)
91 : {
92 0 : op.MultTranspose(x.Real(), y.Real());
93 0 : op.MultTranspose(x.Imag(), y.Imag());
94 0 : }
95 :
96 : } // namespace
97 :
98 : template <typename OperType>
99 0 : void DistRelaxationSmoother<OperType>::Mult2(const VecType &x, VecType &y, VecType &r) const
100 : {
101 : // Apply smoother.
102 0 : for (int it = 0; it < pc_it; it++)
103 : {
104 : // y = y + B (x - A y)
105 0 : B->SetInitialGuess(this->initial_guess || it > 0);
106 0 : B->Mult2(x, y, r);
107 :
108 : // y = y + G B_G Gᵀ (x - A y)
109 0 : A->Mult(y, r);
110 0 : linalg::AXPBY(1.0, x, -1.0, r);
111 0 : RealMultTranspose(*G, r, x_G);
112 0 : if (dbc_tdof_list_G)
113 : {
114 0 : linalg::SetSubVector(x_G, *dbc_tdof_list_G, 0.0);
115 : }
116 0 : B_G->Mult2(x_G, y_G, r_G);
117 0 : RealAddMult(*G, y_G, y);
118 : }
119 0 : }
120 :
121 : template <typename OperType>
122 0 : void DistRelaxationSmoother<OperType>::MultTranspose2(const VecType &x, VecType &y,
123 : VecType &r) const
124 : {
125 : // Apply transpose.
126 0 : B->SetInitialGuess(true);
127 0 : for (int it = 0; it < pc_it; it++)
128 : {
129 : // y = y + G B_Gᵀ Gᵀ (x - A y)
130 0 : if (this->initial_guess || it > 0)
131 : {
132 0 : A->Mult(y, r);
133 0 : linalg::AXPBY(1.0, x, -1.0, r);
134 0 : RealMultTranspose(*G, r, x_G);
135 : }
136 : else
137 : {
138 0 : y = 0.0;
139 0 : RealMultTranspose(*G, x, x_G);
140 : }
141 0 : if (dbc_tdof_list_G)
142 : {
143 0 : linalg::SetSubVector(x_G, *dbc_tdof_list_G, 0.0);
144 : }
145 0 : B_G->MultTranspose2(x_G, y_G, r_G);
146 0 : RealAddMult(*G, y_G, y);
147 :
148 : // y = y + Bᵀ (x - A y)
149 0 : B->MultTranspose2(x, y, r);
150 : }
151 0 : }
152 :
153 : template class DistRelaxationSmoother<Operator>;
154 : template class DistRelaxationSmoother<ComplexOperator>;
155 :
156 : } // namespace palace
|
