Skip to content

Commit

Permalink
changed A Matrix to ublas matrix in Algloop class
Browse files Browse the repository at this point in the history
removed template functions in AlgLoop classes to calc A matrix, only one A matrix type can be used now(sparse matrix can be used with a flag later)
adapted non linear solvers on new getSystemMatrix method
  • Loading branch information
niklwors authored and OpenModelica-Hudson committed Jul 24, 2015
1 parent 12822bc commit 5ff3cc0
Show file tree
Hide file tree
Showing 22 changed files with 1,606 additions and 1,247 deletions.
168 changes: 121 additions & 47 deletions Compiler/Template/CodegenCpp.tpl

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Core/Math/SparseMatrix.cpp
@@ -1,6 +1,6 @@
#include <Core/ModelicaDefine.h>
#include <Core/Modelica.h>
#include <Core/Math/SparseMatrix.h>
#include <Core/Math/matrix_t.h>
#ifdef USE_UMFPACK
#include "umfpack.h"
#endif
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/cpp/Core/Modelica/ModelicaSystem.cpp
Expand Up @@ -181,12 +181,12 @@ bool Modelica::checkConditions()
throw ModelicaSimulationError(MODEL_EQ_SYSTEM,"checkConditions is not yet implemented");
}

void Modelica::getJacobian(SparseMatrix& matrix)
void Modelica::getJacobian(matrix_t& matrix)
{
throw ModelicaSimulationError(MODEL_EQ_SYSTEM,"giveJacobian is not yet implemented");
}

