Equation solvers from IML library

bicg.h

+//*****************************************************************
+// 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

bicgstab.h

+//*****************************************************************
+// 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

cg.h

+//*****************************************************************
+// 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

cgs.h

+//*****************************************************************
+// 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

cheby.h

+//*****************************************************************
+// 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

ir.h

+//*****************************************************************
+// 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
+
+

qmr.h

+//*****************************************************************
+// 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