diff --git a/MathLib/ODE/CVodeSolver.cpp b/MathLib/ODE/CVodeSolver.cpp
index c4cbde8da151ae751ac86da0b6e99e3ecd652ca7..6b9cb14f2971ec01525980755f4ff825790dea64 100644
--- a/MathLib/ODE/CVodeSolver.cpp
+++ b/MathLib/ODE/CVodeSolver.cpp
@@ -257,7 +257,7 @@ void CVodeSolverImpl::preSolve()
static_cast<FunctionHandles*>(function_handles)
->callJacobian(t, NV_DATA_S(y), NV_DATA_S(ydot),
DENSE_COL(jac, 0),
- BaseLib::StorageOrder::ColumnMajor);
+ MathLib::StorageOrder::ColumnMajor);
return successful ? 0 : 1;
};
@@ -308,7 +308,7 @@ CVodeSolverImpl::~CVodeSolverImpl()
}
CVodeSolverInternal::CVodeSolverInternal(BaseLib::ConfigTree const& config)
- : _impl.reset(std::unique_ptr<CVodeSolverImpl>{new CVodeSolverImpl{config}})
+ : _impl{new CVodeSolverImpl{config}}
{
}
diff --git a/MathLib/ODE/OdeSolver.h b/MathLib/ODE/OdeSolver.h
index 184cd1cf0bbabc202a02b2cfcb17df3be946ef29..227d32d14f0b0c9eece1366599bb8d461c009a71 100644
--- a/MathLib/ODE/OdeSolver.h
+++ b/MathLib/ODE/OdeSolver.h
@@ -28,7 +28,8 @@ class OdeSolver
{
public:
using Arr = std::array<double, NumEquations>;
- using ConstArrRef = BaseLib::ArrayRef<const double, NumEquations>;
+ using ConstArrRef =
+ Eigen::Map<const Eigen::Matrix<double, NumEquations, 1>>;
using Function = MathLib::Function<NumEquations, FunctionArguments...>;
using JacobianFunction =
MathLib::JacobianFunction<NumEquations, FunctionArguments...>;
diff --git a/MathLib/ODE/OdeSolverFactory.h b/MathLib/ODE/OdeSolverFactory.h
index a6ec30b68f21e0e20dd7c87d5a293b3b94811c0c..5b9d7c7beec3f179125914c93c3e2979e3f18559 100644
--- a/MathLib/ODE/OdeSolverFactory.h
+++ b/MathLib/ODE/OdeSolverFactory.h
@@ -40,21 +40,21 @@ struct Handles<N, FunctionArgument> : public MathLib::FunctionHandles
// alike
if (f)
return f(t,
- BaseLib::ArrayRef<const double, N>{y},
- BaseLib::ArrayRef<double, N>{ydot},
+ Eigen::Map<const Eigen::Matrix<double, N, 1>>{y},
+ Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot},
*_data);
return true;
}
bool callJacobian(const double t, const double* const y,
const double* const ydot, double* const jac,
- BaseLib::StorageOrder order) override
+ MathLib::StorageOrder order) override
{
if (df)
return df(t,
- BaseLib::ArrayRef<const double, N>{y},
- BaseLib::ArrayRef<const double, N>{ydot},
- BaseLib::MatrixRef<double, N, N>{jac, order},
+ Eigen::Map<const Eigen::Matrix<double, N, 1>>{y},
+ Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot},
+ Eigen::Map<Eigen::Matrix<double, N, N>>{jac /*, order*/},
*_data);
return true;
}
@@ -85,21 +85,25 @@ struct Handles<N> : public MathLib::FunctionHandles
double* const ydot) override
{
if (f)
+ {
+ // auto ydot_ = Eigen::Map<Eigen::Matrix<double, N, 1>>{y};
+ // auto ydot_ = Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot};
return f(t,
- BaseLib::ArrayRef<const double, N>{y},
- BaseLib::ArrayRef<double, N>{ydot});
+ Eigen::Map<const Eigen::Matrix<double, N, 1>>{y},
+ Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot});
+ }
return true;
}
bool callJacobian(const double t, const double* const y,
const double* const ydot, double* const jac,
- BaseLib::StorageOrder order) override
+ MathLib::StorageOrder order) override
{
if (df)
return df(t,
- BaseLib::ArrayRef<const double, N>{y},
- BaseLib::ArrayRef<const double, N>{ydot},
- BaseLib::MatrixRef<double, N, N>{jac, order});
+ Eigen::Map<const Eigen::Matrix<double, N, 1>>{y},
+ Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot},
+ Eigen::Map<Eigen::Matrix<double, N, N>>{jac /*, order*/});
return true;
}
diff --git a/MathLib/ODE/declarations.h b/MathLib/ODE/declarations.h
index 473e883d414866d32a79106327dc2d47e7f271b1..000f4246a4a6e10350ff82aaa259c47090252d2d 100644
--- a/MathLib/ODE/declarations.h
+++ b/MathLib/ODE/declarations.h
@@ -10,25 +10,29 @@
#ifndef MATHLIB_ODE_DECLARATIONS_H
#define MATHLIB_ODE_DECLARATIONS_H
-#include <cassert>
-
-#include "BaseLib/ArrayRef.h"
-#include "BaseLib/MatrixRef.h"
+#include <Eigen/Core>
namespace MathLib
{
+enum class StorageOrder
+{
+ ColumnMajor,
+ RowMajor
+};
+
template <unsigned N, typename... FunctionArguments>
using Function = bool (*)(const double t,
- BaseLib::ArrayRef<const double, N> y,
- BaseLib::ArrayRef<double, N> ydot,
+ Eigen::Map<const Eigen::Matrix<double, N, 1>> const y,
+ Eigen::Map<Eigen::Matrix<double, N, 1>> ydot,
FunctionArguments&... arg);
template <unsigned N, typename... FunctionArguments>
-using JacobianFunction = bool (*)(const double t,
- BaseLib::ArrayRef<const double, N> y,
- BaseLib::ArrayRef<const double, N> ydot,
- BaseLib::MatrixRef<double, N, N> jac,
- FunctionArguments&... arg);
+using JacobianFunction =
+ bool (*)(const double t,
+ Eigen::Map<const Eigen::Matrix<double, N, 1>> const y,
+ Eigen::Map<Eigen::Matrix<double, N, 1>> ydot,
+ Eigen::Map<Eigen::Matrix<double, N, N>> jac,
+ FunctionArguments&... arg);
// This is an internal detail
class FunctionHandles
@@ -38,9 +42,9 @@ public:
double* const ydot) = 0;
virtual bool callJacobian(const double t,
double const* const y,
- double const* const ydot,
+ double* const ydot,
double* const jac,
- BaseLib::StorageOrder order) = 0;
+ StorageOrder order) = 0;
virtual bool hasJacobian() const = 0;
diff --git a/Tests/MathLib/TestCVode.cpp b/Tests/MathLib/TestCVode.cpp
index 7d601b24fb6463a6ccafdbf2db8521c6a61194c3..8b2ed6ea589f0e056f42e6f71e71494bda07dd09 100644
--- a/Tests/MathLib/TestCVode.cpp
+++ b/Tests/MathLib/TestCVode.cpp
@@ -16,8 +16,8 @@
#include <cstdio>
bool f(const double,
- BaseLib::ArrayRef<const double, 1> y,
- BaseLib::ArrayRef<double, 1> ydot)
+ Eigen::Map<const Eigen::Matrix<double, 1, 1>> const y,
+ Eigen::Map<Eigen::Matrix<double, 1, 1>> ydot)
{
if (y[0] <= 0.0) return false;
@@ -26,9 +26,9 @@ bool f(const double,
}
bool df(const double /*t*/,
- BaseLib::ArrayRef<const double, 1> y,
- BaseLib::ArrayRef<const double, 1> /*ydot*/,
- BaseLib::MatrixRef<double, 1, 1> jac)
+ Eigen::Map<const Eigen::Matrix<double, 1, 1>> const y,
+ Eigen::Map<Eigen::Matrix<double, 1, 1>> /*ydot*/,
+ Eigen::Map<Eigen::Matrix<double, 1, 1>> jac)
{
if (y[0] <= 0.0) return false;
@@ -42,8 +42,8 @@ struct ExtraData
};
bool f_extra(const double,
- BaseLib::ArrayRef<const double, 1> y,
- BaseLib::ArrayRef<double, 1> ydot,
+ Eigen::Map<const Eigen::Matrix<double, 1, 1>> const y,
+ Eigen::Map<Eigen::Matrix<double, 1, 1>> ydot,
ExtraData& data)
{
if (y[0] <= 0.0) return false;