From b68259720a37afdf26cf693f4783945970535ef1 Mon Sep 17 00:00:00 2001 From: Erik Schultheis Date: Fri, 30 Oct 2020 14:44:55 +0200 Subject: [PATCH 1/4] reduce scope of local variables --- newton.cpp | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/newton.cpp b/newton.cpp index 3c1dd98..505f1a6 100644 --- a/newton.cpp +++ b/newton.cpp @@ -95,11 +95,8 @@ NEWTON::~NEWTON() void NEWTON::newton(double *w) { int n = fun_obj->get_nr_variable(); - int i, cg_iter; - double step_size; - double f, fold, actred; double init_step_size = 1; - int search = 1, iter = 1, inc = 1; + int search = 1, inc = 1; double *s = new double[n]; double *r = new double[n]; double *g = new double[n]; @@ -109,14 +106,14 @@ void NEWTON::newton(double *w) // calculate gradient norm at w=0 for stopping condition. double *w0 = new double[n]; - for (i=0; ifun(w0); fun_obj->grad(w0, g); double gnorm0 = dnrm2_(&n, g, &inc); delete [] w0; - f = fun_obj->fun(w); + double f = fun_obj->fun(w); info("init f %5.3e\n", f); fun_obj->grad(w, g); double gnorm = dnrm2_(&n, g, &inc); @@ -125,15 +122,15 @@ void NEWTON::newton(double *w) search = 0; double *w_new = new double[n]; - while (iter <= max_iter && search) + for (int iter=1; iter <= max_iter && search; ) { fun_obj->get_diag_preconditioner(M); - for(i=0; ilinesearch_and_update(w, s, & f, g, init_step_size); + double fold = f; + double step_size = fun_obj->linesearch_and_update(w, s, & f, g, init_step_size); if (step_size == 0) { @@ -143,7 +140,7 @@ void NEWTON::newton(double *w) info("iter %2d f %5.3e |g| %5.3e CG %3d step_size %4.2e \n", iter, f, gnorm, cg_iter, step_size); - actred = fold - f; + double actred = fold - f; iter++; fun_obj->grad(w, g); From bd027d59aafb692b7d92dddef6c64893c3920b5f Mon Sep 17 00:00:00 2001 From: Erik Schultheis Date: Fri, 30 Oct 2020 14:48:24 +0200 Subject: [PATCH 2/4] automatic memory management using std::vector --- newton.cpp | 88 ++++++++++++++++++++++++------------------------------ 1 file changed, 39 insertions(+), 49 deletions(-) diff --git a/newton.cpp b/newton.cpp index 505f1a6..b0165ab 100644 --- a/newton.cpp +++ b/newton.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include "newton.h" #ifndef min @@ -39,7 +40,7 @@ double function::linesearch_and_update(double *w, double *s, double *f, double * double eta = 0.01; int n = get_nr_variable(); int max_num_linesearch = 20; - double *w_new = new double[n]; + std::vector w_new(n); double fold = *f; for (int i=0;iget_nr_variable(); double init_step_size = 1; - int search = 1, inc = 1; - double *s = new double[n]; - double *r = new double[n]; - double *g = new double[n]; + int inc = 1; + std::vector s(n); + std::vector r(n); + std::vector g(n); const double alpha_pcg = 0.01; - double *M = new double[n]; + std::vector M(n, 0.0); // calculate gradient norm at w=0 for stopping condition. - double *w0 = new double[n]; - for (int i=0; ifun(w0); - fun_obj->grad(w0, g); - double gnorm0 = dnrm2_(&n, g, &inc); - delete [] w0; + // the vector M has not been used yet, and has been filled + // with zeros, so we have M == x_0 == 0 here. + fun_obj->fun(M.data()); + fun_obj->grad(M.data(), g.data()); + + double gnorm0 = dnrm2_(&n, g.data(), &inc); + double f = fun_obj->fun(w); info("init f %5.3e\n", f); - fun_obj->grad(w, g); - double gnorm = dnrm2_(&n, g, &inc); + fun_obj->grad(w, g.data()); + double gnorm = dnrm2_(&n, g.data(), &inc); if (gnorm <= eps*gnorm0) - search = 0; + return; + - double *w_new = new double[n]; - for (int iter=1; iter <= max_iter && search; ) + for (int iter=1; iter <= max_iter; ) { - fun_obj->get_diag_preconditioner(M); + fun_obj->get_diag_preconditioner(M.data()); for(int i=0; ilinesearch_and_update(w, s, & f, g, init_step_size); + double step_size = fun_obj->linesearch_and_update(w, s.data(), &f, g.data(), init_step_size); if (step_size == 0) { @@ -143,9 +143,9 @@ void NEWTON::newton(double *w) double actred = fold - f; iter++; - fun_obj->grad(w, g); + fun_obj->grad(w, g.data()); - gnorm = dnrm2_(&n, g, &inc); + gnorm = dnrm2_(&n, g.data(), &inc); if (gnorm <= eps*gnorm0) break; if (f < -1.0e+32) @@ -159,12 +159,6 @@ void NEWTON::newton(double *w) break; } } - - delete[] g; - delete[] r; - delete[] w_new; - delete[] s; - delete[] M; } int NEWTON::pcg(double *g, double *M, double *s, double *r) @@ -172,10 +166,10 @@ int NEWTON::pcg(double *g, double *M, double *s, double *r) int i, inc = 1; int n = fun_obj->get_nr_variable(); double one = 1; - double *d = new double[n]; - double *Hd = new double[n]; + std::vector d(n); + std::vector Hd(n); + std::vector z(n); double zTr, znewTrnew, alpha, beta, cgtol; - double *z = new double[n]; double Q = 0, newQ, Qdiff; for (i=0; iHv(d, Hd); + fun_obj->Hv(d.data(), Hd.data()); - alpha = zTr/ddot_(&n, d, &inc, Hd, &inc); - daxpy_(&n, &alpha, d, &inc, s, &inc); + alpha = zTr/ddot_(&n, d.data(), &inc, Hd.data(), &inc); + daxpy_(&n, &alpha, d.data(), &inc, s, &inc); alpha = -alpha; - daxpy_(&n, &alpha, Hd, &inc, r, &inc); + daxpy_(&n, &alpha, Hd.data(), &inc, r, &inc); // Using quadratic approximation as CG stopping criterion newQ = -0.5*(ddot_(&n, s, &inc, r, &inc) - ddot_(&n, s, &inc, g, &inc)); @@ -213,27 +207,23 @@ int NEWTON::pcg(double *g, double *M, double *s, double *r) else { info("WARNING: quadratic approximation > 0 or increasing in CG\n"); - break; + return cg_iter; } Q = newQ; for (i=0; i Date: Fri, 30 Oct 2020 14:56:56 +0200 Subject: [PATCH 3/4] convert indents to tabs --- newton.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/newton.cpp b/newton.cpp index b0165ab..b17ad57 100644 --- a/newton.cpp +++ b/newton.cpp @@ -97,18 +97,18 @@ void NEWTON::newton(double *w) int n = fun_obj->get_nr_variable(); double init_step_size = 1; int inc = 1; - std::vector s(n); - std::vector r(n); - std::vector g(n); + std::vector s(n); + std::vector r(n); + std::vector g(n); const double alpha_pcg = 0.01; - std::vector M(n, 0.0); + std::vector M(n, 0.0); // calculate gradient norm at w=0 for stopping condition. // the vector M has not been used yet, and has been filled // with zeros, so we have M == x_0 == 0 here. - fun_obj->fun(M.data()); - fun_obj->grad(M.data(), g.data()); + fun_obj->fun(M.data()); + fun_obj->grad(M.data(), g.data()); double gnorm0 = dnrm2_(&n, g.data(), &inc); @@ -130,7 +130,7 @@ void NEWTON::newton(double *w) int cg_iter = pcg(g.data(), M.data(), s.data(), r.data()); double fold = f; - double step_size = fun_obj->linesearch_and_update(w, s.data(), &f, g.data(), init_step_size); + double step_size = fun_obj->linesearch_and_update(w, s.data(), &f, g.data(), init_step_size); if (step_size == 0) { From 27b2629c24dfd151f137ff02ab22f78e323215f5 Mon Sep 17 00:00:00 2001 From: Erik Schultheis Date: Fri, 30 Oct 2020 14:59:01 +0200 Subject: [PATCH 4/4] simplify function returns --- newton.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/newton.cpp b/newton.cpp index b17ad57..011cfa5 100644 --- a/newton.cpp +++ b/newton.cpp @@ -202,7 +202,7 @@ int NEWTON::pcg(double *g, double *M, double *s, double *r) if (newQ <= 0 && Qdiff <= 0) { if (cg_iter * Qdiff >= cgtol * newQ) - break; + return cg_iter; } else { @@ -220,8 +220,7 @@ int NEWTON::pcg(double *g, double *M, double *s, double *r) zTr = znewTrnew; } - if (cg_iter == max_cg_iter) - info("WARNING: reaching maximal number of CG steps\n"); + info("WARNING: reaching maximal number of CG steps\n"); return cg_iter; }