28 #include <shylu_internal_gmres_tools.h> 35 void ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
42 template <
typename Scalar>
43 void GeneratePlaneRotation(
const Scalar &dx,
const Scalar &dy, Scalar &cs,
49 }
else if (std::abs(dy) > std::abs(dx)) {
50 Scalar temp = dx / dy;
51 sn = 1.0 / std::sqrt( 1.0 + temp * temp );
54 Scalar temp = dy / dx;
55 cs = 1.0 / std::sqrt( 1.0 + temp * temp );
62 template <
typename Scalar>
63 void ApplyPlaneRotation(Scalar &dx, Scalar &dy,
const Scalar &cs,
66 Scalar temp = cs * dx + sn * dy;
67 dy = -sn * dx + cs * dy;
73 template <
typename LocalMatrix,
typename LocalVector,
typename MultiVector >
74 void Update(MultiVector &x,
const int k,
const LocalMatrix &h,
75 const LocalVector &s,
const MultiVector &v)
80 for (
int i = k; i >= 0; i--) {
82 for (
int j = i - 1; j >= 0; j--) {
83 y[j] -= h[j][i] * y[i];
87 for (
int j = 0; j <= k; j++) {
88 x.Update(y[j], *v(j), 1.0);
92 template <
typename Operator,
typename MultiVector,
typename LeftPrec,
93 typename RightPrec,
typename GMRESManager,
typename LocalVector,
95 int GMRES(
const Operator &A, MultiVector &x,
const MultiVector &b,
96 LeftPrec *L, RightPrec *M, GMRESManager &G,
int &max_iter,
104 int i(0), j(1), k(0);
109 LocalVector s(G.restart + 1);
110 MultiVector w(b.Map(), 1,
true);
113 L->ApplyInverse(b, w);
116 MultiVector t(b.Map(), 1,
true);
118 w.Update(1.0, b, -1.0, t, 0.0);
120 MultiVector r(b.Map(), 1,
true);
121 L->ApplyInverse(w, r);
129 if ((resid = beta / normb) <= tol) {
135 while (j <= max_iter) {
136 MultiVector* v0 = (*G.v)(0);
137 v0->Update(1.0 / beta, r, 0.0);
138 s.assign(G.restart + 1, 0.0);
141 for (i = 0; i < G.restart && j <= max_iter; i++, j++) {
142 M->ApplyInverse(*((*G.v)(i)), t);
144 L->ApplyInverse(r, w);
146 for (k = 0; k <= i; k++) {
147 MultiVector* vk = (*G.v)(k);
148 w.Dot(*vk, &(G.H[k][i]));
149 w.Update(-G.H[k][i], *vk, 1.0);
151 w.Norm2(&(G.H[i + 1][i]));
152 MultiVector* vi1 = (*G.v)(i + 1);
154 vi1->Scale(1.0 / G.H[i + 1][i], w);
156 for (k = 0; k < i; k++) {
157 ApplyPlaneRotation(G.H[k][i], G.H[k + 1][i], G.cs[k], G.sn[k]);
161 GeneratePlaneRotation(G.H[i][i], G.H[i + 1][i], G.cs[i], G.sn[i]);
163 ApplyPlaneRotation(G.H[i][i], G.H[i + 1][i], G.cs[i], G.sn[i]);
164 ApplyPlaneRotation(s[i], s[i + 1], G.cs[i], G.sn[i]);
168 if ((resid = abs(s[i + 1]) / normb) < tol) {
169 MultiVector y(b.Map(), 1,
true);
170 Update(y, i, G.H, s, *(G.v));
171 M->ApplyInverse(y, t);
172 x.Update(1.0, t, 1.0);
180 MultiVector y(b.Map(), 1,
true);
181 Update(y, i - 1, G.H, s, *(G.v));
182 M->ApplyInverse(y, t);
183 x.Update(1.0, t, 1.0);
185 w.Update(1.0, b, -1.0, t, 0.0);
186 L->ApplyInverse(w, r);
189 if ((resid = beta / normb) < tol) {
204 #endif // IQR_GMRES_H