Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
krylov_dh.c
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Euclid_dh.h"
44#include "krylov_dh.h"
45#include "Mem_dh.h"
46#include "Parser_dh.h"
47#include "Mat_dh.h"
48
49
50#undef __FUNC__
51#define __FUNC__ "bicgstab_euclid"
52void
53bicgstab_euclid (Mat_dh A, Euclid_dh ctx, double *x, double *b, int *itsOUT)
54{
55 START_FUNC_DH int its, m = ctx->m;
56 bool monitor;
57 int maxIts = ctx->maxIts;
58 double atol = ctx->atol, rtol = ctx->rtol;
59
60 /* scalars */
61 double alpha, alpha_1,
62 beta_1,
63 widget, widget_1, rho_1, rho_2, s_norm, eps, exit_a, b_iprod, r_iprod;
64
65 /* vectors */
66 double *t, *s, *s_hat, *v, *p, *p_hat, *r, *r_hat;
67
68 monitor = Parser_dhHasSwitch (parser_dh, "-monitor");
69
70 /* allocate working space */
71 t = (double *) MALLOC_DH (m * sizeof (double));
72 s = (double *) MALLOC_DH (m * sizeof (double));
73 s_hat = (double *) MALLOC_DH (m * sizeof (double));
74 v = (double *) MALLOC_DH (m * sizeof (double));
75 p = (double *) MALLOC_DH (m * sizeof (double));
76 p_hat = (double *) MALLOC_DH (m * sizeof (double));
77 r = (double *) MALLOC_DH (m * sizeof (double));
78 r_hat = (double *) MALLOC_DH (m * sizeof (double));
79
80 /* r = b - Ax */
81 Mat_dhMatVec (A, x, s); /* s = Ax */
82 CopyVec (m, b, r); /* r = b */
83 Axpy (m, -1.0, s, r); /* r = b-Ax */
84 CopyVec (m, r, r_hat); /* r_hat = r */
85
86 /* compute stopping criteria */
87 b_iprod = InnerProd (m, b, b);
89 exit_a = atol * atol * b_iprod;
90 CHECK_V_ERROR; /* absolute stopping criteria */
91 eps = rtol * rtol * b_iprod; /* relative stoping criteria (residual reduction) */
92
93 its = 0;
94 while (1)
95 {
96 ++its;
97 rho_1 = InnerProd (m, r_hat, r);
98 if (rho_1 == 0)
99 {
100 SET_V_ERROR ("(r_hat . r) = 0; method fails");
101 }
102
103 if (its == 1)
104 {
105 CopyVec (m, r, p); /* p = r_0 */
107 }
108 else
109 {
110 beta_1 = (rho_1 / rho_2) * (alpha_1 / widget_1);
111
112 /* p_i = r_(i-1) + beta_(i-1)*( p_(i-1) - w_(i-1)*v_(i-1) ) */
113 Axpy (m, -widget_1, v, p);
115 ScaleVec (m, beta_1, p);
117 Axpy (m, 1.0, r, p);
119 }
120
121 /* solve M*p_hat = p_i */
122 Euclid_dhApply (ctx, p, p_hat);
124
125 /* v_i = A*p_hat */
126 Mat_dhMatVec (A, p_hat, v);
128
129 /* alpha_i = rho_(i-1) / (r_hat^T . v_i ) */
130 {
131 double tmp = InnerProd (m, r_hat, v);
133 alpha = rho_1 / tmp;
134 }
135
136 /* s = r_(i-1) - alpha_i*v_i */
137 CopyVec (m, r, s);
139 Axpy (m, -alpha, v, s);
141
142 /* check norm of s; if small enough:
143 * set x_i = x_(i-1) + alpha_i*p_i and stop.
144 * (Actually, we use the square of the norm)
145 */
146 s_norm = InnerProd (m, s, s);
147 if (s_norm < exit_a)
148 {
149 SET_INFO ("reached absolute stopping criteria");
150 break;
151 }
152
153 /* solve M*s_hat = s */
154 Euclid_dhApply (ctx, s, s_hat);
156
157 /* t = A*s_hat */
158 Mat_dhMatVec (A, s_hat, t);
160
161 /* w_i = (t . s)/(t . t) */
162 {
163 double tmp1, tmp2;
164 tmp1 = InnerProd (m, t, s);
166 tmp2 = InnerProd (m, t, t);
168 widget = tmp1 / tmp2;
169 }
170
171 /* x_i = x_(i-1) + alpha_i*p_hat + w_i*s_hat */
172 Axpy (m, alpha, p_hat, x);
174 Axpy (m, widget, s_hat, x);
176
177 /* r_i = s - w_i*t */
178 CopyVec (m, s, r);
180 Axpy (m, -widget, t, r);
182
183 /* check convergence; continue if necessary;
184 * for continuation it is necessary thea w != 0.
185 */
186 r_iprod = InnerProd (m, r, r);
188 if (r_iprod < eps)
189 {
190 SET_INFO ("stipulated residual reduction achieved");
191 break;
192 }
193
194 /* monitor convergence */
195 if (monitor && myid_dh == 0)
196 {
197 fprintf (stderr, "[it = %i] %e\n", its, sqrt (r_iprod / b_iprod));
198 }
199
200 /* prepare for next iteration */
201 rho_2 = rho_1;
202 widget_1 = widget;
203 alpha_1 = alpha;
204
205 if (its >= maxIts)
206 {
207 its = -its;
208 break;
209 }
210 }
211
212 *itsOUT = its;
213
214 FREE_DH (t);
215 FREE_DH (s);
216 FREE_DH (s_hat);
217 FREE_DH (v);
218 FREE_DH (p);
219 FREE_DH (p_hat);
220 FREE_DH (r);
221 FREE_DH (r_hat);
223
224#undef __FUNC__
225#define __FUNC__ "cg_euclid"
226void
227cg_euclid (Mat_dh A, Euclid_dh ctx, double *x, double *b, int *itsOUT)
228{
229 START_FUNC_DH int its, m = A->m;
230 double *p, *r, *s;
231 double alpha, beta, gamma, gamma_old, eps, bi_prod, i_prod;
232 bool monitor;
233 int maxIts = ctx->maxIts;
234 /* double atol = ctx->atol */
235 double rtol = ctx->rtol;
236
237 monitor = Parser_dhHasSwitch (parser_dh, "-monitor");
238
239 /* compute square of absolute stopping threshold */
240 /* bi_prod = <b,b> */
241 bi_prod = InnerProd (m, b, b);
243 eps = (rtol * rtol) * bi_prod;
244
245 p = (double *) MALLOC_DH (m * sizeof (double));
246 s = (double *) MALLOC_DH (m * sizeof (double));
247 r = (double *) MALLOC_DH (m * sizeof (double));
248
249 /* r = b - Ax */
250 Mat_dhMatVec (A, x, r); /* r = Ax */
252 ScaleVec (m, -1.0, r); /* r = b */
254 Axpy (m, 1.0, b, r); /* r = r + b */
256
257 /* solve Mp = r */
258 Euclid_dhApply (ctx, r, p);
260
261 /* gamma = <r,p> */
262 gamma = InnerProd (m, r, p);
264
265 its = 0;
266 while (1)
267 {
268 ++its;
269
270 /* s = A*p */
271 Mat_dhMatVec (A, p, s);
273
274 /* alpha = gamma / <s,p> */
275 {
276 double tmp = InnerProd (m, s, p);
278 alpha = gamma / tmp;
279 gamma_old = gamma;
280 }
281
282 /* x = x + alpha*p */
283 Axpy (m, alpha, p, x);
285
286 /* r = r - alpha*s */
287 Axpy (m, -alpha, s, r);
289
290 /* solve Ms = r */
291 Euclid_dhApply (ctx, r, s);
293
294 /* gamma = <r,s> */
295 gamma = InnerProd (m, r, s);
297
298 /* set i_prod for convergence test */
299 i_prod = InnerProd (m, r, r);
301
302 if (monitor && myid_dh == 0)
303 {
304 fprintf (stderr, "iter = %i rel. resid. norm: %e\n", its,
305 sqrt (i_prod / bi_prod));
306 }
307
308 /* check for convergence */
309 if (i_prod < eps)
310 break;
311
312 /* beta = gamma / gamma_old */
313 beta = gamma / gamma_old;
314
315 /* p = s + beta p */
316 ScaleVec (m, beta, p);
318 Axpy (m, 1.0, s, p);
320
321 if (its >= maxIts)
322 {
323 its = -its;
324 break;
325 }
326 }
327
328 *itsOUT = its;
329
330 FREE_DH (p);
331 FREE_DH (s);
332 FREE_DH (r);
void Euclid_dhApply(Euclid_dh ctx, double *rhs, double *lhs)
void Mat_dhMatVec(Mat_dh mat, double *x, double *b)
Definition Mat_dh.c:476
bool Parser_dhHasSwitch(Parser_dh p, char *s)
Definition Parser_dh.c:213
void ScaleVec(int n, double alpha, double *x)
Definition blas_dh.c:105
void CopyVec(int n, double *xIN, double *yOUT)
Definition blas_dh.c:91
double InnerProd(int n, double *x, double *y)
Definition blas_dh.c:118
void Axpy(int n, double alpha, double *x, double *y)
Definition blas_dh.c:77
int myid_dh
Parser_dh parser_dh
#define MALLOC_DH(s)
#define FREE_DH(p)
void bicgstab_euclid(Mat_dh A, Euclid_dh ctx, double *x, double *b, int *itsOUT)
Definition krylov_dh.c:53
void cg_euclid(Mat_dh A, Euclid_dh ctx, double *x, double *b, int *itsOUT)
Definition krylov_dh.c:227
#define SET_V_ERROR(msg)
Definition macros_dh.h:126
#define START_FUNC_DH
Definition macros_dh.h:181
#define CHECK_V_ERROR
Definition macros_dh.h:138
#define SET_INFO(msg)
Definition macros_dh.h:156
#define END_FUNC_DH
Definition macros_dh.h:187
int m
Definition Mat_dh.h:64