IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
krylov_dh.c
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"
52 void
53 bicgstab_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);
88  CHECK_V_ERROR;
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 */
106  CHECK_V_ERROR;
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);
114  CHECK_V_ERROR;
115  ScaleVec (m, beta_1, p);
116  CHECK_V_ERROR;
117  Axpy (m, 1.0, r, p);
118  CHECK_V_ERROR;
119  }
120 
121  /* solve M*p_hat = p_i */
122  Euclid_dhApply (ctx, p, p_hat);
123  CHECK_V_ERROR;
124 
125  /* v_i = A*p_hat */
126  Mat_dhMatVec (A, p_hat, v);
127  CHECK_V_ERROR;
128 
129  /* alpha_i = rho_(i-1) / (r_hat^T . v_i ) */
130  {
131  double tmp = InnerProd (m, r_hat, v);
132  CHECK_V_ERROR;
133  alpha = rho_1 / tmp;
134  }
135 
136  /* s = r_(i-1) - alpha_i*v_i */
137  CopyVec (m, r, s);
138  CHECK_V_ERROR;
139  Axpy (m, -alpha, v, s);
140  CHECK_V_ERROR;
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);
155  CHECK_V_ERROR;
156 
157  /* t = A*s_hat */
158  Mat_dhMatVec (A, s_hat, t);
159  CHECK_V_ERROR;
160 
161  /* w_i = (t . s)/(t . t) */
162  {
163  double tmp1, tmp2;
164  tmp1 = InnerProd (m, t, s);
165  CHECK_V_ERROR;
166  tmp2 = InnerProd (m, t, t);
167  CHECK_V_ERROR;
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);
173  CHECK_V_ERROR;
174  Axpy (m, widget, s_hat, x);
175  CHECK_V_ERROR;
176 
177  /* r_i = s - w_i*t */
178  CopyVec (m, s, r);
179  CHECK_V_ERROR;
180  Axpy (m, -widget, t, r);
181  CHECK_V_ERROR;
182 
183  /* check convergence; continue if necessary;
184  * for continuation it is necessary thea w != 0.
185  */
186  r_iprod = InnerProd (m, r, r);
187  CHECK_V_ERROR;
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);
222 END_FUNC_DH}
223 
224 #undef __FUNC__
225 #define __FUNC__ "cg_euclid"
226 void
227 cg_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);
242  CHECK_V_ERROR;
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 */
251  CHECK_V_ERROR;
252  ScaleVec (m, -1.0, r); /* r = b */
253  CHECK_V_ERROR;
254  Axpy (m, 1.0, b, r); /* r = r + b */
255  CHECK_V_ERROR;
256 
257  /* solve Mp = r */
258  Euclid_dhApply (ctx, r, p);
259  CHECK_V_ERROR;
260 
261  /* gamma = <r,p> */
262  gamma = InnerProd (m, r, p);
263  CHECK_V_ERROR;
264 
265  its = 0;
266  while (1)
267  {
268  ++its;
269 
270  /* s = A*p */
271  Mat_dhMatVec (A, p, s);
272  CHECK_V_ERROR;
273 
274  /* alpha = gamma / <s,p> */
275  {
276  double tmp = InnerProd (m, s, p);
277  CHECK_V_ERROR;
278  alpha = gamma / tmp;
279  gamma_old = gamma;
280  }
281 
282  /* x = x + alpha*p */
283  Axpy (m, alpha, p, x);
284  CHECK_V_ERROR;
285 
286  /* r = r - alpha*s */
287  Axpy (m, -alpha, s, r);
288  CHECK_V_ERROR;
289 
290  /* solve Ms = r */
291  Euclid_dhApply (ctx, r, s);
292  CHECK_V_ERROR;
293 
294  /* gamma = <r,s> */
295  gamma = InnerProd (m, r, s);
296  CHECK_V_ERROR;
297 
298  /* set i_prod for convergence test */
299  i_prod = InnerProd (m, r, r);
300  CHECK_V_ERROR;
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);
317  CHECK_V_ERROR;
318  Axpy (m, 1.0, s, p);
319  CHECK_V_ERROR;
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);
333 END_FUNC_DH}
Definition: Mat_dh.h:62