Skip to content

remove remaining eigen to gsl conversions #40331

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
SilkeSchomann merged 15 commits into from
Jan 15, 2026
Merged

remove remaining eigen to gsl conversions #40331

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
7ff30c0
remove gsl from LMMD
MialLewis Nov 12, 2025
63871d1
switch out gsl for eigen in LST
MialLewis Nov 12, 2025
7d54e0c
add eigen covar_from_jacobian
MialLewis Nov 13, 2025
d32c3f3
remove gsl to eigen conversion from more files
MialLewis Nov 13, 2025
930610d
remove eigen to gsl helper classes
MialLewis Nov 13, 2025
368fc19
resolve cppcheck suppression
MialLewis Nov 13, 2025
c2454d6
resolve another cpp suppression
MialLewis Nov 13, 2025
79483e3
move new fn to own header
MialLewis Jan 5, 2026
f9d4c88
update eigen fn impl and add equivalence tests
MialLewis Jan 7, 2026
b567e40
ammend doyxgen
MialLewis Jan 7, 2026
9df6cb7
add DLL export
MialLewis Jan 7, 2026
425fa49
clarify comments and simplify eigen logic (at the slight cost of a copy)
MialLewis Jan 7, 2026
5ef5ea0
update reference file
MialLewis Jan 8, 2026
000e4df
change reference file format back to GSA
MialLewis Jan 9, 2026
1a43872
resolve gcc compiler warning
MialLewis Jan 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion Framework/CurveFitting/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ set(SRC_FILES
src/FitMW.cpp
src/EigenComplexMatrix.cpp
src/EigenComplexVector.cpp
src/EigenFunctions.cpp
src/EigenMatrix.cpp
src/EigenMatrixView.cpp
src/EigenVector.cpp
Expand Down Expand Up @@ -347,6 +348,7 @@ set(INC_FILES
inc/MantidCurveFitting/Functions/VesuvioResolution.h
inc/MantidCurveFitting/Functions/Voigt.h
inc/MantidCurveFitting/Functions/ConvTempCorrection.h
inc/MantidCurveFitting/EigenFunctions.h
inc/MantidCurveFitting/GSLFunctions.h
inc/MantidCurveFitting/GeneralDomainCreator.h
inc/MantidCurveFitting/HalfComplex.h
Expand Down Expand Up @@ -409,9 +411,9 @@ set(TEST_FILES
EigenComplexVectorTest.h
EigenFortranMatrixTest.h
EigenFortranVectorTest.h
EigenFunctionsTest.h
EigenJacobianTest.h
EigenMatrixTest.h
EigenToGSLTest.h
EigenVectorTest.h
EigenViewTest.h
FitMWTest.h
Expand Down
15 changes: 15 additions & 0 deletions Framework/CurveFitting/inc/MantidCurveFitting/EigenFunctions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2026 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once

#include "MantidCurveFitting/EigenMatrix.h"

namespace Mantid::CurveFitting {

MANTID_CURVEFITTING_DLL Eigen::MatrixXd covar_from_jacobian(const map_type &J, double epsrel);

} // namespace Mantid::CurveFitting
14 changes: 0 additions & 14 deletions Framework/CurveFitting/inc/MantidCurveFitting/GSLFunctions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,19 +50,5 @@ int gsl_f(const gsl_vector *x, void *params, gsl_vector *f);
int gsl_df(const gsl_vector *x, void *params, gsl_matrix *J);
int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J);

/// take data from Eigen Vector and take a gsl view
inline gsl_vector_view getGSLVectorView(vec_map_type &v) { return gsl_vector_view_array(v.data(), v.size()); }
/// take data from an Eigen Matrix and return a transposed a gsl view.
inline gsl_matrix_view getGSLMatrixView(map_type &tr) { return gsl_matrix_view_array(tr.data(), tr.cols(), tr.rows()); }

/// take const data from Eigen Vector and take a gsl view
inline gsl_vector_const_view const getGSLVectorView_const(const vec_map_type v) {
return gsl_vector_const_view_array(v.data(), v.size());
}
/// take data from a constEigen Matrix and return a transposed gsl view.
inline gsl_matrix_const_view const getGSLMatrixView_const(const map_type m) {
return gsl_matrix_const_view_array(m.data(), m.cols(), m.rows());
}

} // namespace CurveFitting
} // namespace Mantid
14 changes: 3 additions & 11 deletions Framework/CurveFitting/src/CostFunctions/CostFuncFitting.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,11 @@
//----------------------------------------------------------------------
#include "MantidCurveFitting/CostFunctions/CostFuncFitting.h"
#include "MantidAPI/IConstraint.h"
#include "MantidCurveFitting/EigenFunctions.h"
#include "MantidCurveFitting/EigenJacobian.h"
#include "MantidCurveFitting/GSLFunctions.h"
#include "MantidCurveFitting/SeqDomain.h"
#include "MantidKernel/Exception.h"

