Rewrite of CellStepper class hierarchy, CGOptimizer, LineMinimizer

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

src/CGCellStepper.C

@@ -17,234 +17,155 @@
 ////////////////////////////////////////////////////////////////////////////////
 
 #include "CGCellStepper.h"
-#include <valarray>
+#include "CGOptimizer.h"
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
-CGCellStepper::CGCellStepper(Sample& s) : CellStepper(s), first_step_(true),
-  sigma1_(0.1), sigma2_(0.5)
+CGCellStepper::CGCellStepper(Sample& s) : CellStepper(s), 
+  cgopt_(CGOptimizer(3*s.atoms.size()+9)), cell0(s_.atoms.cell())
 {
   nat_ = atoms_.size();
+  cgopt_.set_alpha_start(0.5);
+  cgopt_.set_beta_max(5.0);
+#ifdef DEBUG
+  cgopt_.set_debug_print();
+#endif
 
-  xc_.resize(3*nat_+9);
-  pc_.resize(3*nat_+9);
-  fc_.resize(3*nat_+9);
+  rp_.resize(s.atoms.nsp());
+  for ( int is = 0; is < rp_.size(); is++ )
+    rp_[is].resize(3*atoms_.na(is));
 
-  rc_.resize(atoms_.nsp());
-  rp_.resize(atoms_.nsp());
-  for ( int is = 0; is < rc_.size(); is++ )
+  // store full strain tensor u for consistency of dot products
+  u_.resize(9);
+  up_.resize(9);
+  for ( int i = 0; i < u_.size(); i++ )
   {
-    rc_[is].resize(3*atoms_.na(is));
-    rp_[is].resize(3*atoms_.na(is));
+    u_[i] = 0.0;
+    up_[i] = 0.0;
   }
-
-  linmin_.set_sigma1(sigma1_);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-void CGCellStepper::compute_new_cell(double e0, const valarray<double>& sigma,
+void CGCellStepper::compute_new_cell(double e, const valarray<double>& sigma,
   const std::vector<std::vector< double> >& fion)
 {
   // compute new cell and ionic positions using the stress tensor sigma
   // and the forces on ions fion
-  const UnitCell& cell = s_.wf.cell();
+  const UnitCell cell = s_.wf.cell();
 
   // total number of dofs: 3* natoms + cell parameters
-  valarray<double> f0(3*nat_+9), x0(3*nat_+9), xp(3*nat_+9);
+  valarray<double> x(3*nat_+9), xp(3*nat_+9), g(3*nat_+9);
 
-  // copy current positions into x0
+  // copy current positions into x
   vector<vector<double> > r0;
   atoms_.get_positions(r0);
   double tmp3[3];
   // convert position from r0 to tau (relative) coordinates: tau = A^-1 R
-  for ( int i = 0, is = 0; is < r0.size(); is++ )
+  for ( int is = 0, i = 0; is < r0.size(); is++ )
   {
     for ( int ia = 0; ia < r0[is].size()/3; ia++ )
     {
       cell.vecmult3x3(cell.amat_inv(),&r0[is][3*ia],tmp3);
-      x0[i++]=tmp3[0];
-      x0[i++]=tmp3[1];
-      x0[i++]=tmp3[2];
+      x[i++]=tmp3[0];
+      x[i++]=tmp3[1];
+      x[i++]=tmp3[2];
     }
   }
 
-  // convert forces on positions to forces on tau coordinates, store in f0
+  // copy current strain tensor into x
+  for ( int i = 0; i < 9; i++ )
+    x[3*nat_+i] = u_[i];
+
+  // convert forces on positions to forces  on tau coordinates, store -f in g
   // f = A^-1 * fion
   for ( int i = 0, is = 0; is < fion.size(); is++ )
   {
     for ( int ia = 0; ia < fion[is].size()/3; ia++ )
     {
       cell.vecmult3x3(cell.amat_inv(),&fion[is][3*ia],tmp3);
-      f0[i++]=tmp3[0];
-      f0[i++]=tmp3[1];
-      f0[i++]=tmp3[2];
+      g[i++]=-tmp3[0];
+      g[i++]=-tmp3[1];
+      g[i++]=-tmp3[2];
     }
   }
 
-  // convert descent direction dcell to cell space from the stress tensor
-  // dcell = sigma * A
-  valarray<double> dcell(9); // descent direction in cell space
-  assert(sigma.size()==6);
-  // next line: local copy of sigma to circumvent compiler error
-  valarray<double> sigma_loc(sigma);
-  cell.smatmult3x3(&sigma_loc[0],cell.amat(),&dcell[0]);
-  for ( int i = 0; i < dcell.size(); i++ )
-    f0[3*nat_+i] = dcell[i];
-
-  // cout << "dcell = " << endl;
-  // cout << dcell[0] << " " << dcell[1] << " " << dcell[2] << endl;
-  // cout << dcell[3] << " " << dcell[4] << " " << dcell[5] << endl;
-  // cout << dcell[6] << " " << dcell[7] << " " << dcell[8] << endl;
-
-  // the vector f0 now contains the derivatives of the energy in tau+cell space
-#ifdef DEBUG
-  if ( s_.ctxt_.onpe0() )
-  {
-    cout << " f0:" << endl;
-    for ( int i = 0; i < f0.size(); i++ )
-      cout << f0[i] << endl;
-  }
-#endif
+  // compute descent direction in strain space
+  // gradient g = - sigma * volume
+  g[3*nat_+0] = -sigma[0] * cell.volume();
+  g[3*nat_+1] = -sigma[3] * cell.volume();
+  g[3*nat_+2] = -sigma[5] * cell.volume();
 
-  double fp0; // derivative f'(x) at x=xp
-  bool wolfe1, wolfe2; // First and second Wolfe conditions in line search
+  g[3*nat_+3] = -sigma[3] * cell.volume();
+  g[3*nat_+4] = -sigma[1] * cell.volume();
+  g[3*nat_+5] = -sigma[4] * cell.volume();
 
-  // CG algorithm
+  g[3*nat_+6] = -sigma[5] * cell.volume();
+  g[3*nat_+7] = -sigma[4] * cell.volume();
+  g[3*nat_+8] = -sigma[2] * cell.volume();
 
-  if ( !first_step_ )
-  {
-    wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_;
-    // compute fp0: projection of forces on direction pc_
-    fp0 = 0.0;
-    for ( int i = 0; i < f0.size(); i++ )
-      fp0 -= f0[i] * pc_[i];
-
-    wolfe2 = fabs(fp0) < sigma2_ * fabs(fpc_);
-    if ( s_.ctxt_.onpe0() )
-    {
-      cout << "  CGCellStepper: fpc = " << fpc_ << endl;
-      cout << "  CGCellStepper: fp0 = " << fp0 << endl;
-      cout << "  CGCellStepper: ec = " << ec_ << " e0 = " << e0 <<  endl;
-      cout << "  CGCellStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
-           << ec_ + fpc_ * sigma1_ * alpha_ << endl;
-      cout << "  CGCellStepper: wolfe1/wolfe2 = "
-           << wolfe1 << "/" << wolfe2 << endl;
-    }
-  }
+  // the vector g now contains the gradient of the energy in tau+strain space
 
-  if ( first_step_ || (wolfe1 && wolfe2) )
-  {
-    // define new descent direction pc_
-    if ( first_step_ )
-    {
-      pc_ = f0;
-    }
-    else
-    {
-      // Polak-Ribiere definition
-      double num = 0.0, den = 0.0;
-      for ( int i = 0; i < f0.size(); i++ )
-      {
-        const double fctmp = fc_[i];
-        const double f0tmp = f0[i];
-        num += f0tmp * ( f0tmp - fctmp );
-        den += fctmp * fctmp;
-      }
-      double beta = den > 0.0 ? num/den : 0.0;
-      beta = max(beta,0.0);
-      if ( s_.ctxt_.onpe0() )
-        cout << "  CGCellStepper: beta = " << beta << endl;
-      for ( int i = 0; i < f0.size(); i++ )
-        pc_[i] = beta * pc_[i] + f0[i];
-    }
+  // project the gradient in a direction compatible with constraints
+ 
+  enforce_constraints(&g[3*nat_]);
 
-    fc_ = f0;
-    // fpc = d_e / d_alpha in direction pc
-    fpc_ = 0.0;
-    for ( int i = 0; i < f0.size(); i++ )
-      fpc_ -= fc_[i] * pc_[i];
-    ec_ = e0;
-    xc_ = x0;
-    fp0 = fpc_;
-    // reset line minimizer
-    linmin_.reset();
-  }
+  // CG algorithm
 
-#ifdef DEBUG
-  if ( s_.ctxt_.onpe0() )
-  {
-    cout << " pc:" << endl;
-    for ( int i = 0; i < pc_.size(); i++ )
-      cout << pc_[i] << endl;
-  }
+  cgopt_.compute_xp(x,e,g,xp);
 
-  // find largest component of pc_
-  double pcmax = 0.0;
-  for ( int i = 0; i < pc_.size(); i++ )
-    pcmax = max(fabs(pc_[i]),pcmax);
   if ( s_.ctxt_.onpe0() )
   {
-    cout << "CGCellStepper: largest component of pc_: "
-         << pcmax << endl;
+    cout << "  CGCellStepper: alpha = " << cgopt_.alpha() << endl;
   }
-#endif
-
-  alpha_ = linmin_.newalpha(alpha_,e0,fp0);
-
-  if ( s_.ctxt_.onpe0() )
-    cout << "  CGCellStepper: alpha = " << alpha_ << endl;
-
-  // update
-  xp = xc_ + alpha_ * pc_;
 
-  // save the change of cell parameters into dcell
   for ( int i = 0; i < 9; i++ )
-    dcell[i] = xp[3*nat_+i]-xc_[3*nat_+i];
-
-  // cellp = cell + dcell
-  D3vector a0p = cell.a(0) + D3vector(dcell[0],dcell[1],dcell[2]);
-  D3vector a1p = cell.a(1) + D3vector(dcell[3],dcell[4],dcell[5]);
-  D3vector a2p = cell.a(2) + D3vector(dcell[6],dcell[7],dcell[8]);
-  // cout << "dcell = " << endl;
-  // cout << dcell[0] << " " << dcell[1] << " " << dcell[2] << endl;
-  // cout << dcell[3] << " " << dcell[4] << " " << dcell[5] << endl;
-  // cout << dcell[6] << " " << dcell[7] << " " << dcell[8] << endl;
-  cellp = UnitCell(a0p,a1p,a2p);
-
-  // check for cell_lock constraints and modify cellp if needed
-  enforce_constraints(cell, cellp);
-
-  // compute atomic positions rc_[is][3*ia+j] from xc_
+    up_[i] = xp[3*nat_+i]; 
+
+  // enforce cell_lock constraints
+  enforce_constraints(&up_[0]);
+
+  // compute cellp = ( I + Up ) * A0
+  // symmetric matrix I+U stored in iumat: xx, yy, zz, xy, yz, xz
+  double iupmat[6];
+  iupmat[0] = 1.0 + up_[0];
+  iupmat[1] = 1.0 + up_[4];
+  iupmat[2] = 1.0 + up_[8];
+  iupmat[3] = 0.5 * ( up_[1] + up_[3] );
+  iupmat[4] = 0.5 * ( up_[5] + up_[7] );
+  iupmat[5] = 0.5 * ( up_[2] + up_[6] );
+
+  const double *a0mat = cell0.amat();
+  double apmat[9];
+  cell0.smatmult3x3(&iupmat[0],a0mat,apmat);
+
+  D3vector a0p(apmat[0],apmat[1],apmat[2]);
+  D3vector a1p(apmat[3],apmat[4],apmat[5]);
+  D3vector a2p(apmat[6],apmat[7],apmat[8]);
+  cellp.set(a0p,a1p,a2p);
+  
   // compute new atomic positions rp_[is][3*ia+j] from xp and cellp
   for ( int is = 0, i = 0; is < fion.size(); is++ )
   {
     for ( int ia = 0; ia < fion[is].size()/3; ia++ )
     {
-      tmp3[0] = xc_[i+0];
-      tmp3[1] = xc_[i+1];
-      tmp3[2] = xc_[i+2];
-      cell.vecmult3x3(cell.amat(),tmp3,&rc_[is][3*ia]);
       tmp3[0] = xp[i+0];
       tmp3[1] = xp[i+1];
       tmp3[2] = xp[i+2];
-      cell.vecmult3x3(cellp.amat(),tmp3,&rp_[is][3*ia]);
+      cellp.vecmult3x3(cellp.amat(),tmp3,&rp_[is][3*ia]);
       i+=3;
     }
   }
 
   // enforce constraints on atomic positions
-  s_.constraints.enforce_r(rc_,rp_);
-
-  first_step_ = false;
+  s_.constraints.enforce_r(r0,rp_);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 void CGCellStepper::update_cell(void)
 {
-  const UnitCell& cell = s_.wf.cell();
-
   s_.atoms.set_positions(rp_);
   s_.atoms.set_cell(cellp);
+  u_ = up_;
 
   // resize wavefunction and basis sets
   s_.wf.resize(cellp,s_.wf.refcell(),s_.wf.ecut());

src/CGCellStepper.h

@@ -20,19 +20,17 @@
 #define CGCELLSTEPPER_H
 
 #include "CellStepper.h"
-#include "LineMinimizer.h"
+#include "CGOptimizer.h"
 
 class CGCellStepper : public CellStepper
 {
   private:
 
-  bool first_step_;
+  CGOptimizer cgopt_;
+  UnitCell cell0; 
+  std::vector<std::vector<double> > rp_;
+
   int nat_;
-  std::valarray<double> xc_, pc_, fc_;
-  std::vector<std::vector<double> > rc_, rp_;
-  double ec_, fpc_;
-  double alpha_, sigma1_, sigma2_;
-  LineMinimizer linmin_;
 
   public:
 

src/CGIonicStepper.C

@@ -15,148 +15,60 @@
 // CGIonicStepper.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: CGIonicStepper.C,v 1.4 2008-09-08 15:56:18 fgygi Exp $
 
 #include "CGIonicStepper.h"
-#include <cassert>
+#include "CGOptimizer.h"
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
-void CGIonicStepper::compute_r(double e0, const vector<vector< double> >& f0)
+CGIonicStepper::CGIonicStepper(Sample& s) : IonicStepper(s), 
+  cgopt_(CGOptimizer(3*natoms_))
 {
-  double fp0;
-  bool wolfe1, wolfe2;
-
-  // check if the largest component of the force f0 is larger than max_force
-  const double max_force = 0.1;
-  double largest_force = 0.0;
-  for ( int is = 0; is < r0_.size(); is++ )
-  {
-    for ( int i = 0; i < r0_[is].size(); i++ )
-    {
-      largest_force = max(largest_force,fabs(f0[is][i]));
-    }
-  }
+  cgopt_.set_beta_max(5.0);
+}
+////////////////////////////////////////////////////////////////////////////////
+void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
+{
+  // CG algorithm
 
-  if ( largest_force > max_force )
-  {
-    if ( s_.ctxt_.onpe0() )
-      cout << "  CGIonicStepper: force exceeds limit, taking SD step " << endl;
-    // take a steepest descent step with limited displacement and exit
-    const double alpha_sd = max_force/largest_force;
-    // SD step
-    for ( int is = 0; is < r0_.size(); is++ )
+  valarray<double> x(3*natoms_),xp(3*natoms_),g(3*natoms_);
+  for ( int is = 0, i = 0; is < r0_.size(); is++ )
+    for ( int j = 0; j < r0_[is].size(); j++ )
     {
-      for ( int i = 0; i < r0_[is].size(); i++ )
-      {
-        rp_[is][i] = r0_[is][i] + alpha_sd * f0[is][i];
-      }
+      x[i] = r0_[is][j];
+      g[i] = -f0[is][j];
+      i++;
     }
-    constraints_.enforce_r(r0_,rp_);
-    rm_ = r0_;
-    r0_ = rp_;
-    atoms_.set_positions(r0_);
-    // reset the CG algorithm
-    first_step_ = true;
-    return;
-  }
 
-  // CG algorithm
+  cgopt_.compute_xp(x,e0,g,xp);
 
-  // define the descent direction
-  if ( first_step_ )
+  // check largest displacement
+  // max_disp: largest acceptable displacement
+  const double max_disp = 0.05;
+  double largest_disp = 0.0;
+  for ( int i = 0; i < xp.size(); i++ )
+    largest_disp = max(largest_disp,fabs(xp[i]-x[i]));
+  if ( largest_disp > max_disp )
   {
-    // first step: set descent direction pc_ = f0
-    pc_ = f0;
-    // fp0: dE/dalpha in the direction -pc_
-    fp0 = 0.0;
-    for ( int is = 0; is < r0_.size(); is++ )
-      for ( int i = 0; i < r0_[is].size(); i++ )
-        fp0 -= f0[is][i] * pc_[is][i];
-    // initialize fc_, fpc_, ec_, rc_
-    fc_ = f0;
-    // fpc = d_e / d_alpha in direction pc
-    fpc_ = fp0;
-    ec_ = e0;
-    rc_ = r0_;
-  }
-  else
-  {
-    // check if the Wolfe conditions are satisfied
-    wolfe1 = e0 < ec_ + fpc_ * sigma1_ * alpha_;
-    // fp0 = -proj(f0,pc)
-    fp0 = 0.0;
-    for ( int is = 0; is < r0_.size(); is++ )
-      for ( int i = 0; i < r0_[is].size(); i++ )
-        fp0 -= f0[is][i] * pc_[is][i];
-    wolfe2 = fabs(fp0) < sigma2_ * fabs(fpc_);
     if ( s_.ctxt_.onpe0() )
-    {
-      cout << "  CGIonicStepper: fpc = " << fpc_ << endl;
-      cout << "  CGIonicStepper: fp0 = " << fp0 << endl;
-      cout << "  CGIonicStepper: ec = " << ec_ << " e0 = " << e0 <<  endl;
-      cout << "  CGIonicStepper: ec_ + fpc_ * sigma1_ * alpha_ ="
-           << ec_ + fpc_ * sigma1_ * alpha_ << endl;
-      cout << "  CGIonicStepper: wolfe1/wolfe2 = "
-           << wolfe1 << "/" << wolfe2 << endl;
-    }
-
-    if ( wolfe1 && wolfe2 )
-    {
-      // done with this descent direction
-      // define new descent direction pc_ using the Fletcher-Reeves formula
-      double num = 0.0, den = 0.0;
-      for ( int is = 0; is < r0_.size(); is++ )
-      {
-        for ( int i = 0; i < r0_[is].size(); i++ )
-        {
-          const double fctmp = fc_[is][i];
-          const double f0tmp = f0[is][i];
-          // Fletcher-Reeves definition
-          num += f0tmp * f0tmp;
-          den += fctmp * fctmp;
-        }
-      }
-      double beta = den > 0.0 ? num/den : 0.0;
-      if ( s_.ctxt_.onpe0() )
-        cout << "  CGIonicStepper: beta = " << beta << endl;
-      for ( int is = 0; is < r0_.size(); is++ )
-      {
-        for ( int i = 0; i < r0_[is].size(); i++ )
-          pc_[is][i] = beta * pc_[is][i] + f0[is][i];
-      }
-      // initialize fc_, fpc_, ec_, rc_
-      // fp0 = d_e / d_alpha in direction pc_
-      fp0 = 0.0;
-      for ( int is = 0; is < r0_.size(); is++ )
-        for ( int i = 0; i < r0_[is].size(); i++ )
-          fp0 -= f0[is][i] * pc_[is][i];
-      fc_ = f0;
-      fpc_ = fp0;
-      ec_ = e0;
-      rc_ = r0_;
-
-      // reset the line minimizer
-      linmin_.reset();
-    }
+      cout << "  CGIonicStepper: displacement exceeds limit, rescaling" << endl;
+    // rescale displacement and reset the CG optimizer
+    xp = x + (max_disp/largest_disp) * (xp - x);
+    
+    cgopt_.reset();
   }
 
-  // the descent direction pc_ is defined
-
-  alpha_ = linmin_.newalpha(alpha_,e0,fp0);
-
   if ( s_.ctxt_.onpe0() )
-    cout << "  CGIonicStepper: alpha = " << alpha_ << endl;
+  {
+    cout << "  CGIonicStepper: alpha = " << cgopt_.alpha() << endl;
+  }
 
-  // rp = rc + alpha * pc
-  for ( int is = 0; is < r0_.size(); is++ )
-    for ( int i = 0; i < r0_[is].size(); i++ )
-      rp_[is][i] = rc_[is][i] + alpha_ * pc_[is][i];
+  for ( int is = 0, i = 0; is < r0_.size(); is++ )
+    for ( int j = 0; j < r0_[is].size(); j++ )
+      rp_[is][j] = xp[i++];
 
   constraints_.enforce_r(r0_,rp_);
   rm_ = r0_;
   r0_ = rp_;
   atoms_.set_positions(r0_);
-
-  first_step_ = false;
 }

src/CGIonicStepper.h

@@ -15,34 +15,23 @@
 // CGIonicStepper.h:
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: CGIonicStepper.h,v 1.4 2008-09-08 15:56:18 fgygi Exp $
 
 #ifndef CGIONICSTEPPER_H
 #define CGIONICSTEPPER_H
 
 #include "IonicStepper.h"
-#include "LineMinimizer.h"
-#include <vector>
+#include "CGOptimizer.h"
 
 class CGIonicStepper : public IonicStepper
 {
   private:
 
-  bool first_step_;
-  std::vector<std::vector< double> > rc_;
-  std::vector<std::vector< double> > pc_;
-  std::vector<std::vector< double> > fc_;
-  double ec_, fpc_;
-  double alpha_, sigma1_, sigma2_;
-  LineMinimizer linmin_;
+  CGOptimizer cgopt_;
 
   public:
 
-  CGIonicStepper(Sample& s) : IonicStepper(s), first_step_(true),
-    sigma1_(0.1), sigma2_(0.5) { linmin_.set_sigma1(sigma1_); }
-
+  CGIonicStepper(Sample& s);
   void compute_r(double e0, const std::vector<std::vector< double> >& f0);
   void compute_v(double e0, const std::vector<std::vector< double> >& f0) {}
 };
-
 #endif

src/CGOptimizer.C

+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2011 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// CGOptimizer.C
+//
+////////////////////////////////////////////////////////////////////////////////
+#include "CGOptimizer.h"
+#include <iostream>
+#include <cassert>
+using namespace std;
+////////////////////////////////////////////////////////////////////////////////
+double CGOptimizer::norm2(const std::valarray<double>& v)
+{
+  return inner_product(&v[0],&v[v.size()],&v[0],0.0);
+}
+////////////////////////////////////////////////////////////////////////////////
+void CGOptimizer::compute_xp(const valarray<double>& x, const double f, 
+                             const valarray<double>& g, valarray<double>& xp)
+{
+  // Use the function value f and the gradient g at x to generate a new point xp
+  // using the Fletcher-Powell CG algorithm
+  // return xp=x if the 2-norm of g is smaller than tol
+  const double tol = 1.0e-18;
+
+  assert(x.size()==n_ && g.size()==n_ && xp.size()==n_);
+
+  double fp;
+  // define the descent direction
+  if ( first_step_ )
+  {
+    p_ = -g;
+
+    x0_ = x;
+    f0_ = f;
+
+    g0norm2_ = norm2(g);
+    if ( g0norm2_ < tol )
+    {
+      xp = x;
+      return;
+    }
+    fp = -g0norm2_;
+    fp0_ = fp;
+    linmin_.reset();
+    alpha_ = linmin_.next_alpha(alpha_,f,fp);
+    if ( debug_print )
+      cout << "  CGOptimizer: first_step: alpha=" << alpha_ 
+           << " f=" << f << " fp=" << fp << endl;
+    
+    xp = x0_ + alpha_ * p_;
+    first_step_ = false;
+  }
+  else
+  {
+    // fp: derivative along the current descent direction p_
+    // fp = df(x0+alpha*p)/dalpha at x
+    fp = inner_product(&g[0],&g[g.size()],&p_[0],0.0);
+    double new_alpha = linmin_.next_alpha(alpha_,f,fp);
+    if ( debug_print )
+      cout << "  CGOptimizer: alpha=" << alpha_ 
+           << " f=" << f << " fp=" << fp << endl;
+
+    if ( linmin_.fail() )
+    {
+      // line minimization failed
+      if ( debug_print )
+        cout << "CGOptimizer: line minimization failed" << endl;
+
+      // restart from current point
+      p_ = -g;
+
+      x0_ = x;
+      f0_ = f;
+
+      g0norm2_ = norm2(g);
+      if ( g0norm2_ < tol )
+      {
+        xp = x;
+        return;
+      }
+      fp = -g0norm2_;
+      fp0_ = fp;
+      linmin_.reset();
+      alpha_ = linmin_.next_alpha(alpha_,f,fp);
+      if ( debug_print )
+        cout << "  CGOptimizer: restart after fail: alpha=" << alpha_ 
+           << " f=" << f << " fp=" << fp << endl;
+    
+      xp = x0_ + alpha_ * p_;
+      first_step_ = false;
+    }
+
+    if ( linmin_.done() )
+    {
+      // wolfe1_ && wolfe2_ are true at alpha_
+      if ( debug_print )
+        cout << "  CGOptimizer: done with current descent direction" << endl;
+      // define a new descent direction p_ using the Fletcher-Reeves formula
+      assert(g0norm2_ > 0.0);
+      double beta = norm2(g) / g0norm2_;
+      if ( beta_max_ > 0.0 && beta > beta_max_ )
+      {
+        if ( debug_print )
+          cout << "  CGOptimizer: beta exceeds beta_max " << endl;
+        beta = min(beta, beta_max_);
+      }
+      if ( debug_print )
+        cout << "  CGOptimizer: beta = " << beta << endl;
+      p_ = beta * p_ - g;
+
+      x0_ = x;
+      f0_ = f;
+      // recalculate f0, fp0
+      // fp0 = d_e / d_alpha in direction pc_
+      fp0_ = inner_product(&g[0],&g[g.size()],&p_[0],0.0);
+      g0norm2_ = norm2(g);
+      fp = fp0_;
+      assert(fp<0.0 && "CGOptimizer: p_ not a descent direction");
+
+      // reset the line minimizer
+      linmin_.reset();
+      alpha_ = linmin_.next_alpha(alpha_,f,fp);
+      if ( debug_print )
+        cout << "  CGOptimizer: restart: alpha=" << alpha_ 
+             << " f=" << f << " fp=" << fp << endl;
+    }
+    else
+    {
+      // not done with the current descent direction
+      alpha_ = new_alpha;
+    }
+    xp = x0_ + alpha_ * p_;
+  }
+}