20 #ifndef JDFTX_CORE_MINIMIZE_LBFGS_H 21 #define JDFTX_CORE_MINIMIZE_LBFGS_H 23 #include <core/Minimize_linmin.h> 33 double E = sync(compute(&g));
34 Vector Kg = precondition(g);
39 double linminTest = 0.;
47 std::list< std::shared_ptr<History> > history;
51 Linmin linmin = getLinmin(p);
58 { E = sync(compute(&g));
60 fprintf(p.
fpLog,
"%s\tState modified externally: resetting history.\n", p.
linePrefix);
65 double gKnorm = sync(
dot(g,Kg));
69 if(alpha) fprintf(p.
fpLog,
" alpha: %10.3le", alpha);
70 if(linminTest) fprintf(p.
fpLog,
" linmin: %10.3le", linminTest);
77 fflush(p.
fpLog);
return E;
79 if(ediffCheck.checkConvergence(E))
80 { fprintf(p.
fpLog,
"%sConverged (|Delta %s|<%le for %d iters).\n",
82 fflush(p.
fpLog);
return E;
84 if(!std::isfinite(gKnorm))
85 { fprintf(p.
fpLog,
"%s|grad|_K=%le. Stopping ...\n", p.
linePrefix, gKnorm);
86 fflush(p.
fpLog);
return E;
90 fflush(p.
fpLog);
return E;
95 auto h = std::make_shared<History>();
100 std::stack<double> a;
101 for(
auto h=history.rbegin(); h!=history.rend(); h++)
102 { a.push( (*h)->rho * sync(
dot((*h)->s, d)) );
103 axpy(-a.top(), (*h)->Ky, d);
105 if(gamma) d *= gamma;
106 for(
auto h=history.begin(); h!=history.end(); h++)
107 {
double b = (*h)->rho * sync(
dot((*h)->Ky, d));
108 axpy(a.top()-b, (*h)->s, d);
112 if((
int)history.size() == p.
history) history.pop_front();
116 if(!linmin(*
this, p, d, p.
alphaTstart, alpha, E, g))
120 E = sync(compute(&g));
121 Kg = precondition(g);
124 fprintf(p.
fpLog,
"%s\tStep failed: resetting history.\n", p.
linePrefix);
133 fprintf(p.
fpLog,
"%s\tStep failed along negative gradient direction.\n", p.
linePrefix);
134 fprintf(p.
fpLog,
"%sProbably at roundoff error limit. (Stopping)\n", p.
linePrefix);
141 linminTest = sync(
dot(g,d))/
sqrt(sync(
dot(g,g))*sync(
dot(d,d)));
143 h->Ky = Kg; Kg = precondition(g);
144 h->Ky *= -1;
axpy(1., Kg, h->Ky);
145 y *= -1;
axpy(1., g, y);
146 double ydots = sync(
dot(y, h->s));
148 gamma = ydots / sync(
dot(y, h->Ky));
149 history.push_back(h);
151 fprintf(p.
fpLog,
"%sNone of the convergence criteria satisfied after %d iterations.\n", p.
linePrefix, iter);
155 #endif //JDFTX_CORE_MINIMIZE_LBFGS_H Definition: Minimize_linmin.h:27
ScalarField sqrt(const ScalarField &)
Elementwise square root (preserve input)
int nIterations
Maximum number of iterations (default 100)
Definition: MinimizeParams.h:48
complex dot(const Tptr &X, const Tptr &Y)
Definition: Operators.h:196
int history
Number of past variables and residuals to store (BFGS only)
Definition: MinimizeParams.h:50
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
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
double knormThreshold
stop when norm of residual against preconditioner falls below this (default: 0)
Definition: MinimizeParams.h:55
Tptr clone(const Tptr &X)
Clone (NOTE: operator= is by reference for the ScalarField classes)
Definition: Operators.h:111
double clock_sec()
Elapsed time in microseconds (from start of program)
bool killFlag
Flag set by signal handlers - all compute loops should quit cleanly when this is set.
const char * linePrefix
prefix for each output line of minimizer, useful for nested minimizations (default "CG\t") ...
Definition: MinimizeParams.h:52
Nonlinear minimization templates.
const char * energyLabel
Label for the minimized quantity (default "E")
Definition: MinimizeParams.h:53
int nDim
Dimension of optimization space; used only for knormThreshold (default 1)
Definition: MinimizeParams.h:49
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
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