9
9
#include " MantidAPI/Axis.h"
10
10
#include " MantidAPI/ITableWorkspace.h"
11
11
#include " MantidAPI/MatrixWorkspaceValidator.h"
12
- #include " MantidCurveFitting/EigenMatrix.h"
13
12
#include " MantidDataObjects/Workspace2D.h"
14
13
#include " MantidDataObjects/WorkspaceCreation.h"
15
14
#include " MantidHistogramData/Histogram.h"
16
15
#include " MantidKernel/PhysicalConstants.h"
17
16
#include " MantidKernel/Unit.h"
18
17
18
+ #include " Eigen/Dense"
19
+
20
+ // Use of a `long double` datatype is required on osx-arm64 to bring the precision of eigen vector-matrix multiplication
21
+ // inline with the other operating systems.
22
+ #if defined(__APPLE__) && defined(__arm64__)
23
+ typedef long double eigenDataType;
24
+ #else
25
+ typedef double eigenDataType;
26
+ #endif
27
+
19
28
using namespace Mantid ::DataObjects;
20
29
using namespace Mantid ::HistogramData;
21
30
@@ -238,7 +247,7 @@ API::MatrixWorkspace_sptr PhaseQuadMuon::squash(const API::MatrixWorkspace_sptr
238
247
}
239
248
std::vector<bool > emptySpectrum;
240
249
emptySpectrum.reserve (nspec);
241
- std::vector<CurveFitting::EigenVector > n0Vectors (nspec);
250
+ std::vector<Eigen::Vector<eigenDataType, 2 > > n0Vectors (nspec);
242
251
243
252
// Calculate coefficients aj, bj
244
253
@@ -254,23 +263,21 @@ API::MatrixWorkspace_sptr PhaseQuadMuon::squash(const API::MatrixWorkspace_sptr
254
263
const double phi = phase->Double (h, phaseIndex);
255
264
const double X = n0[h] * asym * cos (phi);
256
265
const double Y = n0[h] * asym * sin (phi);
257
- n0Vectors[h] = CurveFitting::EigenVector ({X, Y});
266
+ Eigen::Vector<eigenDataType, 2 > n0vec;
267
+ n0Vectors[h] = {X, Y};
258
268
sxx += X * X;
259
269
syy += Y * Y;
260
270
sxy += X * Y;
261
271
} else {
262
- n0Vectors[h] = CurveFitting::EigenVector ({ 0.0 , 0.0 } );
272
+ n0Vectors[h] = Eigen::Vector<eigenDataType, 2 >:: Zero ( );
263
273
}
264
274
}
265
275
266
- CurveFitting::EigenMatrix muLamMatrix (2 , 2 );
267
- muLamMatrix.set (0 , 0 , sxx);
268
- muLamMatrix.set (0 , 1 , sxy);
269
- muLamMatrix.set (1 , 0 , sxy);
270
- muLamMatrix.set (1 , 1 , syy);
271
- muLamMatrix.invert ();
276
+ Eigen::Matrix<eigenDataType, 2 , 2 > muLamMatrix;
277
+ muLamMatrix << sxx, sxy, sxy, syy;
278
+ muLamMatrix = Eigen::PartialPivLU<Eigen::Matrix<eigenDataType, 2 , 2 >>(muLamMatrix).inverse ();
272
279
273
- std::vector<double > aj (nspec), bj (nspec);
280
+ std::vector<eigenDataType > aj (nspec), bj (nspec);
274
281
for (size_t h = 0 ; h < nspec; h++) {
275
282
aj[h] = bj[h] = 0 ;
276
283
if (!emptySpectrum[h]) {
0 commit comments