From 0afbd390d6e9202043c7af7cae3b35c3c7c93a4b Mon Sep 17 00:00:00 2001
From: Karl Palmen <karl.palmen@stfc.ac.uk>
Date: Mon, 30 Oct 2017 15:46:20 +0000
Subject: [PATCH] Remove calculateStep workspace argument re #16748
Signed-off-by: Karl Palmen <karl.palmen@stfc.ac.uk>
---
.../FuncMinimizers/DTRSMinimizer.h | 14 +-
.../FuncMinimizers/MoreSorensenMinimizer.h | 17 ++-
.../FuncMinimizers/TrustRegionMinimizer.h | 3 +-
.../MantidCurveFitting/RalNlls/Workspaces.h | 27 +---
.../src/FuncMinimizers/DTRSMinimizer.cpp | 48 ++++---
.../FuncMinimizers/MoreSorensenMinimizer.cpp | 121 +++++++++---------
.../FuncMinimizers/TrustRegionMinimizer.cpp | 3 +-
7 files changed, 113 insertions(+), 120 deletions(-)
diff --git a/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/DTRSMinimizer.h b/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/DTRSMinimizer.h
index b3f39a763dc..2b07e383166 100644
--- a/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/DTRSMinimizer.h
+++ b/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/DTRSMinimizer.h
@@ -42,8 +42,18 @@ private:
const DoubleFortranMatrix &hf,
const DoubleFortranVector &g, double Delta,
DoubleFortranVector &d, double &normd,
- const NLLS::nlls_options &options,
- NLLS::calculate_step_work &w) override;
+ const NLLS::nlls_options &options) override;
+
+ void solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
+ const DoubleFortranMatrix &hf, double Delta,
+ DoubleFortranVector &d, double &normd,
+ const NLLS::nlls_options &options);
+
+ // Used for calculating step
+ DoubleFortranMatrix m_A, m_ev;
+ DoubleFortranVector m_ew, m_v, m_v_trans, m_d_trans;
+ NLLS::all_eig_symm_work m_all_eig_symm_ws;
+ NLLS::apply_scaling_work m_apply_scaling_ws;
};
} // namespace FuncMinimisers
diff --git a/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/MoreSorensenMinimizer.h b/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/MoreSorensenMinimizer.h
index a21b05c94fd..0b101a8c9d4 100644
--- a/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/MoreSorensenMinimizer.h
+++ b/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/MoreSorensenMinimizer.h
@@ -45,8 +45,21 @@ private:
const DoubleFortranMatrix &hf,
const DoubleFortranVector &g, double Delta,
DoubleFortranVector &d, double &normd,
- const NLLS::nlls_options &options,
- NLLS::calculate_step_work &w) override;
+ const NLLS::nlls_options &options) override;
+
+ void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
+ const DoubleFortranMatrix &hf, double Delta,
+ DoubleFortranVector &d, double &nd,
+ const NLLS::nlls_options &options);
+
+ bool getPdShift(double &sigma, DoubleFortranVector &d,
+ const NLLS::nlls_options &options);
+
+ // Used for calculating step
+ DoubleFortranMatrix m_A, m_LtL, m_AplusSigma;
+ DoubleFortranVector m_v, m_q, m_y1;
+ NLLS::min_eig_symm_work m_min_eig_symm_ws;
+ NLLS::apply_scaling_work m_apply_scaling_ws;
};
} // namespace FuncMinimisers
diff --git a/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/TrustRegionMinimizer.h b/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/TrustRegionMinimizer.h
index 047be9f88d7..4b1d3453c1c 100644
--- a/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/TrustRegionMinimizer.h
+++ b/Framework/CurveFitting/inc/MantidCurveFitting/FuncMinimizers/TrustRegionMinimizer.h
@@ -62,8 +62,7 @@ private:
calculateStep(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
const DoubleFortranMatrix &hf, const DoubleFortranVector &g,
double Delta, DoubleFortranVector &d, double &normd,
- const NLLS::nlls_options &options,
- NLLS::calculate_step_work &w) = 0;
+ const NLLS::nlls_options &options) = 0;
/// Stored cost function
boost::shared_ptr<CostFunctions::CostFuncLeastSquares> m_leastSquares;
diff --git a/Framework/CurveFitting/inc/MantidCurveFitting/RalNlls/Workspaces.h b/Framework/CurveFitting/inc/MantidCurveFitting/RalNlls/Workspaces.h
index 275c06eaa12..eff98a6bf6f 100644
--- a/Framework/CurveFitting/inc/MantidCurveFitting/RalNlls/Workspaces.h
+++ b/Framework/CurveFitting/inc/MantidCurveFitting/RalNlls/Workspaces.h
@@ -240,31 +240,6 @@ struct all_eig_symm_work {
/// workspace for subroutine applyScaling
struct apply_scaling_work {
DoubleFortranVector diag;
- DoubleFortranMatrix ev;
- DoubleFortranVector tempvec;
- all_eig_symm_work all_eig_symm_ws;
-};
-
-/// workspace for subroutine dtrs_work
-struct solve_dtrs_work {
- DoubleFortranMatrix A, ev;
- DoubleFortranVector ew, v, v_trans, d_trans;
- all_eig_symm_work all_eig_symm_ws;
- apply_scaling_work apply_scaling_ws;
-};
-
-/// workspace for subroutine moreSorensen
-struct more_sorensen_work {
- DoubleFortranMatrix A, LtL, AplusSigma;
- DoubleFortranVector v, q, y1;
- min_eig_symm_work min_eig_symm_ws;
- apply_scaling_work apply_scaling_ws;
-};
-
-/// workspace for subroutine calculateStep
-struct calculate_step_work {
- more_sorensen_work more_sorensen_ws;
- solve_dtrs_work solve_dtrs_ws;
};
/// workspace for subroutine getSvdJ
@@ -295,7 +270,7 @@ struct NLLS_workspace {
DoubleFortranVector resvec, gradvec;
DoubleFortranVector largest_sv, smallest_sv;
get_svd_J_work get_svd_J_ws;
- calculate_step_work calculate_step_ws;
+ //calculate_step_work calculate_step_ws;
evaluate_model_work evaluate_model_ws;
NLLS_workspace();
void initialize(int n, int m, const nlls_options &options);
diff --git a/Framework/CurveFitting/src/FuncMinimizers/DTRSMinimizer.cpp b/Framework/CurveFitting/src/FuncMinimizers/DTRSMinimizer.cpp
index 04fa69f1dd2..21e0a49ffbd 100644
--- a/Framework/CurveFitting/src/FuncMinimizers/DTRSMinimizer.cpp
+++ b/Framework/CurveFitting/src/FuncMinimizers/DTRSMinimizer.cpp
@@ -885,6 +885,8 @@ void dtrsSolve(int n, double radius, double f, const DoubleFortranVector &c,
}
}
+} // namespace
+
/** Solve the trust-region subproblem using
* the DTRS method from Galahad
*
@@ -903,11 +905,10 @@ void dtrsSolve(int n, double radius, double f, const DoubleFortranVector &c,
* @param inform :: The inform struct.
* @param w :: The work struct.
*/
-void solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
- const DoubleFortranMatrix &hf, double Delta,
- DoubleFortranVector &d, double &normd,
- const NLLS::nlls_options &options,
- NLLS::solve_dtrs_work &w) {
+void DTRSMinimizer::solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
+ const DoubleFortranMatrix &hf, double Delta,
+ DoubleFortranVector &d, double &normd,
+ const NLLS::nlls_options &options) {
dtrs_control_type dtrs_options;
dtrs_inform_type dtrs_inform;
@@ -925,23 +926,23 @@ void solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
//
// first, find the matrix H and vector v
// Set A = J^T J
- NLLS::matmultInner(J, w.A);
+ NLLS::matmultInner(J, m_A);
// add any second order information...
// so A = J^T J + HF
- w.A += hf;
+ m_A += hf;
// now form v = J^T f
- NLLS::multJt(J, f, w.v);
+ NLLS::multJt(J, f, m_v);
// if scaling needed, do it
if (options.scale != 0) {
- applyScaling(J, w.A, w.v, w.apply_scaling_ws, options);
+ applyScaling(J, m_A, m_v, m_apply_scaling_ws, options);
}
// Now that we have the unprocessed matrices, we need to get an
// eigendecomposition to make A diagonal
//
- NLLS::allEigSymm(w.A, w.ew, w.ev);
+ NLLS::allEigSymm(m_A, m_ew, m_ev);
// We can now change variables, setting y = Vp, getting
// Vd = arg min_(Vx) v^T p + 0.5 * (Vp)^T D (Vp)
@@ -951,51 +952,48 @@ void solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
// s.t. ||x|| \leq Delta
// <=>
// we need to get the transformed vector v
- NLLS::multJt(w.ev, w.v, w.v_trans);
+ NLLS::multJt(m_ev, m_v, m_v_trans);
// we've now got the vectors we need, pass to dtrsSolve
dtrsInitialize(dtrs_options);
auto n = J.len2();
- if (w.v_trans.len() != n) {
- w.v_trans.allocate(n);
+ if (m_v_trans.len() != n) {
+ m_v_trans.allocate(n);
}
for (int ii = 1; ii <= n; ++ii) { // for_do(ii, 1,n)
- if (fabs(w.v_trans(ii)) < EPSILON_MCH) {
- w.v_trans(ii) = ZERO;
+ if (fabs(m_v_trans(ii)) < EPSILON_MCH) {
+ m_v_trans(ii) = ZERO;
}
- if (fabs(w.ew(ii)) < EPSILON_MCH) {
- w.ew(ii) = ZERO;
+ if (fabs(m_ew(ii)) < EPSILON_MCH) {
+ m_ew(ii) = ZERO;
}
}
- dtrsSolve(n, Delta, ZERO, w.v_trans, w.ew, w.d_trans, dtrs_options,
+ dtrsSolve(n, Delta, ZERO, m_v_trans, m_ew, m_d_trans, dtrs_options,
dtrs_inform);
// and return the un-transformed vector
- NLLS::multJ(w.ev, w.d_trans, d);
+ NLLS::multJ(m_ev, m_d_trans, d);
normd = NLLS::norm2(d); // ! ||d||_D
if (options.scale != 0) {
for (int ii = 1; ii <= n; ++ii) { // for_do(ii, 1, n)
- d(ii) = d(ii) / w.apply_scaling_ws.diag(ii);
+ d(ii) = d(ii) / m_apply_scaling_ws.diag(ii);
}
}
} // solveDtrs
-} // namespace
-
/** Implements the abstarct method of TrustRegionMinimizer.
*/
void DTRSMinimizer::calculateStep(
const DoubleFortranMatrix &J, const DoubleFortranVector &f,
const DoubleFortranMatrix &hf, const DoubleFortranVector &, double Delta,
- DoubleFortranVector &d, double &normd, const NLLS::nlls_options &options,
- NLLS::calculate_step_work &w) {
- solveDtrs(J, f, hf, Delta, d, normd, options, w.solve_dtrs_ws);
+ DoubleFortranVector &d, double &normd, const NLLS::nlls_options &options) {
+ solveDtrs(J, f, hf, Delta, d, normd, options);
}
} // namespace FuncMinimisers
diff --git a/Framework/CurveFitting/src/FuncMinimizers/MoreSorensenMinimizer.cpp b/Framework/CurveFitting/src/FuncMinimizers/MoreSorensenMinimizer.cpp
index 429ddb78933..9106e71a93e 100644
--- a/Framework/CurveFitting/src/FuncMinimizers/MoreSorensenMinimizer.cpp
+++ b/Framework/CurveFitting/src/FuncMinimizers/MoreSorensenMinimizer.cpp
@@ -93,35 +93,6 @@ DoubleFortranVector negative(const DoubleFortranVector &v) {
return neg;
}
-/** Given an indefinite matrix w.A, find a shift sigma
- * such that (A + sigma I) is positive definite.
- * @param sigma :: The result (shift).
- * @param d :: A solution vector to the system of linear equations
- * with the found positive defimnite matrix. The RHS vector is -w.v.
- * @param options :: The options.
- * @param inform :: The inform struct.
- * @param w :: The work struct.
- * @return true if successful
- */
-bool getPdShift(double &sigma, DoubleFortranVector &d,
- const NLLS::nlls_options &options,
- NLLS::more_sorensen_work &w) {
- int no_shifts = 0;
- bool successful_shift = false;
- while (!successful_shift) {
- shiftMatrix(w.A, sigma, w.AplusSigma);
- successful_shift = solveSpd(w.AplusSigma, negative(w.v), w.LtL, d);
- if (!successful_shift) {
- // We try again with a shifted sigma, but no too many times.
- no_shifts = no_shifts + 1;
- if (no_shifts == 10) {
- return false;
- }
- sigma = sigma + (pow(10.0, no_shifts)) * options.more_sorensen_shift;
- }
- }
- return true;
-}
/** A subroutine to find the optimal beta such that
* || d || = Delta, where d = a + beta * b
@@ -156,6 +127,37 @@ bool findBeta(const DoubleFortranVector &a, const DoubleFortranVector &b,
return true;
}
+} // namespace
+
+ /** Given an indefinite matrix m_A, find a shift sigma
+ * such that (A + sigma I) is positive definite.
+ * @param sigma :: The result (shift).
+ * @param d :: A solution vector to the system of linear equations
+ * with the found positive defimnite matrix. The RHS vector is -m_v.
+ * @param options :: The options.
+ * @param inform :: The inform struct.
+ * @param w :: The work struct.
+ * @return true if successful
+ */
+bool MoreSorensenMinimizer::getPdShift(double &sigma, DoubleFortranVector &d,
+ const NLLS::nlls_options &options) {
+ int no_shifts = 0;
+ bool successful_shift = false;
+ while (!successful_shift) {
+ shiftMatrix(m_A, sigma, m_AplusSigma);
+ successful_shift = solveSpd(m_AplusSigma, negative(m_v), m_LtL, d);
+ if (!successful_shift) {
+ // We try again with a shifted sigma, but no too many times.
+ no_shifts = no_shifts + 1;
+ if (no_shifts == 10) {
+ return false;
+ }
+ sigma = sigma + (pow(10.0, no_shifts)) * options.more_sorensen_shift;
+ }
+ }
+ return true;
+}
+
/** Solve the trust-region subproblem using
* the method of More and Sorensen
*
@@ -174,11 +176,10 @@ bool findBeta(const DoubleFortranVector &a, const DoubleFortranVector &b,
* @param inform :: The inform struct.
* @param w :: The work struct.
*/
-void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
+void MoreSorensenMinimizer::moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
const DoubleFortranMatrix &hf, double Delta,
DoubleFortranVector &d, double &nd,
- const NLLS::nlls_options &options,
- NLLS::more_sorensen_work &w) {
+ const NLLS::nlls_options &options) {
// The code finds
// d = arg min_p v^T p + 0.5 * p^T A p
@@ -187,42 +188,42 @@ void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
// set A and v for the model being considered here...
// Set A = J^T J
- NLLS::matmultInner(J, w.A);
+ NLLS::matmultInner(J, m_A);
// add any second order information...
// so A = J^T J + HF
- w.A += hf;
+ m_A += hf;
// now form v = J^T f
- NLLS::multJt(J, f, w.v);
+ NLLS::multJt(J, f, m_v);
// if scaling needed, do it
if (options.scale != 0) {
- applyScaling(J, w.A, w.v, w.apply_scaling_ws, options);
+ applyScaling(J, m_A, m_v, m_apply_scaling_ws, options);
}
auto n = J.len2();
- auto scaleBack = [n, &d, &options, &w]() {
+ auto scaleBack = [n, &d, &options, this]() {
if (options.scale != 0) {
for (int i = 1; i <= n; ++i) {
- d(i) = d(i) / w.apply_scaling_ws.diag(i);
+ d(i) = d(i) / m_apply_scaling_ws.diag(i);
}
}
};
auto local_ms_shift = options.more_sorensen_shift;
// d = -A\v
- DoubleFortranVector negv = w.v;
+ DoubleFortranVector negv = m_v;
negv *= -1.0;
- bool matrix_ok = solveSpd(w.A, negv, w.LtL, d);
+ bool matrix_ok = solveSpd(m_A, negv, m_LtL, d);
double sigma = 0.0;
if (matrix_ok) {
// A is symmetric positive definite....
sigma = NLLS::ZERO;
} else {
// shift and try again
- minEigSymm(w.A, sigma, w.y1);
+ minEigSymm(m_A, sigma, m_y1);
sigma = -(sigma - local_ms_shift);
// find a shift that makes (A + sigma I) positive definite
- bool ok = getPdShift(sigma, d, options, w);
+ bool ok = getPdShift(sigma, d, options);
if (!ok) {
scaleBack();
return;
@@ -250,10 +251,10 @@ void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
// we're good....exit
break;
}
- if (w.y1.len() == n) {
+ if (m_y1.len() == n) {
double alpha = 0.0;
- if(findBeta(d, w.y1, Delta, alpha)) {
- DoubleFortranVector tmp = w.y1;
+ if(findBeta(d, m_y1, Delta, alpha)) {
+ DoubleFortranVector tmp = m_y1;
tmp *= alpha;
d += tmp;
}
@@ -262,22 +263,22 @@ void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
break;
}
- // w.q = R'\d
- // DTRSM( "Left", "Lower", "No Transpose", "Non-unit", n, 1, one, w.LtL, n,
- // w.q, n );
- for (int j = 1; j <= w.LtL.len1(); ++j) {
- for (int k = j + 1; k <= w.LtL.len1(); ++k) {
- w.LtL(j, k) = 0.0;
+ // m_q = R'\d
+ // DTRSM( "Left", "Lower", "No Transpose", "Non-unit", n, 1, one, m_LtL, n,
+ // m_q, n );
+ for (int j = 1; j <= m_LtL.len1(); ++j) {
+ for (int k = j + 1; k <= m_LtL.len1(); ++k) {
+ m_LtL(j, k) = 0.0;
}
}
- w.LtL.solve(d, w.q);
+ m_LtL.solve(d, m_q);
- auto nq = NLLS::norm2(w.q);
+ auto nq = NLLS::norm2(m_q);
sigma_shift = (pow((nd / nq), 2)) * ((nd - Delta) / Delta);
if (fabs(sigma_shift) < options.more_sorensen_tiny * fabs(sigma)) {
if (no_restarts < 1) {
// find a shift that makes (A + sigma I) positive definite
- bool ok = getPdShift(sigma, d, options, w);
+ bool ok = getPdShift(sigma, d, options);
if (!ok) {
break;
}
@@ -291,10 +292,10 @@ void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
sigma = sigma + sigma_shift;
}
- shiftMatrix(w.A, sigma, w.AplusSigma);
- DoubleFortranVector negv = w.v;
+ shiftMatrix(m_A, sigma, m_AplusSigma);
+ DoubleFortranVector negv = m_v;
negv *= -1.0;
- bool matrix_ok = solveSpd(w.AplusSigma, negv, w.LtL, d);
+ bool matrix_ok = solveSpd(m_AplusSigma, negv, m_LtL, d);
if (!matrix_ok) {
break;
}
@@ -309,16 +310,14 @@ void moreSorensen(const DoubleFortranMatrix &J, const DoubleFortranVector &f,
scaleBack();
}
-} // namespace
/** Implements the abstract method of TrustRegionMinimizer.
*/
void MoreSorensenMinimizer::calculateStep(
const DoubleFortranMatrix &J, const DoubleFortranVector &f,
const DoubleFortranMatrix &hf, const DoubleFortranVector &, double Delta,
- DoubleFortranVector &d, double &normd, const NLLS::nlls_options &options,
- NLLS::calculate_step_work &w) {
- moreSorensen(J, f, hf, Delta, d, normd, options, w.more_sorensen_ws);
+ DoubleFortranVector &d, double &normd, const NLLS::nlls_options &options) {
+ moreSorensen(J, f, hf, Delta, d, normd, options);
}
} // namespace FuncMinimisers
diff --git a/Framework/CurveFitting/src/FuncMinimizers/TrustRegionMinimizer.cpp b/Framework/CurveFitting/src/FuncMinimizers/TrustRegionMinimizer.cpp
index 1d5b487058b..5971e11a567 100644
--- a/Framework/CurveFitting/src/FuncMinimizers/TrustRegionMinimizer.cpp
+++ b/Framework/CurveFitting/src/FuncMinimizers/TrustRegionMinimizer.cpp
@@ -241,8 +241,7 @@ bool TrustRegionMinimizer::iterate(size_t) {
return true;
}
// Calculate the step d that the model thinks we should take next
- calculateStep(w.J, w.f, w.hf, w.g, w.Delta, w.d, w.normd, options,
- w.calculate_step_ws);
+ calculateStep(w.J, w.f, w.hf, w.g, w.Delta, w.d, w.normd, options);
// Accept the step?
w.Xnew = X;
--
GitLab