1 #ifndef IQR_GMRES_TOOLS_H 2 #define IQR_GMRES_TOOLS_H 5 #include <shylu_internal_gmres.h> 7 #include <Epetra_Operator.h> 25 template <
typename Map,
typename MultiVector,
typename LocalMatrix,
typename LocalVector >
33 GMRESManager(
const Map& map,
const int &rest,
const int flex = 0,
const bool scaling =
true);
37 int initialGuess(
const MultiVector &b,
38 MultiVector& X)
const;
40 int ApplyInverse(
const MultiVector &b, MultiVector& X)
const;
67 template <
typename Map,
typename MultiVector,
typename LocalMatrix,
typename LocalVector >
75 H(rest + 1, LocalVector(rest, 0.0)),
80 std::cout <<
"shylu_internal_gmres_tools: restart must be positive" << std::endl;
85 v =
new MultiVector(map_, rest + 1,
true);
90 w =
new MultiVector(map_, rest + 1,
true);
94 template <
typename Map,
typename MultiVector,
typename LocalMatrix,
typename LocalVector >
95 GMRESManager< Map, MultiVector, LocalMatrix, LocalVector >::~GMRESManager()
104 template <
typename Map,
typename MultiVector,
typename LocalMatrix,
typename LocalVector >
105 int GMRESManager< Map, MultiVector, LocalMatrix, LocalVector >::initialGuess(
const MultiVector &b,
106 MultiVector& X)
const 111 std::cout <<
"GMRES_TOOLS.initialGuess: mm>m, not expected using m" 116 if (P != 0 && isFlex == 0) {
117 this->ApplyInverse(b, X);
131 LocalVector bp(mm + 1, 0.0);
133 for (i = 0; i < mm + 1; i++) {
134 b.Dot(*(*v)(i), &bp[i]);
138 for (i = 0; i < mm; i++) {
139 ApplyPlaneRotation(bp[i], bp[i + 1], cs[i], sn[i]);
143 Update(X, mm-1, H, bp, *w);
150 template <
typename Map,
typename MultiVector,
typename LocalMatrix,
typename LocalVector >
151 int GMRESManager< Map, MultiVector, LocalMatrix, LocalVector >::start()
164 template <
typename Map,
typename MultiVector,
typename LocalMatrix,
typename LocalVector >
165 int GMRESManager< Map, MultiVector,
166 LocalMatrix, LocalVector >::ApplyInverse(
const MultiVector &b,
167 MultiVector& X)
const 172 std::cout <<
"GMRES_TOOLS.solve: mm>m, not expected using m" 178 std::cout <<
"GMRES_TOOLS.solve: no preconditioner" << std::endl;
188 LocalVector bp(mm + 1, 0.0);
189 MultiVector x(map_, 1,
false);
193 for (i = 0; i < mm + 1; i++) {
194 MultiVector* vi = (*v)(i);
198 x.Update(-bp[i], *vi, 1.0);
202 for (i = 0; i < mm; i++) {
203 ApplyPlaneRotation(bp[i], bp[i + 1], cs[i], sn[i]);
206 typedef typename LocalVector::value_type Scalar;
208 for (i = 0; i < mm; i++)
212 scaling =
static_cast<Scalar
>(mm) / scaling;
214 LocalVector bq(mm + 1, 0.0);
218 for (i = mm-1; i >= 0; i--) {
219 ApplyPlaneRotation(bq[i], bq[i + 1], cs[i], -sn[i]);
225 for (i = 0; i <= mm; i++) {
226 MultiVector* vi = (*v)(i);
228 x.Update(bq[i], *vi, 1.0);
236 Update(x, mm-1, H, bp, *w);
239 if (P2 !=0 && isFlex == 0) {
240 P2->ApplyInverse(x, X);
245 if (P != 0 && isFlex == 0) {
246 P->ApplyInverse(x, X);
256 #endif // IQR_GMRES_TOOLS_H