20 #ifndef JDFTX_CORE_MINIMIZE_H 21 #define JDFTX_CORE_MINIMIZE_H 26 #include <core/MinimizeParams.h> 48 virtual void step(
const Vector& dir,
double alpha)=0;
51 virtual double compute(Vector* grad)=0;
59 virtual bool report(
int iter) {
return false; }
66 virtual double sync(
double x)
const {
return x; }
77 Linmin getLinmin(
const MinimizeParams& params)
const;
78 double lBFGS(
const MinimizeParams& params);
91 virtual Vector hessian(
const Vector&)
const=0;
97 virtual double sync(
double x)
const {
return x; }
108 #include <core/Minimize_linmin.h> 109 #include <core/Minimize_lBFGS.h> 120 bool forceGradDirection =
true;
124 double alpha = alphaT;
126 double gKNorm = 0.0, gKNormPrev = 0.0;
129 Linmin linmin = getLinmin(p);
137 fprintf(p.
fpLog,
"%s\tState modified externally: resetting search direction.\n", p.
linePrefix);
139 forceGradDirection =
true;
146 fprintf(p.
fpLog,
" |grad|_K: %10.3le alpha: %10.3le",
sqrt(gKNorm/p.
nDim), alpha);
150 if(!forceGradDirection)
151 {
double dotgd =
sync(
dot(g,d));
152 double dotgPrevKg =
sync(
dot(gPrev, Kg));
155 double cgtest = dotgPrevKg/
sqrt(gKNorm*gKNormPrev);
156 fprintf(p.
fpLog,
" linmin: %10.3le", linmin);
157 fprintf(p.
fpLog,
" cgtest: %10.3le", cgtest);
160 switch(currentDirUpdateScheme)
168 { fprintf(p.
fpLog,
"\n%sEncountered beta<0, resetting CG.", p.
linePrefix);
172 forceGradDirection =
false;
176 fflush(p.
fpLog);
return E;
178 if(ediffCheck.checkConvergence(E))
179 { fprintf(p.
fpLog,
"%sConverged (|Delta %s|<%le for %d iters).\n",
181 fflush(p.
fpLog);
return E;
183 if(!std::isfinite(gKNorm))
184 { fprintf(p.
fpLog,
"%s|grad|_K=%le. Stopping ...\n", p.
linePrefix, gKNorm);
185 fflush(p.
fpLog);
return E;
187 if(!std::isfinite(E))
189 fflush(p.
fpLog);
return E;
196 d *= beta;
axpy(-1.0, Kg, d);
199 if(linmin(*
this, p, d, alphaT, alpha, E, g))
214 fprintf(p.
fpLog,
"%s\tStep failed: resetting search direction.\n", p.
linePrefix);
216 forceGradDirection =
true;
220 fprintf(p.
fpLog,
"%s\tStep failed along negative gradient direction.\n", p.
linePrefix);
221 fprintf(p.
fpLog,
"%sProbably at roundoff error limit. (Stopping)\n", p.
linePrefix);
227 fprintf(p.
fpLog,
"%sNone of the convergence criteria satisfied after %d iterations.\n", p.
linePrefix, iter);
234 const double deltaMin = 1e-9;
235 const double deltaMax = 1e+1;
236 const double deltaScale = 1e+1;
238 const char* fdPrefix = fdPrefixString.c_str();
239 fprintf(p.
fpLog,
"%s--------------------------------------\n", fdPrefix);
252 double dE_ddelta =
sync(
dot(dx, g));
255 for(
double delta=deltaMin; delta<=deltaMax; delta*=deltaScale)
256 {
double dE = dE_ddelta*delta;
257 step(dx, delta-deltaPrev); deltaPrev=delta;
259 fprintf(p.
fpLog,
"%s delta=%le:\n", fdPrefix, delta);
260 fprintf(p.
fpLog,
"%s d%s Ratio: %19.16lf\n", fdPrefix, p.
energyLabel, deltaE/dE);
263 fprintf(p.
fpLog,
"%s--------------------------------------\n", fdPrefix);
264 step(dx, -deltaPrev);
270 Vector r =
clone(rhs);
axpy(-1.0, hessian(state), r);
272 double beta=0.0, rdotzPrev=0.0, rdotz =
sync(
dot(r, z));
275 double rzNorm =
sqrt(fabs(rdotz)/p.
nDim);
284 { beta = rdotz/rdotzPrev;
285 d *= beta;
axpy(1.0, z, d);
289 Vector w = hessian(d);
290 double alpha = rdotz/
sync(
dot(w,d));
291 axpy(alpha, d, state);
297 double rzNorm =
sqrt(fabs(rdotz)/p.
nDim);
298 fprintf(p.
fpLog,
"%sIter: %3d sqrt(|r.z|): %12.6le alpha: %12.6le beta: %13.6le\n",
303 fprintf(p.
fpLog,
"%sGradient did not converge within threshold in %d iterations\n", p.
linePrefix, iter); fflush(p.
fpLog);
308 #endif // JDFTX_CORE_MINIMIZE_H double alphaTmin
minimum value of the test-step size (algorithm gives up when difficulties cause alphaT to fall below ...
Definition: MinimizeParams.h:60
virtual void constrain(Vector &)
Definition: Minimize.h:63
Limited memory version of the BFGS algorithm.
Definition: MinimizeParams.h:36
Definition: Minimize_linmin.h:27
ScalarField sqrt(const ScalarField &)
Elementwise square root (preserve input)
Definition: Minimize.h:86
void randomize(TptrCollection &x)
Initialize to normal random numbers:
Definition: ScalarFieldArray.h:157
int nIterations
Maximum number of iterations (default 100)
Definition: MinimizeParams.h:48
complex dot(const Tptr &X, const Tptr &Y)
Definition: Operators.h:196
Hestenes-Stiefel (preconditioned) conjugate gradients.
Definition: MinimizeParams.h:35
virtual bool report(int iter)
Definition: Minimize.h:59
Parameters to control the minimization algorithm.
Definition: MinimizeParams.h:29
double energyDiffThreshold
stop when energy change is below this for nEnergyDiff successive iterations (default: 0) ...
Definition: MinimizeParams.h:56
virtual double sync(double x) const
Override to synchronize scalars over MPI processes (if the same minimization is happening in sync ove...
Definition: Minimize.h:66
double alphaTstart
initial value for the test-step size (default: 1.0)
Definition: MinimizeParams.h:59
int nEnergyDiff
number of successive iterations for energyDiffThreshold check (default: 2)
Definition: MinimizeParams.h:57
Steepest Descent (always along negative (preconditioned) gradient)
Definition: MinimizeParams.h:37
virtual Vector precondition(const Vector &grad)
Definition: Minimize.h:55
double knormThreshold
stop when norm of residual against preconditioner falls below this (default: 0)
Definition: MinimizeParams.h:55
void fdTest(const MinimizeParams ¶ms)
Tptr clone(const Tptr &X)
Clone (NOTE: operator= is by reference for the ScalarField classes)
Definition: Operators.h:111
virtual Vector precondition(const Vector &v) const
Override to enable preconditioning: return the preconditioned vector, given a vector.
Definition: Minimize.h:94
bool updateTestStepSize
set alphaT=alpha after every iteration if true (default: true)
Definition: MinimizeParams.h:61
bool killFlag
Flag set by signal handlers - all compute loops should quit cleanly when this is set.
virtual double sync(double x) const
Override to synchronize scalars over MPI processes (if the same minimization is happening in sync ove...
Definition: Minimize.h:97
double minimize(const MinimizeParams ¶ms)
Minimize this objective function with algorithm controlled by params and return the minimized value...
Polak-Ribiere (preconditioned) conjugate gradients (default)
Definition: MinimizeParams.h:33
const char * linePrefix
prefix for each output line of minimizer, useful for nested minimizations (default "CG\t") ...
Definition: MinimizeParams.h:52
const char * energyLabel
Label for the minimized quantity (default "E")
Definition: MinimizeParams.h:53
Fletcher-Reeves (preconditioned) conjugate gradients.
Definition: MinimizeParams.h:34
int nDim
Dimension of optimization space; used only for knormThreshold (default 1)
Definition: MinimizeParams.h:49
int solve(const Vector &rhs, const MinimizeParams ¶ms)
const char * energyFormat
printf format for the minimized quantity (default "%22.15le")
Definition: MinimizeParams.h:54
FILE * fpLog
Stream to log iterations to.
Definition: MinimizeParams.h:51
Definition: Minimize.h:45
DirectionUpdateScheme
Search direction update scheme.
Definition: MinimizeParams.h:32
virtual void step(const Vector &dir, double alpha)=0
Move the state in parameter space along direction dir with scale alpha.
void axpy(double alpha, const Tptr &X, Tptr &Y)
Generic axpy for complex data types (Note: null pointers are treated as zero)
Definition: Operators.h:158
std::basic_string< char, ichar_traits > string
Case-insensitive string.
Definition: string.h:42
virtual double compute(Vector *grad)=0
Returns the objective function at the current state and store the gradient in grad, if non-null.
bool fdTest
whether to perform a finite difference test before each minimization (default false) ...
Definition: MinimizeParams.h:70
Vector state
the location of the minimum, obtained by solving hessian * state == rhs
Definition: Minimize.h:88