46 #include <OpenMEEG_Export.h> 57 for (
unsigned i = 0; i < m.nlin(); ++i) {
58 J(i, i) = 1.0 / m(i,i);
80 }
else if (std::abs(dy) > std::abs(dx)) {
81 double temp = dx / dy;
82 sn = 1.0 / sqrt( 1.0 + temp*temp );
85 double temp = dy / dx;
86 cs = 1.0 / sqrt( 1.0 + temp*temp );
93 double temp = cs * dx + sn * dy;
94 dy = -sn * dx + cs * dy;
103 for (
int i = k; i >= 0; i--) {
105 for (
int j = i - 1; j >= 0; j--)
106 y(j) -= h(j,i) * y(i);
108 for (
int j = 0; j <= k; j++) {
114 template<
class T,
class P>
115 unsigned GMRes(
const T& A,
const P& M,
Vector &x,
const Vector& b,
int max_iter,
double tol,
unsigned m) {
123 Vector s(m+1), cs(m+1), sn(m+1), w;
125 double normb = (M(b)).norm();
127 double beta = r.
norm();
132 if ((resid = r.
norm() / normb) <= tol) {
139 while (j <= max_iter) {
140 v[0] = r * (1.0 / beta);
144 for (i = 0; i < m && j <= max_iter; i++, j++) {
146 for (k = 0; k <= i; k++) {
150 H(i+1, i) = w.
norm();
151 v[i+1] = (w / H(i+1, i));
153 for (k = 0; k < i; k++)
160 if ((resid = std::abs(s(i+1)) / normb) < tol) {
169 Update(x, i - 1, H, s, v);
172 if ((resid = beta / normb) < tol) {
void ApplyPlaneRotation(double &dx, double &dy, double &cs, double &sn)
void GeneratePlaneRotation(double &dx, double &dy, double &cs, double &sn)
Vector operator()(const Vector &g) const
void Update(Vector &x, int k, T &h, Vector &s, Vector v[])
unsigned GMRes(const T &A, const P &M, Vector &x, const Vector &b, int max_iter, double tol, unsigned m)