43 #include "Euclid_dh.h"
44 #include "krylov_dh.h"
46 #include "Parser_dh.h"
51 #define __FUNC__ "bicgstab_euclid"
53 bicgstab_euclid (
Mat_dh A,
Euclid_dh ctx,
double *x,
double *b,
int *itsOUT)
55 START_FUNC_DH
int its, m = ctx->m;
57 int maxIts = ctx->maxIts;
58 double atol = ctx->atol, rtol = ctx->rtol;
61 double alpha, alpha_1,
63 widget, widget_1, rho_1, rho_2, s_norm, eps, exit_a, b_iprod, r_iprod;
66 double *t, *s, *s_hat, *v, *p, *p_hat, *r, *r_hat;
68 monitor = Parser_dhHasSwitch (parser_dh,
"-monitor");
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));
81 Mat_dhMatVec (A, x, s);
84 CopyVec (m, r, r_hat);
87 b_iprod = InnerProd (m, b, b);
89 exit_a = atol * atol * b_iprod;
91 eps = rtol * rtol * b_iprod;
97 rho_1 = InnerProd (m, r_hat, r);
100 SET_V_ERROR (
"(r_hat . r) = 0; method fails");
110 beta_1 = (rho_1 / rho_2) * (alpha_1 / widget_1);
113 Axpy (m, -widget_1, v, p);
115 ScaleVec (m, beta_1, p);
122 Euclid_dhApply (ctx, p, p_hat);
126 Mat_dhMatVec (A, p_hat, v);
131 double tmp = InnerProd (m, r_hat, v);
139 Axpy (m, -alpha, v, s);
146 s_norm = InnerProd (m, s, s);
149 SET_INFO (
"reached absolute stopping criteria");
154 Euclid_dhApply (ctx, s, s_hat);
158 Mat_dhMatVec (A, s_hat, t);
164 tmp1 = InnerProd (m, t, s);
166 tmp2 = InnerProd (m, t, t);
168 widget = tmp1 / tmp2;
172 Axpy (m, alpha, p_hat, x);
174 Axpy (m, widget, s_hat, x);
180 Axpy (m, -widget, t, r);
186 r_iprod = InnerProd (m, r, r);
190 SET_INFO (
"stipulated residual reduction achieved");
195 if (monitor && myid_dh == 0)
197 fprintf (stderr,
"[it = %i] %e\n", its, sqrt (r_iprod / b_iprod));
225 #define __FUNC__ "cg_euclid"
229 START_FUNC_DH
int its, m = A->m;
231 double alpha, beta, gamma, gamma_old, eps, bi_prod, i_prod;
233 int maxIts = ctx->maxIts;
235 double rtol = ctx->rtol;
237 monitor = Parser_dhHasSwitch (parser_dh,
"-monitor");
241 bi_prod = InnerProd (m, b, b);
243 eps = (rtol * rtol) * bi_prod;
245 p = (
double *) MALLOC_DH (m *
sizeof (
double));
246 s = (
double *) MALLOC_DH (m *
sizeof (
double));
247 r = (
double *) MALLOC_DH (m *
sizeof (
double));
250 Mat_dhMatVec (A, x, r);
252 ScaleVec (m, -1.0, r);
258 Euclid_dhApply (ctx, r, p);
262 gamma = InnerProd (m, r, p);
271 Mat_dhMatVec (A, p, s);
276 double tmp = InnerProd (m, s, p);
283 Axpy (m, alpha, p, x);
287 Axpy (m, -alpha, s, r);
291 Euclid_dhApply (ctx, r, s);
295 gamma = InnerProd (m, r, s);
299 i_prod = InnerProd (m, r, r);
302 if (monitor && myid_dh == 0)
304 fprintf (stderr,
"iter = %i rel. resid. norm: %e\n", its,
305 sqrt (i_prod / bi_prod));
313 beta = gamma / gamma_old;
316 ScaleVec (m, beta, p);