#include <gsl/gsl_multifit_nlin.h>
#include <limits>
#include <utility>

Expand Down Expand Up @@ -132,15 +131,8 @@ void CostFuncFitting::calActiveCovarianceMatrix(EigenMatrix &covar, double epsre
// calculate the derivatives
m_function->functionDeriv(*m_domain, J);

// let the GSL to compute the covariance matrix
EigenMatrix covarTr(covar.tr());
EigenMatrix tempJTr;
tempJTr = j.transpose();
gsl_matrix_const_view tempJTrGSL = getGSLMatrixView_const(tempJTr.inspector());
gsl_matrix_view covarTrGSL = getGSLMatrixView(covarTr.mutator());

gsl_multifit_covar(&tempJTrGSL.matrix, epsrel, &covarTrGSL.matrix);
covar = covarTr.tr();
// let Eigen compute the covariance matrix
covar = covar_from_jacobian(j, epsrel);
}

/** Calculates covariance matrix
Expand Down
70 changes: 70 additions & 0 deletions Framework/CurveFitting/src/EigenFunctions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 2026 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/EigenFunctions.h"

namespace Mantid::CurveFitting {

/** Mimics gsl_multifit_covar(J, epsrel, covar)
* @param J :: Jacobian
* @param epsrel :: relative threshold to decide rank deficiency
* @return covariance matrix, rows/cols set to zero for dependent params
*/
Eigen::MatrixXd covar_from_jacobian(const map_type &J, double epsrel) {
if (epsrel < 0.0) {
throw std::invalid_argument("epsrel must be non-negative");
}

const Eigen::Index nr = J.rows(); // nr = num rows
const Eigen::Index nc = J.cols(); // nc = num cols
if ((nc == 0) || (nr == 0))
return Eigen::MatrixXd{};

// Pivoted QR Decomposition: JP = QR
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(J);

// R is (nc x np) in general; the useful part for J^T J is (nc x nc) block
Eigen::MatrixXd R = qr.matrixR().topLeftCorner(nc, nc).template triangularView<Eigen::Upper>().toDenseMatrix();

// Determine rank using the same rule as GSL: compare diag entries to |R_11|
const double r11 = std::abs(R(0, 0));
Eigen::Index rank = 0;
if (r11 > 0.0) {
for (Eigen::Index k = 0; k < nc; ++k) {
if (std::abs(R(k, k)) > epsrel * r11) { // column considered linerally independant
++rank;
} else {
break; // with pivoting, following first dependent subsequent cols should also be dependant
}
}
} else {
// R11 == 0, everything dependent
rank = 0;
}

// Build covariance in the pivoted parameter order
Eigen::MatrixXd cov_pivot = Eigen::MatrixXd::Zero(nc, nc);

if (rank > 0) {
// cov = (R^T R)^{-1} for the independent cols = R^{-1} R^{-T}
const Eigen::MatrixXd R1 = R.topLeftCorner(rank, rank).template triangularView<Eigen::Upper>();
Eigen::MatrixXd invR1 = R1.inverse();
Eigen::MatrixXd cov1 = invR1 * invR1.transpose();
cov_pivot.topLeftCorner(rank, rank) = cov1;
}

// Unpivot back to original parameter order:
// cov = (J^T J)^{-1} = P (R^T R)^{-1} P^T
const Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> P = qr.colsPermutation();
Eigen::MatrixXd cov = P * cov_pivot * P.transpose();

return cov;
}

} // namespace Mantid::CurveFitting
5 changes: 1 addition & 4 deletions Framework/CurveFitting/src/FuncMinimizers/DampedGaussNewtonMinimizer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@
#include "MantidKernel/Logger.h"

