14#include <OpenMEEG_Export.h>
25 for (
unsigned i = 0; i < m.nlin(); ++i) {
26 J(i, i) = 1.0 / m(i,i);
48 }
else if (std::abs(dy) > std::abs(dx)) {
49 double temp = dx / dy;
50 sn = 1.0 / sqrt( 1.0 + temp*temp );
53 double temp = dy / dx;
54 cs = 1.0 / sqrt( 1.0 + temp*temp );
61 double temp = cs * dx + sn * dy;
62 dy = -sn * dx + cs * dy;
71 for (
int i = k; i >= 0; i--) {
73 for (
int j = i - 1; j >= 0; j--)
74 y(j) -= h(j,i) * y(i);
76 for (
int j = 0; j <= k; j++) {
82 template<
class T,
class P>
83 unsigned GMRes(
const T& A,
const P& M,
Vector &x,
const Vector& b,
int max_iter,
double tol,
unsigned m) {
91 Vector s(m+1), cs(m+1), sn(m+1), w;
93 double normb = (M(b)).norm();
95 double beta = r.
norm();
100 if ((resid = r.
norm() / normb) <= tol) {
107 while (j <= max_iter) {
108 v[0] = r * (1.0 / beta);
112 for (i = 0; i < m && j <= max_iter; i++, j++) {
114 for (k = 0; k <= i; k++) {
118 H(i+1, i) = w.
norm();
119 v[i+1] = (w / H(i+1, i));
121 for (k = 0; k < i; k++)
128 if ((resid = std::abs(s(i+1)) / normb) < tol) {
137 Update(x, i - 1, H, s, v);
140 if ((resid = beta / normb) < tol) {
Vector operator()(const Vector &g) const
Matrix class Matrix class.
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
unsigned GMRes(const T &A, const P &M, Vector &x, const Vector &b, int max_iter, double tol, unsigned m)
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
void Update(Vector &x, int k, T &h, Vector &s, Vector v[])