testAndersonMixer.C 2.83 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
// testAndersonMixer.C
//
////////////////////////////////////////////////////////////////////////////////

#include <iostream>
Francois Gygi committed
20
#include <vector>
Francois Gygi committed
21 22 23
using namespace std;
#include "AndersonMixer.h"

Francois Gygi committed
24
// use: testAndersonMixer ndim nmax niter
Francois Gygi committed
25 26
int main(int argc, char** argv)
{
Francois Gygi committed
27 28
  int mype;
  MPI_Comm comm = MPI_COMM_WORLD;
Francois Gygi committed
29
  MPI_Init(&argc,&argv);
Francois Gygi committed
30 31
  MPI_Comm_rank(MPI_COMM_WORLD,&mype);
  if ( argc != 4 )
Francois Gygi committed
32
  {
Francois Gygi committed
33 34 35 36 37 38 39 40 41
    cout << " use: testAndersonMixer ndim nmax niter" << endl;
    return 1;
  }
  // ndim: dimension of vector x
  const int ndim = atoi(argv[1]);
  // nmax: maximum dimension of subspace used for acceleration
  const int nmax = atoi(argv[2]);
  // niter: number of iterations
  const int niter = atoi(argv[3]);
42

Francois Gygi committed
43 44 45 46
  char processor_name[MPI_MAX_PROCESSOR_NAME];
  int namelen;
  PMPI_Get_processor_name(processor_name,&namelen);
  // cout << " Process " << ctxt.mype() << " on " << processor_name << endl;
Francois Gygi committed
47

Francois Gygi committed
48 49 50 51 52 53
  const double alpha = 0.1;
  vector<double> x,f,xbar,fbar;
  x.resize(ndim);
  f.resize(ndim);
  xbar.resize(ndim);
  fbar.resize(ndim);
54

Francois Gygi committed
55
  AndersonMixer mixer(ndim,nmax,&comm);
56

Francois Gygi committed
57 58
  for ( int i = 0; i < ndim; i++ )
    x[i] = (i+5);
59

Francois Gygi committed
60 61 62
  for ( int iter = 0; iter < niter; iter++ )
  {
    // compute gradient
Francois Gygi committed
63
#if 1
Francois Gygi committed
64 65 66 67
    // quadratic form
    for ( int i = 0; i < ndim; i++ )
      f[i] = -(i+1) * (x[i]-i);
      //f[i] = -(0.1*i+1)*(x[i]-i);
Francois Gygi committed
68
#if 1
Francois Gygi committed
69 70 71
    // precondition f
    for ( int i = 0; i < ndim; i++ )
      f[i] *= 1.0/(i+3);
Francois Gygi committed
72 73
#endif

Francois Gygi committed
74
#else
Francois Gygi committed
75 76 77 78 79 80 81 82 83
    // Rosenbrock function
    assert(ndim%2==0);
    // f = sum_i^(n-1) (1-x_i)^2 + 100*(x_i+1 - x_i^2)^2
    f[0] = 0.0;
    for ( int i = 0; i < ndim-1; i++ )
    {
      f[i+1] = 200 * ( x[i+1] - x[i]*x[i] );
      f[i] += -2.0 * ( 1.0 - x[i] ) + 200 * x[i] * ( x[i+1] - x[i]*x[i] );
    }
Francois Gygi committed
84
#endif
Francois Gygi committed
85

Francois Gygi committed
86 87 88 89 90 91 92
    double resnorm = 0.0;
    for ( int i = 0; i < ndim; i++ )
      resnorm += f[i]*f[i];
    double rbuf;
    MPI_Allreduce(&resnorm,&rbuf,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    resnorm = rbuf;
    if ( mype == 0 )
Francois Gygi committed
93
      cout << " resnorm: " << sqrt(resnorm) << endl;
Francois Gygi committed
94
    mixer.update(&x[0],&f[0],&xbar[0],&fbar[0]);
95

Francois Gygi committed
96
#if 0
Francois Gygi committed
97 98 99 100
    // precondition fbar
    for ( int i = 0; i < ndim; i++ )
      fbar[i] /= (i+3);
      //fbar[i] *= 1.0/(i+2);
Francois Gygi committed
101
#endif
102

Francois Gygi committed
103 104
    for ( int i = 0; i < ndim; i++ )
      x[i] = xbar[i] + alpha * fbar[i];
Francois Gygi committed
105
  }
Francois Gygi committed
106

Francois Gygi committed
107 108 109
  MPI_Finalize();
  return 0;
}