#include <cmath>
#include <gsl/gsl_blas.h>

namespace Mantid::CurveFitting::FuncMinimisers {

Expand Down Expand Up @@ -121,9 +120,7 @@ bool DampedGaussNewtonMinimizer::iterate(size_t /*iteration*/) {
// Try the stop condition
EigenVector p(n);
m_leastSquares->getParameters(p);
gsl_vector_view dx_gsl = getGSLVectorView(dx.mutator());
double dx_norm = gsl_blas_dnrm2(&dx_gsl.vector);
return dx_norm >= m_relTol;
return dx.norm() >= m_relTol;
}

/// Return current value of the cost function
Expand Down
15 changes: 5 additions & 10 deletions Framework/CurveFitting/src/FuncMinimizers/LevenbergMarquardtMDMinimizer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@
#include "MantidKernel/Logger.h"

#include <cmath>
#include <gsl/gsl_blas.h>

namespace Mantid::CurveFitting::FuncMinimisers {
namespace {
Expand Down Expand Up @@ -194,15 +193,11 @@ bool LevenbergMarquardtMDMinimizer::iterate(size_t /*iteration*/) {
double dL;
// der -> - der - 0.5 * hessian * dx

EigenMatrix tempHessianTr = m_costFunction->getHessian().tr();
const gsl_matrix_const_view tempHessianTrGSL = getGSLMatrixView_const(tempHessianTr.inspector());
const gsl_vector_const_view dxGSL = getGSLVectorView_const(dx.inspector());
gsl_vector_view ddGSL = getGSLVectorView(dd.mutator());

gsl_blas_dgemv(CblasNoTrans, -0.5, &tempHessianTrGSL.matrix, &dxGSL.vector, 1., &ddGSL.vector);
EigenVector Trdx = m_costFunction->getHessian().tr() * dx;
Trdx *= -0.5;
dd += Trdx;
// calculate the linear part of the change in cost function
// dL = - der * dx - 0.5 * dx * hessian * dx
gsl_blas_ddot(&ddGSL.vector, &dxGSL.vector, &dL);
dL = dd.dot(dx);

double F1 = m_costFunction->val();
if (verbose) {
Expand All @@ -216,7 +211,7 @@ bool LevenbergMarquardtMDMinimizer::iterate(size_t /*iteration*/) {
if (m_rho >= 0) {
EigenVector p(n);
m_costFunction->getParameters(p);
double dx_norm = gsl_blas_dnrm2(&dxGSL.vector);
double dx_norm = dx.norm();
if (dx_norm < absError) {
if (verbose) {
g_log.warning() << "Successful fit, parameters changed by less than " << absError << '\n';
Expand Down
38 changes: 15 additions & 23 deletions Framework/CurveFitting/src/Functions/ChebfunBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,11 @@
#include "MantidCurveFitting/GSLFunctions.h"
#include "MantidCurveFitting/HalfComplex.h"

#include <gsl/gsl_eigen.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_real.h>

#include <Eigen/Eigenvalues>

#include <algorithm>
#include <cassert>
#include <cmath>
Expand Down Expand Up @@ -684,21 +684,16 @@ std::vector<double> ChebfunBase::roots(const std::vector<double> &a) const {
C.set(N + i, lasti, tmp);
}

gsl_vector_complex *eigenval = gsl_vector_complex_alloc(N2);
auto workspace = gsl_eigen_nonsymm_alloc(N2);

EigenMatrix C_tr(C.tr());
gsl_matrix_view C_tr_gsl = getGSLMatrixView(C_tr.mutator());
gsl_eigen_nonsymm(&C_tr_gsl.matrix, eigenval, workspace);
gsl_eigen_nonsymm_free(workspace);
Eigen::EigenSolver<Eigen::MatrixXd> es(C.mutator());
const Eigen::VectorXcd evals = es.eigenvalues();

const double Dx = endX() - startX();
bool isFirst = true;
double firstIm = 0;
for (size_t i = 0; i < N2; ++i) {
auto val = gsl_vector_complex_get(eigenval, i);
double re = GSL_REAL(val);
double im = GSL_IMAG(val);
const std::complex<double> val = evals[static_cast<Eigen::Index>(i)];
double re = val.real();
double im = val.imag();
double ab = re * re + im * im;
if (fabs(ab - 1.0) > 1e-2) {
isFirst = true;
Expand All @@ -715,8 +710,6 @@ std::vector<double> ChebfunBase::roots(const std::vector<double> &a) const {
isFirst = true;
}
}
gsl_vector_complex_free(eigenval);

return r;
}

Expand Down Expand Up @@ -824,25 +817,24 @@ std::vector<double> ChebfunBase::smooth(const std::vector<double> &xvalues, cons
return y;
}

// high frequency filter values: smooth decreasing function
auto ri0f = static_cast<double>(i0 + 1);
double xm = (1.0 + ri0f) / 2;
ym /= ri0f;
double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
double b1 = ym - a1 * xm;

// std::cerr << "(a1,b1) = (" << a1 << ',' << b1 << ')' << '\n';

// calculate coeffs of a cubic c3*i^3 + c2*i^2 + c1*i + c0
// which will replace the linear a1*i + b1 in building the
// second part of the filter
double c0, c1, c2, c3;
{
// high frequency filter values: smooth decreasing function
auto ri0f = static_cast<double>(i0 + 1);
double xm = (1.0 + ri0f) / 2;
ym /= ri0f;

auto x0 = double(i0 + 1);
auto x1 = double(n + 1);
double sigma = g_tolerance / noise / 10;
double s = sigma / (1.0 - sigma);
double m1 = log(s);
double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
double b1 = ym - a1 * xm;
// std::cerr << "(a1,b1) = (" << a1 << ',' << b1 << ')' << '\n';
double m0 = a1 * x0 + b1;
if (a1 < 0.0) {
c3 = 0.0;
Expand Down
13 changes: 4 additions & 9 deletions Framework/CurveFitting/test/CostFunctions/LeastSquaresTest.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@
#include "MantidCurveFitting/Functions/UserFunction.h"
#include "MantidCurveFitting/GSLFunctions.h"

#include <gsl/gsl_blas.h>
#include <sstream>

using namespace Mantid;
Expand Down Expand Up @@ -201,15 +200,11 @@ class LeastSquaresTest : public CxxTest::TestSuite {
dx.set(1, -0.2);

double L; // = d*dx + 0.5 * dx * H * dx
EigenVector Trdx = H.tr() * dx;
Trdx *= 0.5;
g += Trdx;
L = g.dot(dx);

EigenMatrix temp_H_tr = H.tr();
const gsl_matrix_const_view temp_H_tr_gsl = getGSLMatrixView_const(temp_H_tr.inspector());
const gsl_vector_const_view dx_gsl = getGSLVectorView_const(dx.inspector());
gsl_vector_view g_gsl = getGSLVectorView(g.mutator());

gsl_blas_dgemv(CblasNoTrans, 0.5, &temp_H_tr_gsl.matrix, &dx_gsl.vector, 1., &g_gsl.vector);

gsl_blas_ddot(&g_gsl.vector, &dx_gsl.vector, &L);
TS_ASSERT_DELTA(L, -0.145, 1e-10); // L + costFun->val() == 0
}

Expand Down
Loading