void Modelica::getStateSetJacobian(SparseMatrix& matrix)
void Modelica::getStateSetJacobian(matrix_t& matrix)
{
throw ModelicaSimulationError(MODEL_EQ_SYSTEM,"giveStateJacobian is not yet implemented");
}
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Core/Solver/SystemStateSelection.cpp
Expand Up @@ -76,7 +76,7 @@ return true;
{
boost::shared_array<int> oldColPivot(new int[_dimStateCanditates[i]]);
boost::shared_array<int> oldRowPivot(new int[_dimDummyStates[i]]);
SparseMatrix stateset_matrix;
matrix_t stateset_matrix;
_system->getStateSetJacobian(i,stateset_matrix);

/* call pivoting function to select the states */
Expand Down
7 changes: 1 addition & 6 deletions SimulationRuntime/cpp/Core/System/SimVars.cpp
Expand Up @@ -235,12 +235,7 @@ bool* SimVars::initBoolArrayVar(size_t size, size_t start_index)
*/
void SimVars::initRealAliasArray(int indices[], size_t n, double* ref_data[])
{
for(int i = 0; i < n; i++)
{
int index = indices[i];
double* refToVar = SimVars::getRealVar(index);
ref_data[i] = refToVar;
}
std::transform(indices,indices+n,ref_data,boost::lambda::bind(&SimVars::getRealVar,this,boost::lambda::_1));
}

void SimVars::initRealAliasArray(std::vector<int> indices, double* ref_data[])
Expand Down
9 changes: 5 additions & 4 deletions SimulationRuntime/cpp/Include/Core/Modelica.h
Expand Up @@ -89,10 +89,11 @@ using std::vector;
//using boost::timer::nanosecond_type;
typedef ublas::shallow_array_adaptor<double> adaptor_t;
typedef ublas::vector<double, adaptor_t> shared_vector_t;
typedef ublas::matrix<double, adaptor_t> shared_matrix_t;
typedef boost::function<bool (unsigned int)> getCondition_type;
typedef boost::function<void (unordered_map<string,unsigned int>&,unordered_map<string,unsigned int>&)> init_prevars_type;
typedef uBlas::compressed_matrix<double, uBlas::row_major, 0, uBlas::unbounded_array<int>, uBlas::unbounded_array<double> > SparseMatrix;
typedef ublas::matrix<double, ublas::column_major,adaptor_t> shared_matrix_t;
//typedef boost::function<bool (unsigned int)> getCondition_type;
//typedef boost::function<void (unordered_map<string,unsigned int>&,unordered_map<string,unsigned int>&)> init_prevars_type;
typedef uBlas::compressed_matrix<double, uBlas::column_major, 0, uBlas::unbounded_array<int>, uBlas::unbounded_array<double> > sparsematrix_t;
typedef ublas::matrix<double, ublas::column_major> matrix_t;
#include <Core/SimulationSettings/ISettingsFactory.h>
#include <SimCoreFactory/Policies/FactoryConfig.h>
#include <Core/Utils/Modelica/ModelicaSimulationError.h>
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/cpp/Include/Core/Modelica/ModelicaSystem.h
Expand Up @@ -83,8 +83,8 @@ class Modelica : public IMixedSystem, public IContinuous, public IEvent, public
virtual void setTime(const double& t);

// Provide Jacobian
virtual void getJacobian(SparseMatrix& matrix);
virtual void getStateSetJacobian(SparseMatrix& matrix);
virtual void getJacobian(matrix_t& matrix);
virtual void getStateSetJacobian(matrix_t& matrix);
// Provide number (dimension) of zero functions
virtual int getDimZeroFunc() ;

Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Core/System/DiscreteEvents.h
Expand Up @@ -41,7 +41,7 @@ class BOOST_EXTENSION_EVENTHANDLING_DECL DiscreteEvents
bool changeDiscreteVar(double& var);
bool changeDiscreteVar(int& var);
bool changeDiscreteVar(bool& var);
getCondition_type getCondition;
//getCondition_type getCondition;

private:
boost::shared_ptr<ISimVars> _sim_vars;
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/cpp/Include/Core/System/IAlgLoop.h
Expand Up @@ -54,8 +54,8 @@ class IAlgLoop
/// Provide the right hand side (according to the index)
virtual void getRHS(double* res) = 0;

virtual void getSystemMatrix(double* A_matrix) = 0;
virtual void getSystemMatrix(SparseMatrix&) = 0;

virtual const matrix_t& getSystemMatrix() = 0;

virtual bool isLinear() = 0;
virtual bool isLinearTearing() = 0;
Expand Down
6 changes: 3 additions & 3 deletions SimulationRuntime/cpp/Include/Core/System/IMixedSystem.h
Expand Up @@ -26,9 +26,9 @@ class IMixedSystem
public:
virtual ~IMixedSystem() {};
/// Provide Jacobian
virtual void getJacobian(SparseMatrix& matrix) = 0;
virtual void getJacobian(SparseMatrix& matrix,unsigned int index) = 0;
virtual void getStateSetJacobian(unsigned int index, SparseMatrix& matrix) = 0;
virtual const matrix_t& getJacobian() = 0;
virtual const matrix_t& getJacobian(unsigned int index) = 0;
virtual void getStateSetJacobian(unsigned int index, matrix_t& matrix) = 0;
/// Called to handle all events occured at same time
virtual bool handleSystemEvents(bool* events) = 0;

Expand Down
204 changes: 204 additions & 0 deletions SimulationRuntime/cpp/Include/Core/Utils/numeric/utils.h
@@ -0,0 +1,204 @@

#ifndef F_UTILS_H
#define F_UTILS_H

#include <stddef.h>
#include <iostream>
#include <Core/Utils/numeric/bindings/size.hpp>
#include <Core/Utils/numeric/bindings/begin.hpp>
#include <Core/Utils/numeric/bindings/value_type.hpp>


///////////////////////////////
// vectors

// element access:
template <typename V>
struct vct_access_traits {
typedef typename
boost::numeric::bindings::value_type<V>::type val_t;
typedef val_t& ref_t;
static ref_t elem (V& v, size_t i) { return v[i]; }
};

template <typename V>
struct vct_access_traits<V const> {
typedef typename
boost::numeric::bindings::value_type<V>::type val_t;
typedef val_t ref_t;
static ref_t elem (V const& v, size_t i) { return v[i]; }
};

template <typename V>
inline
typename vct_access_traits<V>::ref_t elem_v (V& v, size_t i) {
return vct_access_traits<V>::elem (v, i);
}

// initialization:
struct ident {
size_t operator() (size_t i) const { return i; }
};
struct iplus1 {
size_t operator() (size_t i) const { return i + 1; }
};
template <typename T>
struct const_val {
T val;
const_val (T v = T()) : val (v) {}
T operator() (size_t) const { return val; }
T operator() (size_t, size_t) const { return val; }
};
struct kpp {
size_t val;
kpp (size_t v = 0) : val (v) {}
size_t operator() (size_t) {
size_t tmp = val;
++val;
return tmp;
}
size_t operator() (size_t, size_t) {
size_t tmp = val;
++val;
return tmp;
}
};
template <typename T>
struct times_plus {
T s, a, b;
times_plus (T ss, T aa = T(), T bb = T()) : s (ss), a (aa), b (bb) {}
T operator() (size_t i) const {
return s * (T) i + a;
}
T operator() (size_t i, size_t j) const {
return s * ((T) i + a) + (T) j + b;
}
};

template <typename F, typename V>
void init_v (V& v, F f = F()) {
size_t sz
= boost::numeric::bindings::size( v );
for (std::size_t i = 0; i < sz; ++i)
elem_v (v, i) = f (i);
}

// printing:
template <typename V>
void print_v (V const& v, char const* ch = 0) {
if (ch)
std::cout << ch << ": ";
size_t sz
= boost::numeric::bindings::size( v );
for (std::size_t i = 0; i < sz; ++i)
std::cout << elem_v (v, i) << " ";
std::cout << std::endl;
}


/////////////////////////////////////
// matrices

// element access:
template <typename M>
struct matr_access_traits {
typedef M mat_t;
typedef typename boost::numeric::bindings::value_type<mat_t>::type val_t;
typedef typename mat_t::reference ref_t;
static ref_t elem (mat_t& m, size_t i, size_t j) { return m (i, j); }
};

#ifdef EIGEN_FORWARDDECLARATIONS_H
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct matr_access_traits<Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> > {
typedef Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> mat_t;
typedef typename boost::numeric::bindings::value_type<mat_t>::type val_t;
typedef val_t& ref_t;
static ref_t elem (mat_t& m, size_t i, size_t j) { return m (i, j); }
};
#endif

template <typename M>
struct matr_access_traits<M const> {
typedef M const mat_t;
typedef typename boost::numeric::bindings::value_type<mat_t>::type val_t;
typedef val_t ref_t;
static ref_t elem (mat_t& m, size_t i, size_t j) { return m (i, j); }
};

template <typename M>
inline
typename matr_access_traits<M>::ref_t elem_m (M& m, size_t i, size_t j) {
return matr_access_traits<M>::elem (m, i, j);
}

// initialization:
struct rws {
size_t operator() (size_t i, size_t) const { return i; }
};
struct rws1 {
size_t operator() (size_t i, size_t) const { return i + 1; }
};
struct cls {
size_t operator() (size_t, size_t j) const { return j; }
};
struct cls1 {
size_t operator() (size_t, size_t j) const { return j + 1; }
};

template <typename F, typename M>
void init_m (M& m, F f = F()) {
size_t sz1
= boost::numeric::bindings::size_row( m );
size_t sz2
= boost::numeric::bindings::size_column( m );
for (std::size_t i = 0; i < sz1; ++i)
for (std::size_t j = 0; j < sz2; ++j)
elem_m (m, i, j) = f (i, j);
}

template <typename M>
void init_symm (M& m, char uplo = 'f') {
size_t n
= boost::numeric::bindings::size_row( m );
for (size_t i = 0; i < n; ++i) {
elem_m (m, i, i) = n;
for (size_t j = i + 1; j < n; ++j) {
if (uplo == 'u' || uplo == 'U')
elem_m (m, i, j) = n - (j - i);
else if (uplo == 'l' || uplo == 'L')
elem_m (m, j, i) = n - (j - i);
else
elem_m (m, i, j) = elem_m (m, j, i) = n - (j - i);
}
}
}

// printing:
template <typename M>
void print_m (M const& m, char const* ch = 0) {
if (ch)
std::cout << ch << ":\n";
size_t sz1
= boost::numeric::bindings::size_row( m );
size_t sz2
= boost::numeric::bindings::size_column( m );
for (std::size_t i = 0 ; i < sz1 ; ++i) {
for (std::size_t j = 0 ; j < sz2 ; ++j)
std::cout << elem_m (m, i, j) << " ";
std::cout << std::endl;
}
// std::cout << std::endl;
}

template <typename M>
void print_m_data (M const& m, char const* ch = 0) {
if (ch)
std::cout << ch << " data:\n";
using namespace boost::numeric::bindings;
std::copy( begin_value( m ), end_value( m ), std::ostream_iterator
< typename value_type< const M >::type >( std::cout, " " ) );
std::cout << std::endl;
}

#endif
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Solver/CVode/CVode.h
Expand Up @@ -167,7 +167,7 @@ double
// Variables for Coloured Jacobians
int* _colorOfColumn;
int _maxColors;
SparseMatrix _jacobianA;
matrix_t _jacobianA;
int _jacobianANonzeros;
int const* _jacobianAIndex;
int const* _jacobianALeadindex;
Expand Down
5 changes: 3 additions & 2 deletions SimulationRuntime/cpp/Include/Solver/Hybrj/Hybrj.h
Expand Up @@ -55,7 +55,7 @@ class Hybrj : public IAlgLoopSolver
ITERATIONSTATUS
_iterationStatus; ///< Output - Denotes the status of iteration

int
int
_dimSys; ///< Temp - Number of unknowns (=dimension of system of equations)

bool
Expand All @@ -75,7 +75,8 @@ class Hybrj : public IAlgLoopSolver
*_x_nom, //Nominal value of unknown variables
*_x_ex, //extraplated unknown varibales
*_x_scale, //current scale factor of unknown varibales
_t0, //old time
*_zeroVec,
_t0, //old time
_t1, //old time
_t2; //old time
bool _usescale;
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Solver/IDA/IDA.h
Expand Up @@ -162,7 +162,7 @@ double
// Variables for Coloured Jacobians
int* _colorOfColumn;
int _maxColors;
SparseMatrix _jacobianA;
matrix_t _jacobianA;
int _jacobianANonzeros;
int const* _jacobianAIndex;
int const* _jacobianALeadindex;
Expand Down
8 changes: 5 additions & 3 deletions SimulationRuntime/cpp/Include/Solver/Kinsol/Kinsol.h
Expand Up @@ -15,6 +15,9 @@
#include <kinsol/kinsol_spbcgs.h>
#include <kinsol/kinsol_sptfqmr.h>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>

//#include<kinsol_lapack.h>
int kin_fCallback(N_Vector y, N_Vector fval, void *user_data);
class Kinsol : public IAlgLoopSolver
Expand All @@ -36,8 +39,7 @@ class Kinsol : public IAlgLoopSolver
private:
/// Encapsulation of determination of residuals to given unknowns
void calcFunction(const double* y, double* residual);
/// Encapsulation of determination of Jacobian
void calcJacobian(double* f, double* y);


int check_flag(void *flagvalue, char *funcname, int opt);

Expand Down Expand Up @@ -95,7 +97,7 @@ class Kinsol : public IAlgLoopSolver

realtype _fnorm,
_currentIterateNorm;

ublas::permutation_matrix<double>* _P;
int _counter;
};
/** @} */ // end of solverKinsol
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Solver/UmfPack/UmfPack.h
Expand Up @@ -31,7 +31,7 @@ class UmfPack : public IAlgLoopSolver
ITERATIONSTATUS _iterationStatus;
ILinSolverSettings *_umfpackSettings;
IAlgLoop *_algLoop;
boost::shared_ptr<SparseMatrix> _jacs;
boost::shared_ptr<matrix_t> _jacs;
double * _jacd;
double * _rhs;
double * _x;
Expand Down

0 comments on commit 5ff3cc0

Please sign in to comment.