TBCI Numerical high perf. C++ Library: cg.h Source File

00001 
00005 /* $Id: cg.h,v 1.10.2.2 2003/07/13 22:25:07 garloff Exp $ */
00006 
00007 #ifndef TBCI_SOLVER_CG_H
00008 #define TBCI_SOLVER_CG_H
00009 
00010 #include "../basics.h"
00011 #include "precond.h"
00012 
00013 NAMESPACE_TBCI
00014 
00015 
00016 INST3(template <typename T, Matrix<T>, Vector<T> > class NN friend \
00017 int CG (const Matrix<T>&, Vector<T>&, const Vector<T>&,\
00018         const Preconditioner_Sig<T, Matrix<T> >&, unsigned int&, double&);)
00019 INST3(template <typename T, BdMatrix<T>, Vector<T> > class NN friend \
00020 int CG (const BdMatrix<T>&, Vector<T>&, const Vector<T>&,\
00021         const Preconditioner_Sig<T, BdMatrix<T> >&, unsigned int&, double&);)
00022 
00023 
00044 template < typename T, typename SysMatrix, typename SysVector >
00045 int CG (const SysMatrix &A, SysVector &x, const SysVector &b,
00046         const Preconditioner_Sig<T, SysMatrix> &M, unsigned int &max_iter, double &tol)
00047 {
00048   double resid;
00049   unsigned int dim = A.rows();
00050   SysVector p(dim), z(dim), q(dim);
00051   typename SysVector::value_type alpha, beta, rho, rho_1(0);
00052 
00053   double normb = b.fabs();
00054   if (normb == 0.0)   normb = 1.0; else normb = 1.0 / normb;
00055 
00056   SysVector r(b - A*x);
00057 
00058   
00059   if ((resid = r.fabs() * normb) <= tol) 
00060   { 
00061     tol = resid;
00062     max_iter = 0;
00063     return 0;
00064   }
00065 
00066 //Aenderung Marco:
00067   for (unsigned int i = 1; i <= max_iter; i++) 
00068   {
00069     z = M.solve(r);
00070     rho = dot(r,z);
00071     
00072     if (i == 1)
00073       p = z;
00074     else {
00075       beta = rho/ rho_1;
00076       p = z + beta * p;
00077     }
00078     
00079     q = A*p;
00080     alpha = rho/ dot(p,q);
00081     
00082     x += alpha * p;
00083     r -= alpha * q;
00084     
00085     if ((resid = r.fabs() * normb) <= tol) 
00086     {
00087       tol = resid;
00088       max_iter = i;
00089       return 0;     
00090     }
00091 
00092     rho_1 = rho;
00093   }
00094   
00095   tol = resid;
00096   return 1;
00097 }
00098 
00099   
00100 INST3(template <typename T, Matrix<T>, Vector<T> > class NN friend \
00101 int CG2 (const Matrix<T>&, Vector<T>&, const Vector<T>&,\
00102         const Preconditioner_Sig<T, Matrix<T> >&, unsigned int&, double&);)
00103 INST3(template <typename T, BdMatrix<T>, Vector<T> > class NN friend \
00104 int CG2 (const BdMatrix<T>&, Vector<T>&, const Vector<T>&,\
00105         const Preconditioner_Sig<T, BdMatrix<T> >&, unsigned int&, double&);)
00106 
00111 template < typename T, typename SysMatrix, typename SysVector>
00112 int CG2(const SysMatrix &A, SysVector &x, const SysVector &b,
00113         const Preconditioner_Sig<T, SysMatrix> &M, unsigned int &max_iter, double &tol)
00114 {
00115   double resid;
00116   unsigned int dim = A.rows();
00117   SysVector p(dim), z(dim), q(dim);
00118   typename SysVector::value_type alpha, beta, rho, rho_1(0);
00119 
00120   double normb = b.fabs();
00121   if (normb == 0.0)   normb = 1; else normb = 1 / normb;
00122 
00123   SysVector r(b - A*x);
00124 
00125   
00126   if ((resid = r.fabs() * normb) <= tol) 
00127   { 
00128     tol = resid;
00129     max_iter = 0;
00130     return 0;
00131   }
00132 
00133 //Aenderung Marco:
00134   for (unsigned int i = 1; i <= max_iter; i++) 
00135   {
00136     z = M.solve(r);
00137     rho = (r*z);
00138     
00139     if (i == 1)
00140       p = z;
00141     else {
00142       beta = rho/ rho_1;
00143       p = z + beta * p;
00144     }
00145     
00146     q = A*p;
00147     alpha = rho/ (p*q);
00148     
00149     x += alpha * p;
00150     r -= alpha * q;
00151     
00152     if ((resid = r.fabs() * normb) <= tol) 
00153     {
00154       tol = resid;
00155       max_iter = i;
00156       return 0;     
00157     }
00158 
00159     rho_1 = rho;
00160   }
00161   
00162   tol = resid;
00163   return 1;
00164 }
00165 
00166 NAMESPACE_END
00167 
00168 #endif  /* TBCI_SOLVER_CG_H */