Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
solver_CG_alt.cpp
Go to the documentation of this file.
1 
15 #include "solver_CG.h"
16 
17 #include "communicator.h"
18 #include "bridgeIO.h"
19 using Bridge::vout;
20 
21 #include "fopr.h"
22 
23 
24 //====================================================================
25 void Solver_CG::solve(Field& xq, const Field& b,
26  int& Nconv, double& diff)
27 {
28  // vout.general(m_vl, "CG solver start\n");
29 
30  double snorm = 1.0 / b.norm2();
31 
32  reset_field(b);
33  s = b;
34 
35  Nconv = -1;
36  double rr;
37 
38  solve_init(rr);
39  //vout.general(m_vl, " init: %22.15e\n",rr*snorm);
40 
41  for (int it = 0; it < Niter; it++) {
42  solve_step(rr);
43  // vout.general(m_vl, "%6d %22.15e\n",it,rr*snorm);
44 
45  if (rr * snorm < enorm) {
46  Nconv = it;
47  break;
48  }
49  }
50  if (Nconv == -1) { cout << "Not converged." << __FILE__ << "(" << __LINE__ << ")" << endl; abort(); }
51 
52  //p = opr->mult(x);
53  opr->mult(p, x);
54  p -= b;
55  diff = p.norm();
56 
57  xq = x;
58 }
59 
60 
61 //====================================================================
62 void Solver_CG::reset_field(const Field& b)
63 {
64  int Nin = b.nin();
65  int Nvol = b.nvol();
66  int Nex = b.nex();
67 
68  if ((s.nin() != Nin) || (s.nvol() != Nvol) || (s.nex() != Nex)) {
69  s.reset(Nin, Nvol, Nex);
70  r.reset(Nin, Nvol, Nex);
71  x.reset(Nin, Nvol, Nex);
72  p.reset(Nin, Nvol, Nex);
73  v.reset(Nin, Nvol, Nex);
74 
75  vout.general(m_vl, " Solver_CG: field size reset.\n");
76  }
77 }
78 
79 
80 //====================================================================
81 void Solver_CG::solve_init(double& rr)
82 {
83  r = s;
84  x = s;
85  //s = opr->mult(x);
86  opr->mult(s, x);
87  r -= s;
88  p = r;
89  rr = r * r;
90 }
91 
92 
93 //====================================================================
94 void Solver_CG::solve_step(double& rr)
95 {
96  //s = opr->mult(p);
97  opr->mult(s, p);
98 
99  double pap = p * s;
100  double rrp = rr;
101  double cr = rrp / pap;
102 
103  v = p;
104  v *= cr;
105  x += v;
106  // x += cr*p;
107 
108  s *= cr;
109  r -= s;
110 
111  rr = r * r;
112  p *= rr / rrp;
113  p += r;
114 }
115 
116 
117 //====================================================================
118 //============================================================END=====