AndersonMixer.C 1.67 KB
Newer Older
Francois Gygi committed
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// AndersonMixer.C
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: AndersonMixer.C,v 1.4 2007-03-17 01:14:00 fgygi Exp $
Francois Gygi committed
7 8 9

#include "AndersonMixer.h"
#include "blas.h"
10
using namespace std;
Francois Gygi committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

////////////////////////////////////////////////////////////////////////////////
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);
      
38 39 40 41 42 43 44
      if ( pctxt_ != 0 )
      {
        double work[2] = { a, b };
        pctxt_->dsum(2,1,work,2);
        a = work[0];
        b = work[1];
      }
Francois Gygi committed
45 46 47 48 49 50 51 52 53

      if ( b != 0.0 )
        *theta = - a / b;
      else
        *theta = 0.0;
      
      // test if residual has increased
      if ( *theta <= -1.0 )
      {
54
        *theta = theta_nc_;
Francois Gygi committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
      }
      
      *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;
}