Anderson mixer class

git-svn-id: http://qboxcode.org/svn/qb/trunk@306 cba15fb0-1239-40c8-b417-11db7ca47a34

src/AndersonMixer.C

+////////////////////////////////////////////////////////////////////////////////
+//
+// AndersonMixer.C
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: AndersonMixer.C,v 1.1 2004-12-02 22:24:16 fgygi Exp $
+
+#include "AndersonMixer.h"
+#include "blas.h"
+
+////////////////////////////////////////////////////////////////////////////////
+void AndersonMixer::restart(void)
+{
+  extrapolate_ = false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void AndersonMixer::update(const double* f, double* theta, double* fbar)
+{
+    *theta = 0.0;
+    if ( extrapolate_ )
+    {
+      int ione=1;
+      valarray<double> tmp0(n_);
+
+      // compute theta = - a / b
+      // tmp0 = delta_F = f - flast
+      for ( int i = 0; i < n_; i++ )
+        tmp0[i] = f[i] - flast_[i];
+      
+      // a = delta_F * F
+      double a = ddot(&n_, &tmp0[0], &ione, f, &ione);
+        
+      // b = delta_F * delta_F
+      double b = ddot(&n_, &tmp0[0], &ione, &tmp0[0], &ione);
+      
+      double work[2] = { a, b };
+      ctxt_.dsum(2,1,work,2);
+      a = work[0];
+      b = work[1];
+
+      if ( b != 0.0 )
+        *theta = - a / b;
+      else
+        *theta = 0.0;
+      
+      // test if residual has increased
+      if ( *theta <= -1.0 )
+      {
+        *theta = 1.0;
+      }
+      
+      *theta = min(theta_max_,*theta);
+    }
+    
+    // fbar = f + theta * ( f - flast )
+    // flast_ = f_
+    for ( int i = 0; i < n_; i++ )
+    {
+      const double ftmp = f[i];
+      fbar[i] = ftmp + *theta * ( ftmp - flast_[i] );
+      flast_[i] = ftmp;
+    } 
+    extrapolate_ = true;
+}

src/AndersonMixer.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// AndersonMixer.h
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: AndersonMixer.h,v 1.1 2004-12-02 22:24:16 fgygi Exp $
+
+#ifndef ANDERSONMIXER_H
+#define ANDERSONMIXER_H
+
+#include <valarray>
+#include <cassert>
+using namespace std;
+
+#include "Context.h"
+
+class AndersonMixer
+{
+  int     n_;              // size of vectors
+  const   Context ctxt_;
+  double  theta_max_;
+ 
+  valarray<double> flast_; // last residual
+  bool extrapolate_;       // state variable
+
+  public:
+    
+  AndersonMixer(const int n, const Context& ctxt) :
+    n_(n), ctxt_(ctxt), extrapolate_(false), theta_max_(2.0)
+  {
+    assert( n > 0 );
+    flast_.resize(n);
+  }
+
+  void update(const double* f, double* theta, double* fbar);
+  void restart(void);
+  void set_theta_max(double theta_max) { theta_max_ = theta_max; }
+  double theta_max(void) const { return theta_max_; }
+};
+#endif

src/testAndersonMixer.C

+////////////////////////////////////////////////////////////////////////////////
+//
+// testAndersonMixer.C
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: testAndersonMixer.C,v 1.1 2004-12-02 22:24:16 fgygi Exp $
+
+
+#include <iostream>
+using namespace std;
+
+#include "Context.h"
+#include "AndersonMixer.h"
+
+int main(int argc, char** argv)
+{
+#if USE_MPI
+  MPI_Init(&argc,&argv);
+#endif
+  {
+    Context ctxt;
+ 
+    char processor_name[MPI_MAX_PROCESSOR_NAME];
+    int namelen;
+    PMPI_Get_processor_name(processor_name,&namelen);
+    cout << " Process " << ctxt.mype() << " on " << processor_name << endl;
+
+    vector<double> x,f,xlast,fbar;
+    double theta;
+    x.resize(3);
+    xlast.resize(3);
+    f.resize(3);
+    fbar.resize(3);
+    
+    AndersonMixer mixer(x.size(),ctxt);
+    
+    x[0] = 1.0 + 0.2 * ctxt.mype();
+    x[1] = 2.0 + 0.2 * ctxt.mype();
+    x[2] = 3.0 + 0.2 * ctxt.mype();
+    
+    for ( int iter = 0; iter < 100; iter++ )
+    {
+      f[0] = - 2.0 * x[0] - 3.0 * x[0]*x[0]*x[0];
+      f[1] = - 1.0 * x[1] - 1.5 * x[1]*x[1]*x[1];
+      f[2] = - 4.0 * x[2];
+      
+      double resnorm = 0.0;
+      for ( int i = 0; i < x.size(); i++ )
+      {
+        resnorm += f[i]*f[i];
+        f[i] *= 0.1;
+      }
+      ctxt.dsum(1,1,&resnorm,1);
+
+      cout << ctxt.mype() << ": x: " << x[0] << " " << x[1] << " " << x[2] 
+           << " resnorm: " << resnorm << endl;
+    
+      mixer.update(&f[0],&theta,&fbar[0]);
+      
+      for ( int i = 0; i < x.size(); i++ )
+      {
+        // xbar = x + theta * ( x - xlast )
+        // x = xbar + fbar
+        double xtmp = x[i];
+        x[i] += theta * ( x[i] - xlast[i] ) + fbar[i];
+        xlast[i] = xtmp;
+      }
+    }
+  }
+#if USE_MPI
+  MPI_Finalize();
+#endif
+  return 0;
+}