cgmatrix.h

Go to the documentation of this file.
00001 #ifndef __cgmatrix_h
00002 #define __cgmatrix_h
00003 
00004 
00005 namespace Gascoigne
00006 {
00007 template <class MT,class VT,class MEM, class INFO>
00008 inline int CgMatrix(const MT& M, VT& x, const VT& b, MEM& mem, INFO& info)
00009 {
00010   
00011   VT& g  = mem[0];
00012   VT& h  = mem[1];
00013   VT& d  = mem[2];
00014   VT& Ad = mem[3];
00015 
00016   int    it=0,reached = 0;
00017   double gh,alpha,beta;
00018   double     res;
00019 
00020   M.vmulteq(g,x);
00021   g.sequ(1.,-1.,b);
00022   res = g.norm();
00023 
00024   if(res==0.)      return 0;
00025 
00026   d = g;
00027   //M.precondition(d,g);
00028 
00029   gh  =  g*d;
00030   gh *=  -1.;
00031 
00032   for(it=0;!reached;it++)
00033   {
00034     M.vmulteq(Ad,d);
00035 
00036     alpha = d*Ad;
00037     alpha = gh/alpha;
00038 
00039     g.add(alpha,Ad);
00040     x.add(alpha,d );
00041     res = g.norm();
00042 
00043     reached = info.check_residual(it,res);
00044     if (reached) continue;
00045     
00046     h = g;
00047     //M.precondition(h,g);
00048     
00049     beta = gh;
00050     gh   = g*h;
00051     beta = gh/beta;
00052     
00053     d.sequ(beta,-1.,h);
00054   }
00055   if (reached<0) return 1;
00056   return 0;
00057 }
00058 }
00059 
00060 #endif

Generated on Thu Sep 14 10:34:36 2006 for Gascoigne by  doxygen 1.4.7