Skip to content
Snippets Groups Projects
Commit d820b019 authored by Karl Palmen's avatar Karl Palmen
Browse files

Simplify Scaling re #16748


Signed-off-by: default avatarKarl Palmen <karl.palmen@stfc.ac.uk>
parent 0afbd390
No related branches found
No related tags found
No related merge requests found
......@@ -53,7 +53,7 @@ private:
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;
DoubleFortranVector m_scale;
};
} // namespace FuncMinimisers
......
......@@ -59,7 +59,7 @@ private:
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;
DoubleFortranVector m_scale;
};
} // namespace FuncMinimisers
......
......@@ -28,7 +28,7 @@ void rankOneUpdate(DoubleFortranMatrix &hf, NLLS_workspace &w);
void testConvergence(double normF, double normJF, double normF0, double normJF0,
const nlls_options &options, nlls_inform &inform);
void applyScaling(const DoubleFortranMatrix &J, DoubleFortranMatrix &A,
DoubleFortranVector &v, apply_scaling_work &w,
DoubleFortranVector &v, DoubleFortranVector &scale,
const nlls_options &options);
void allEigSymm(const DoubleFortranMatrix &A, DoubleFortranVector &ew,
DoubleFortranMatrix &ev);
......
......@@ -237,11 +237,6 @@ struct all_eig_symm_work {
DoubleFortranVector work;
};
/// workspace for subroutine applyScaling
struct apply_scaling_work {
DoubleFortranVector diag;
};
/// workspace for subroutine getSvdJ
struct get_svd_J_work {
DoubleFortranVector Jcopy, S, work;
......
......@@ -936,7 +936,7 @@ void DTRSMinimizer::solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranV
// if scaling needed, do it
if (options.scale != 0) {
applyScaling(J, m_A, m_v, m_apply_scaling_ws, options);
applyScaling(J, m_A, m_v, m_scale, options);
}
// Now that we have the unprocessed matrices, we need to get an
......@@ -981,7 +981,7 @@ void DTRSMinimizer::solveDtrs(const DoubleFortranMatrix &J, const DoubleFortranV
if (options.scale != 0) {
for (int ii = 1; ii <= n; ++ii) { // for_do(ii, 1, n)
d(ii) = d(ii) / m_apply_scaling_ws.diag(ii);
d(ii) = d(ii) / m_scale(ii);
}
}
......
......@@ -197,14 +197,14 @@ void MoreSorensenMinimizer::moreSorensen(const DoubleFortranMatrix &J, const Dou
// if scaling needed, do it
if (options.scale != 0) {
applyScaling(J, m_A, m_v, m_apply_scaling_ws, options);
applyScaling(J, m_A, m_v, m_scale, options);
}
auto n = J.len2();
auto scaleBack = [n, &d, &options, this]() {
if (options.scale != 0) {
for (int i = 1; i <= n; ++i) {
d(i) = d(i) / m_apply_scaling_ws.diag(i);
d(i) = d(i) / m_scale(i);
}
}
};
......
......@@ -289,12 +289,12 @@ void testConvergence(double normF, double normJF, double normF0, double normJF0,
/// @param options :: The options.
/// @param inform :: The information.
void applyScaling(const DoubleFortranMatrix &J, DoubleFortranMatrix &A,
DoubleFortranVector &v, apply_scaling_work &w,
DoubleFortranVector &v, DoubleFortranVector &scale,
const nlls_options &options) {
auto m = J.len1();
auto n = J.len2();
if (w.diag.len() != n) {
w.diag.allocate(n);
if (scale.len() != n) {
scale.allocate(n);
}
switch (options.scale) {
......@@ -330,9 +330,9 @@ void applyScaling(const DoubleFortranMatrix &J, DoubleFortranMatrix &A,
}
temp = sqrt(temp);
if (options.scale_require_increase) {
w.diag(ii) = std::max(temp, w.diag(ii));
scale(ii) = std::max(temp, scale(ii));
} else {
w.diag(ii) = temp;
scale(ii) = temp;
}
}
break;
......@@ -344,7 +344,7 @@ void applyScaling(const DoubleFortranMatrix &J, DoubleFortranMatrix &A,
// Now we have the w.diagonal scaling matrix, actually scale the
// Hessian approximation and J^Tf
for (int ii = 1; ii <= n; ++ii) { // for_do(ii, 1,n)
double temp = w.diag(ii);
double temp = scale(ii);
v(ii) = v(ii) / temp;
for (int jj = 1; jj <= n; ++jj) { // for_do(jj,1,n)
A(ii, jj) = A(ii, jj) / temp;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment