AndersonMixer.C 7.44 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20
// AndersonMixer.C
//
////////////////////////////////////////////////////////////////////////////////

#include "AndersonMixer.h"
#include "blas.h"
21
#include <iostream>
22 23 24 25 26
#if USE_MPI
#include <mpi.h>
#else
typedef int MPI_Comm;
#endif
27
using namespace std;
Francois Gygi committed
28

29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
////////////////////////////////////////////////////////////////////////////////
AndersonMixer::AndersonMixer(const int m, const int nmax,
  const MPI_Comm* const pcomm) : m_(m), nmax_(nmax), pcomm_(pcomm)
{
#if USE_MPI
  if ( pcomm_ != 0 )
  {
    MPI_Comm_rank(*pcomm_,&mype_);
    MPI_Comm_size(*pcomm_,&npes_);
  }
#endif

  assert( nmax >= 0 );
  x_.resize(nmax_+1);
  f_.resize(nmax_+1);
  for ( int n = 0; n < nmax_+1; n++ )
  {
    x_[n].resize(m_);
    f_[n].resize(m_);
  }
  restart();
}

Francois Gygi committed
52 53 54
////////////////////////////////////////////////////////////////////////////////
void AndersonMixer::restart(void)
{
55 56
  n_ = -1;
  k_ = -1;
Francois Gygi committed
57 58 59
}

