From 2adc07a3e1ba8d4346557f8761901835e82715ac Mon Sep 17 00:00:00 2001
From: Jesper Andersson <jesper@fmrib.ox.ac.uk>
Date: Mon, 4 Jun 2007 13:02:33 +0000
Subject: [PATCH] Equation solvers from IML library
---
bicg.h | 93 +++++++++++++++++++++++++++++
bicgstab.h | 99 +++++++++++++++++++++++++++++++
cg.h | 85 +++++++++++++++++++++++++++
cgs.h | 90 +++++++++++++++++++++++++++++
cheby.h | 88 ++++++++++++++++++++++++++++
ir.h | 69 ++++++++++++++++++++++
qmr.h | 167 +++++++++++++++++++++++++++++++++++++++++++++++++++++
7 files changed, 691 insertions(+)
create mode 100644 bicg.h
create mode 100644 bicgstab.h
create mode 100644 cg.h
create mode 100644 cgs.h
create mode 100644 cheby.h
create mode 100644 ir.h
create mode 100644 qmr.h
diff --git a/bicg.h b/bicg.h
new file mode 100644
index 0000000..10f2d37
--- /dev/null
+++ b/bicg.h
@@ -0,0 +1,93 @@
+//*****************************************************************
+// Iterative template routine -- BiCG
+//
+// BiCG solves the unsymmetric linear system Ax = b
+// using the Preconditioned BiConjugate Gradient method
+//
+// BiCG follows the algorithm described on p. 22 of the
+// SIAM Templates book.
+//
+// The return value indicates convergence within max_iter (input)
+// iterations (0), or no convergence within max_iter iterations (1).
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax = b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef bicg_h
+#define bicg_h
+
+namespace MISCMATHS {
+
+
+template < class Matrix, class Vector, class Preconditioner, class Real >
+int
+BiCG(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &max_iter, Real &tol)
+{
+ Real resid;
+ Vector rho_1(1), rho_2(1), alpha(1), beta(1);
+ Vector z, ztilde, p, ptilde, q, qtilde;
+
+ Real normb = b.NormFrobenius();
+ Vector r = b - A * x;
+ Vector rtilde = r;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ for (int i = 1; i <= max_iter; i++) {
+ z = M.solve(r);
+ ztilde = M.trans_solve(rtilde);
+ rho_1(1) = DotProduct(z, rtilde);
+ if (rho_1(1) == 0) {
+ tol = r.NormFrobenius() / normb;
+ max_iter = i;
+ return 2;
+ }
+ if (i == 1) {
+ p = z;
+ ptilde = ztilde;
+ } else {
+ beta(1) = rho_1(1) / rho_2(1);
+ p = z + beta(1) * p;
+ ptilde = ztilde + beta(1) * ptilde;
+ }
+ q = A * p;
+ qtilde = A.trans_mult(ptilde);
+ alpha(1) = rho_1(1) / DotProduct(ptilde, q);
+ x += alpha(1) * p;
+ r -= alpha(1) * q;
+ rtilde -= alpha(1) * qtilde;
+
+ rho_2(1) = rho_1(1);
+ if ((resid = r.NormFrobenius() / normb) < tol) {
+ tol = resid;
+ max_iter = i;
+ return 0;
+ }
+ }
+
+ tol = resid;
+ return 1;
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef cg_h
diff --git a/bicgstab.h b/bicgstab.h
new file mode 100644
index 0000000..fbaac51
--- /dev/null
+++ b/bicgstab.h
@@ -0,0 +1,99 @@
+//*****************************************************************
+// Iterative template routine -- BiCGSTAB
+//
+// BiCGSTAB solves the unsymmetric linear system Ax = b
+// using the Preconditioned BiConjugate Gradient Stabilized method
+//
+// BiCGSTAB follows the algorithm described on p. 27 of the
+// SIAM Templates book.
+//
+// The return value indicates convergence within max_iter (input)
+// iterations (0), or no convergence within max_iter iterations (1).
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax = b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef bicgstab_h
+#define bicgstab_h
+
+namespace MISCMATHS {
+
+template < class Matrix, class Vector, class Preconditioner, class Real >
+int
+BiCGSTAB(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &max_iter, Real &tol)
+{
+ Real resid;
+ Vector rho_1(1), rho_2(1), alpha(1), beta(1), omega(1);
+ Vector p, phat, s, shat, t, v;
+
+ Real normb = b.NormFrobenius();
+ Vector r = b - A * x;
+ Vector rtilde = r;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ for (int i = 1; i <= max_iter; i++) {
+ rho_1(1) = DotProduct(rtilde, r);
+ if (rho_1(1) == 0) {
+ tol = r.NormFrobenius() / normb;
+ return 2;
+ }
+ if (i == 1)
+ p = r;
+ else {
+ beta(1) = (rho_1(1)/rho_2(1)) * (alpha(1)/omega(1));
+ p = r + beta(1) * (p - omega(1) * v);
+ }
+ phat = M.solve(p);
+ v = A * phat;
+ alpha(1) = rho_1(1) / DotProduct(rtilde, v);
+ s = r - alpha(1) * v;
+ if ((resid = s.NormFrobenius()/normb) < tol) {
+ x += alpha(1) * phat;
+ tol = resid;
+ return 0;
+ }
+ shat = M.solve(s);
+ t = A * shat;
+ omega = DotProduct(t,s) / DotProduct(t,t);
+ x += alpha(1) * phat + omega(1) * shat;
+ r = s - omega(1) * t;
+
+ rho_2(1) = rho_1(1);
+ if ((resid = r.NormFrobenius() / normb) < tol) {
+ tol = resid;
+ max_iter = i;
+ return 0;
+ }
+ if (omega(1) == 0) {
+ tol = r.NormFrobenius() / normb;
+ return 3;
+ }
+ }
+
+ tol = resid;
+ return 1;
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef bicgstab_h
diff --git a/cg.h b/cg.h
new file mode 100644
index 0000000..42094e1
--- /dev/null
+++ b/cg.h
@@ -0,0 +1,85 @@
+//*****************************************************************
+// Iterative template routine -- CG
+//
+// CG solves the symmetric positive definite linear
+// system Ax=b using the Conjugate Gradient method.
+//
+// CG follows the algorithm described on p. 15 in the
+// SIAM Templates book.
+//
+// The return value indicates convergence within max_iter (input)
+// iterations (0), or no convergence within max_iter iterations (1).
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax = b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef cg_h
+#define cg_h
+
+namespace MISCMATHS {
+
+template < class Matrix, class Vector, class Preconditioner, class Real >
+int
+CG(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &max_iter, Real &tol)
+{
+ Real resid;
+ Vector p, z, q;
+ Vector alpha(1), beta(1), rho(1), rho_1(1);
+
+ Real normb = b.NormFrobenius();
+ Vector r = b - A*x;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ for (int i = 1; i <= max_iter; i++) {
+ z = M.solve(r);
+ rho(1) = DotProduct(r, z);
+
+ if (i == 1)
+ p = z;
+ else {
+ beta(1) = rho(1) / rho_1(1);
+ p = z + beta(1) * p;
+ }
+
+ q = A*p;
+ alpha(1) = rho(1) / DotProduct(p, q);
+
+ x += alpha(1) * p;
+ r -= alpha(1) * q;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = i;
+ return 0;
+ }
+
+ rho_1(1) = rho(1);
+ }
+
+ tol = resid;
+ return 1;
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef cg_h
diff --git a/cgs.h b/cgs.h
new file mode 100644
index 0000000..1765336
--- /dev/null
+++ b/cgs.h
@@ -0,0 +1,90 @@
+//*****************************************************************
+// Iterative template routine -- CGS
+//
+// CGS solves the unsymmetric linear system Ax = b
+// using the Conjugate Gradient Squared method
+//
+// CGS follows the algorithm described on p. 26 of the
+// SIAM Templates book.
+//
+// The return value indicates convergence within max_iter (input)
+// iterations (0), or no convergence within max_iter iterations (1).
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax = b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef cgs_h
+#define cgs_h
+
+namespace MISCMATHS {
+
+template < class Matrix, class Vector, class Preconditioner, class Real >
+int
+CGS(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &max_iter, Real &tol)
+{
+ Real resid;
+ Vector rho_1(1), rho_2(1), alpha(1), beta(1);
+ Vector p, phat, q, qhat, vhat, u, uhat;
+
+ Real normb = b.NormFrobenius();
+ Vector r = b - A*x;
+ Vector rtilde = r;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ for (int i = 1; i <= max_iter; i++) {
+ rho_1(1) = DotProduct(rtilde, r);
+ if (rho_1(1) == 0) {
+ tol = r.NormFrobenius() / normb;
+ return 2;
+ }
+ if (i == 1) {
+ u = r;
+ p = u;
+ } else {
+ beta(1) = rho_1(1) / rho_2(1);
+ u = r + beta(1) * q;
+ p = u + beta(1) * (q + beta(1) * p);
+ }
+ phat = M.solve(p);
+ vhat = A*phat;
+ alpha(1) = rho_1(1) / DotProduct(rtilde, vhat);
+ q = u - alpha(1) * vhat;
+ uhat = M.solve(u + q);
+ x += alpha(1) * uhat;
+ qhat = A * uhat;
+ r -= alpha(1) * qhat;
+ rho_2(1) = rho_1(1);
+ if ((resid = r.NormFrobenius() / normb) < tol) {
+ tol = resid;
+ max_iter = i;
+ return 0;
+ }
+ }
+
+ tol = resid;
+ return 1;
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef cgs_h
diff --git a/cheby.h b/cheby.h
new file mode 100644
index 0000000..67be8d5
--- /dev/null
+++ b/cheby.h
@@ -0,0 +1,88 @@
+//*****************************************************************
+// Iterative template routine -- CHEBY
+//
+// CHEBY solves the symmetric positive definite linear
+// system Ax = b using the Preconditioned Chebyshev Method
+//
+// CHEBY follows the algorithm described on p. 30 of the
+// SIAM Templates book.
+//
+// The return value indicates convergence within max_iter (input)
+// iterations (0), or no convergence within max_iter iterations (1).
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax = b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef cheby_h
+#define cheby_h
+
+namespace MISCMATHS {
+
+template < class Matrix, class Vector, class Preconditioner, class Real,
+ class Type >
+int
+CHEBY(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &max_iter, Real &tol,
+ Type eigmin, Type eigmax)
+{
+ Real resid;
+ Type alpha, beta, c, d;
+ Vector p, q, z;
+
+ Real normb = b.NormFrobenius();
+ Vector r = b - A * x;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ c = (eigmax - eigmin) / 2.0;
+ d = (eigmax + eigmin) / 2.0;
+
+ for (int i = 1; i <= max_iter; i++) {
+ z = M.solve(r); // apply preconditioner
+
+ if (i == 1) {
+ p = z;
+ alpha = 2.0 / d;
+ } else {
+ beta = c * alpha / 2.0; // calculate new beta
+ beta = beta * beta;
+ alpha = 1.0 / (d - beta); // calculate new alpha
+ p = z + beta * p; // update search direction
+ }
+
+ q = A * p;
+ x += alpha * p; // update approximation vector
+ r -= alpha * q; // compute residual
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = i;
+ return 0; // convergence
+ }
+ }
+
+ tol = resid;
+ return 1; // no convergence
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef cheby_h
diff --git a/ir.h b/ir.h
new file mode 100644
index 0000000..c6f6ee2
--- /dev/null
+++ b/ir.h
@@ -0,0 +1,69 @@
+//*****************************************************************
+// Iterative template routine -- Preconditioned Richardson
+//
+// IR solves the unsymmetric linear system Ax = b using
+// Iterative Refinement (preconditioned Richardson iteration).
+//
+// The return value indicates convergence within max_iter (input)
+// iterations (0), or no convergence within max_iter iterations (1).
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax = b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef ir_h
+#define ir_h
+
+namespace MISCMATHS {
+
+template < class Matrix, class Vector, class Preconditioner, class Real >
+int
+IR(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &max_iter, Real &tol)
+{
+ Real resid;
+ Vector z;
+
+ Real normb = b.NormFrobenius();
+ Vector r = b - A*x;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ for (int i = 1; i <= max_iter; i++) {
+ z = M.solve(r);
+ x += z;
+ r = b - A * x;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = i;
+ return 0;
+ }
+ }
+
+ tol = resid;
+ return 1;
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef ir_h
+
+
diff --git a/qmr.h b/qmr.h
new file mode 100644
index 0000000..608d3f3
--- /dev/null
+++ b/qmr.h
@@ -0,0 +1,167 @@
+//*****************************************************************
+// Iterative template routine -- QMR
+//
+// QMR.h solves the unsymmetric linear system Ax = b using the
+// Quasi-Minimal Residual method following the algorithm as described
+// on p. 24 in the SIAM Templates book.
+//
+// -------------------------------------------------------------
+// return value indicates
+// ------------ ---------------------
+// 0 convergence within max_iter iterations
+// 1 no convergence after max_iter iterations
+// breakdown in:
+// 2 rho
+// 3 beta
+// 4 gamma
+// 5 delta
+// 6 ep
+// 7 xi
+// -------------------------------------------------------------
+//
+// Upon successful return, output arguments have the following values:
+//
+// x -- approximate solution to Ax=b
+// max_iter -- the number of iterations performed before the
+// tolerance was reached
+// tol -- the residual after the final iteration
+//
+//*****************************************************************
+//
+// Slightly modified version of IML++ template. See ReadMe file.
+//
+// Jesper Andersson
+//
+
+#ifndef qmr_h
+#define qmr_h
+
+#include <math.h>
+
+namespace MISCMATHS {
+
+template < class Matrix, class Vector, class Preconditioner1,
+ class Preconditioner2, class Real >
+int
+QMR(const Matrix &A, Vector &x, const Vector &b, const Preconditioner1 &M1,
+ const Preconditioner2 &M2, int &max_iter, Real &tol)
+{
+ Real resid;
+
+ Vector rho(1), rho_1(1), xi(1), gamma(1), gamma_1(1), theta(1), theta_1(1);
+ Vector eta(1), delta(1), ep(1), beta(1);
+
+ Vector r, v_tld, y, w_tld, z;
+ Vector v, w, y_tld, z_tld;
+ Vector p, q, p_tld, d, s;
+
+ Real normb = b.NormFrobenius();
+
+ r = b - A * x;
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ v_tld = r;
+ y = M1.solve(v_tld);
+ rho(1) = y.NormFrobenius();
+
+ w_tld = r;
+ z = M2.trans_solve(w_tld);
+ xi(1) = z.NormFrobenius();
+
+ gamma(1) = 1.0;
+ eta(1) = -1.0;
+ theta(1) = 0.0;
+
+ for (int i = 1; i <= max_iter; i++) {
+
+ if (rho(1) == 0.0)
+ return 2; // return on breakdown
+
+ if (xi(1) == 0.0)
+ return 7; // return on breakdown
+
+ v = (1. / rho(1)) * v_tld;
+ y = (1. / rho(1)) * y;
+
+ w = (1. / xi(1)) * w_tld;
+ z = (1. / xi(1)) * z;
+
+ delta(1) = DotProduct(z, y);
+ if (delta(1) == 0.0)
+ return 5; // return on breakdown
+
+ y_tld = M2.solve(y); // apply preconditioners
+ z_tld = M1.trans_solve(z);
+
+ if (i > 1) {
+ p = y_tld - (xi(1) * delta(1) / ep(1)) * p;
+ q = z_tld - (rho(1) * delta(1) / ep(1)) * q;
+ } else {
+ p = y_tld;
+ q = z_tld;
+ }
+
+ p_tld = A * p;
+ ep(1) = DotProduct(q, p_tld);
+ if (ep(1) == 0.0)
+ return 6; // return on breakdown
+
+ beta(1) = ep(1) / delta(1);
+ if (beta(1) == 0.0)
+ return 3; // return on breakdown
+
+ v_tld = p_tld - beta(1) * v;
+ y = M1.solve(v_tld);
+
+ rho_1(1) = rho(1);
+ rho(1) = y.NormFrobenius();
+ w_tld = A.trans_mult(q) - beta(1) * w;
+ z = M2.trans_solve(w_tld);
+
+ xi(1) = z.NormFrobenius();
+
+ gamma_1(1) = gamma(1);
+ theta_1(1) = theta(1);
+
+ theta(1) = rho(1) / (gamma_1(1) * beta(1));
+ gamma(1) = 1.0 / sqrt(1.0 + theta(1) * theta(1));
+
+ if (gamma(1) == 0.0)
+ return 4; // return on breakdown
+
+ eta(1) = -eta(1) * rho_1(1) * gamma(1) * gamma(1) /
+ (beta(1) * gamma_1(1) * gamma_1(1));
+
+ if (i > 1) {
+ d = eta(1) * p + (theta_1(1) * theta_1(1) * gamma(1) * gamma(1)) * d;
+ s = eta(1) * p_tld + (theta_1(1) * theta_1(1) * gamma(1) * gamma(1)) * s;
+ } else {
+ d = eta(1) * p;
+ s = eta(1) * p_tld;
+ }
+
+ x += d; // update approximation vector
+ r -= s; // compute residual
+
+ if ((resid = r.NormFrobenius() / normb) <= tol) {
+ tol = resid;
+ max_iter = i;
+ return 0;
+ }
+ }
+
+ tol = resid;
+ return 1; // no convergence
+}
+
+} // End namespace MISCMATHS
+
+#endif // End #ifndef qmr_h
--
GitLab