////////////////////////////////////////////////////////////////////////////////
60
void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
Francois Gygi committed
61
{
62 63 64 65 66 67
  // update:
  // input: x, f
  // output: xbar, fbar
  //
  // Computes the pair (xbar,fbar) using pairs (x,f) used
  // in previous updates, according to the Anderson algorithm.
Francois Gygi committed
68

69 70 71 72 73 74
  // increment index of current vector
  k_ = ( k_ + 1 ) % ( nmax_ + 1 );
  // increment current number of vectors
  if ( n_ < nmax_ ) n_++;

  // save vectors
Francois Gygi committed
75

76 77 78 79 80 81
  for ( int i = 0; i < m_; i++ )
  {
    x_[k_][i] = x[i];
    f_[k_][i] = f[i];
  }

82 83
  valarray<double> a,atmp;
  valarray<double> b,btmp;
84 85 86 87 88 89
  valarray<double> theta;
  if ( n_ > 0 )
  {
    // compute matrix A = F^T F and rhs b = F^T f
    // compute the lower part of A only (i>=j)
    a.resize(n_*n_);
90
    atmp.resize(n_*n_);
91
    b.resize(n_);
92
    btmp.resize(n_);
93 94
    theta.resize(n_);
    for ( int i = 0; i < n_; i++ )
Francois Gygi committed
95
    {
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
      const int kmi = ( k_ - i + nmax_ ) % ( nmax_ + 1 );
      assert(kmi>=0);
      assert(kmi<nmax_+1);
      // cout << "k=" << k_ << " i=" << i << " kmi=" << kmi << endl;
      for ( int j = 0; j <= i; j++ )
      {
        const int kmj = ( k_ - j + nmax_ ) % ( nmax_ + 1 );
        assert(kmj>=0);
        assert(kmj<nmax_+1);
        // cout << "k=" << k_ << " j=" << j << " kmj=" << kmj << endl;
        double sum = 0.0;
        for ( int l = 0; l < m_; l++ )
          sum += (f_[k_][l] - f_[kmi][l]) * (f_[k_][l] - f_[kmj][l]);
        a[i+j*n_] = sum;
      }
      double bsum = 0.0;
      for ( int l = 0; l < m_; l++ )
        bsum += ( f_[k_][l] - f_[kmi][l] ) * f_[k_][l];
      b[i] = bsum;
    }
116

117 118
#if USE_MPI
    if ( pcomm_ != 0 )
119
    {
120 121 122 123
      MPI_Allreduce(&a[0],&atmp[0],n_*n_,MPI_DOUBLE,MPI_SUM,*pcomm_);
      a = atmp;
      MPI_Allreduce(&b[0],&btmp[0],n_,MPI_DOUBLE,MPI_SUM,*pcomm_);
      b = btmp;
124
    }
Francois Gygi committed
125 126 127 128
#endif

    // solve the linear system a * theta = b
    // solve on task 0 and bcast result
129
    if ( pcomm_ == 0 || mype_ == 0 )
130
    {
Francois Gygi committed
131 132
      const bool diag = false;
      if ( diag )
133
      {
Francois Gygi committed
134 135 136 137 138 139 140 141 142
        // solve the linear system using eigenvalues and eigenvectors
        // compute eigenvalues of a
        char jobz = 'V';
        char uplo = 'L';
        valarray<double> w(n_);
        int lwork = 3*n_;
        valarray<double> work(lwork);
        int info;
        dsyev(&jobz,&uplo,&n_,&a[0],&n_,&w[0],&work[0],&lwork,&info);
143 144 145 146 147

        cout << "AndersonMixer: eigenvalues: ";
        for ( int i = 0; i < n_; i++ )
          cout << w[i] << "  ";
        cout << endl;
Francois Gygi committed
148 149 150 151 152 153
        if ( info != 0 )
        {
          cerr << " AndersonMixer: Error in dsyev" << endl;
          cerr << " info = " << info << endl;
          exit(1);
        }
154

Francois Gygi committed
155 156 157 158 159 160
        // solve for theta
        // theta_i = sum_j
        for ( int k = 0; k < n_; k++ )
        {
          // correct only if eigenvalue w[k] is large enough compared to the
          // largest eigenvalue
161 162
          const double eig_ratio = 1.e-14;
          if ( w[k] > eig_ratio * w[n_-1] )
Francois Gygi committed
163 164 165 166 167 168 169
          {
            const double fac = 1.0 / w[k];
            for ( int i = 0; i < n_; i++ )
              for ( int j = 0; j < n_; j++ )
                theta[i] += fac * a[i+k*n_] * a[j+k*n_] * b[j];
          }
        }
170
      }
Francois Gygi committed
171
      else
172
      {
Francois Gygi committed
173
        // solve the linear system directly
174 175

        // Tikhonov regularization parameter
176 177 178 179 180 181 182 183
        // adjust the parameter until the norm of theta is < 1.0
        double tikhonov_parameter = 1.e-12;
        bool norm_ok = false;
        valarray<double> asave(a);
        valarray<double> bsave(b);
        int iter = 0;
        const int maxiter = 100;
        while ( !norm_ok && iter < maxiter )
184
        {
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231
          a = asave;
          b = bsave;
          for ( int i = 0; i < n_; i++ )
            a[i+i*n_] += tikhonov_parameter;

          char uplo = 'L';
          int nrhs = 1;
          valarray<int> ipiv(n_);
          valarray<double> work(n_);
          int info;
          dsysv(&uplo,&n_,&nrhs,&a[0],&n_,&ipiv[0],
                &b[0],&n_,&work[0],&n_,&info);
          if ( info != 0 )
          {
            cerr << " AndersonMixer: Error in dsysv" << endl;
            cerr << " info = " << info << endl;
            exit(1);
          }
          // the vector b now contains the solution
          theta = b;

          // check condition on the norm of theta
          norm_ok = true;
#if 0
          // unit simplex criterion
          double theta_sum = 0.0;
          for ( int i = 0; i < theta.size(); i++ )
          {
            theta_sum += theta[i];
            norm_ok &= theta[i] >= 0.0;
          }
          norm_ok &= fabs(theta_sum) <= 1.0;
#endif
#if 0
          // infinity norm criterion
          for ( int i = 0; i < theta.size(); i++ )
            norm_ok &= fabs(theta[i]) <  3.0;
#endif
#if 1
          // 2-norm criterion
          double theta_sum = 0.0;
          for ( int i = 0; i < theta.size(); i++ )
          {
            theta_sum += theta[i] * theta[i];
          }
          norm_ok = theta_sum <= 1.0;
#endif
Francois Gygi committed
232 233 234 235 236
          // cout << " tp = " << tikhonov_parameter
          //      << " AndersonMixer: theta = ";
          // for ( int i = 0; i < theta.size(); i++ )
          //   cout << theta[i] << " ";
          // cout << endl;
237 238 239

          tikhonov_parameter *= 2.0;
          iter++;
240
        }
Francois Gygi committed
241
      }
Francois Gygi committed
242 243 244 245 246

      cout << " AndersonMixer: theta = ";
      for ( int i = 0; i < theta.size(); i++ )
        cout << theta[i] << " ";
      cout << endl;
247

Francois Gygi committed
248
    }
Francois Gygi committed
249

250
#if USE_MPI
Francois Gygi committed
251
    // broadcast theta from task 0
252
    if ( pcomm_ != 0 )
Francois Gygi committed
253
    {
254
      MPI_Bcast(&theta[0],n_,MPI_DOUBLE,0,*pcomm_);
Francois Gygi committed
255 256 257
    }
#endif

258
  } // if n_ > 0
259

260 261 262 263 264 265 266 267 268 269 270 271 272
  // fbar = f[k] + sum_i theta_i * ( f[k] - f[kmi] )
  // xbar = x[k] + sum_i theta_i * ( x[k] - x[kmi] )
  for ( int l = 0; l < m_; l++ )
  {
    fbar[l] = f_[k_][l];
    xbar[l] = x_[k_][l];
  }
  for ( int i = 0; i < n_; i++ )
  {
    const int kmi = ( k_ - i + nmax_ ) % ( nmax_ + 1 );
    assert(kmi>=0);
    assert(kmi<nmax_+1);
    for ( int l = 0; l < m_; l++ )
Francois Gygi committed
273
    {
274 275
      fbar[l] -= theta[i] * ( f_[k_][l] - f_[kmi][l] );
      xbar[l] -= theta[i] * ( x_[k_][l] - x_[kmi][l] );
276
    }
277
  }
Francois Gygi committed
278
}