changed A Matrix to ublas matrix in Algloop class

niklwors · niklwors · commit af9687a1c39d · 2015-07-24T08:33:50.000+02:00
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
diff --git a/Compiler/Template/CodegenCpp.tpl b/Compiler/Template/CodegenCpp.tpl
diff --git a/SimulationRuntime/cpp/Core/Math/SparseMatrix.cpp b/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
diff --git a/SimulationRuntime/cpp/Core/Modelica/ModelicaSystem.cpp b/SimulationRuntime/cpp/Core/Modelica/ModelicaSystem.cpp
@@ -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");
 }
diff --git a/SimulationRuntime/cpp/Core/Solver/SystemStateSelection.cpp b/SimulationRuntime/cpp/Core/Solver/SystemStateSelection.cpp
@@ -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 */
diff --git a/SimulationRuntime/cpp/Core/System/SimVars.cpp b/SimulationRuntime/cpp/Core/System/SimVars.cpp
@@ -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[])
diff --git a/SimulationRuntime/cpp/Include/Core/Modelica.h b/SimulationRuntime/cpp/Include/Core/Modelica.h
@@ -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>
diff --git a/SimulationRuntime/cpp/Include/Core/Modelica/ModelicaSystem.h b/SimulationRuntime/cpp/Include/Core/Modelica/ModelicaSystem.h
@@ -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() ;
 
diff --git a/SimulationRuntime/cpp/Include/Core/System/DiscreteEvents.h b/SimulationRuntime/cpp/Include/Core/System/DiscreteEvents.h
@@ -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;
diff --git a/SimulationRuntime/cpp/Include/Core/System/IAlgLoop.h b/SimulationRuntime/cpp/Include/Core/System/IAlgLoop.h
@@ -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;
diff --git a/SimulationRuntime/cpp/Include/Core/System/IMixedSystem.h b/SimulationRuntime/cpp/Include/Core/System/IMixedSystem.h
@@ -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;
 
diff --git a/SimulationRuntime/cpp/Include/Core/Utils/numeric/utils.h b/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
diff --git a/SimulationRuntime/cpp/Include/Solver/CVode/CVode.h b/SimulationRuntime/cpp/Include/Solver/CVode/CVode.h
@@ -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;
diff --git a/SimulationRuntime/cpp/Include/Solver/Hybrj/Hybrj.h b/SimulationRuntime/cpp/Include/Solver/Hybrj/Hybrj.h
@@ -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
@@ -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;
diff --git a/SimulationRuntime/cpp/Include/Solver/IDA/IDA.h b/SimulationRuntime/cpp/Include/Solver/IDA/IDA.h
@@ -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;
diff --git a/SimulationRuntime/cpp/Include/Solver/Kinsol/Kinsol.h b/SimulationRuntime/cpp/Include/Solver/Kinsol/Kinsol.h
@@ -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
@@ -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);
 
@@ -95,7 +97,7 @@ class Kinsol : public IAlgLoopSolver
 
   realtype _fnorm,
     _currentIterateNorm;
-
+   ublas::permutation_matrix<double>* _P;
    int _counter;
 };
 /** @} */ // end of solverKinsol
diff --git a/SimulationRuntime/cpp/Include/Solver/UmfPack/UmfPack.h b/SimulationRuntime/cpp/Include/Solver/UmfPack/UmfPack.h
@@ -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;
diff --git a/SimulationRuntime/cpp/Solver/CVode/CVode.cpp b/SimulationRuntime/cpp/Solver/CVode/CVode.cpp
diff --git a/SimulationRuntime/cpp/Solver/Hybrj/Hybrj.cpp b/SimulationRuntime/cpp/Solver/Hybrj/Hybrj.cpp
diff --git a/SimulationRuntime/cpp/Solver/IDA/IDA.cpp b/SimulationRuntime/cpp/Solver/IDA/IDA.cpp
diff --git a/SimulationRuntime/cpp/Solver/Kinsol/Kinsol.cpp b/SimulationRuntime/cpp/Solver/Kinsol/Kinsol.cpp
diff --git a/SimulationRuntime/cpp/Solver/Newton/Newton.cpp b/SimulationRuntime/cpp/Solver/Newton/Newton.cpp
diff --git a/SimulationRuntime/cpp/Solver/UmfPack/UmfPack.cpp b/SimulationRuntime/cpp/Solver/UmfPack/UmfPack.cpp

Original file line number	Diff line number	Diff line change
`@@ -181,12 +181,12 @@ bool Modelica::checkConditions()`
`181`	`181`	`throw ModelicaSimulationError(MODEL_EQ_SYSTEM,"checkConditions is not yet implemented");`
`182`	`182`	`}`
`183`	`183`
`184`		`-void Modelica::getJacobian(SparseMatrix& matrix)`
	`184`	`+void Modelica::getJacobian(matrix_t& matrix)`
`185`	`185`	`{`
`186`	`186`	`throw ModelicaSimulationError(MODEL_EQ_SYSTEM,"giveJacobian is not yet implemented");`
`187`	`187`	`}`
`188`	`188`
`189`		`-void Modelica::getStateSetJacobian(SparseMatrix& matrix)`
	`189`	`+void Modelica::getStateSetJacobian(matrix_t& matrix)`
`190`	`190`	`{`
`191`	`191`	`throw ModelicaSimulationError(MODEL_EQ_SYSTEM,"giveStateJacobian is not yet implemented");`
`192`	`192`	`}`
Original file line number	Diff line number	Diff line change
`@@ -76,7 +76,7 @@ return true;`
`76`	`76`	`{`
`77`	`77`	`boost::shared_array<int> oldColPivot(new int[_dimStateCanditates[i]]);`
`78`	`78`	`boost::shared_array<int> oldRowPivot(new int[_dimDummyStates[i]]);`
`79`		`- SparseMatrix stateset_matrix;`
	`79`	`+ matrix_t stateset_matrix;`
`80`	`80`	`_system->getStateSetJacobian(i,stateset_matrix);`
`81`	`81`
`82`	`82`	`/* call pivoting function to select the states */`
Original file line number	Diff line number	Diff line change
`@@ -235,12 +235,7 @@ bool* SimVars::initBoolArrayVar(size_t size, size_t start_index)`
`235`	`235`	`*/`
`236`	`236`	`void SimVars::initRealAliasArray(int indices[], size_t n, double* ref_data[])`
`237`	`237`	`{`
`238`		`- for(int i = 0; i < n; i++)`
`239`		`- {`
`240`		`- int index = indices[i];`
`241`		`- double* refToVar = SimVars::getRealVar(index);`
`242`		`- ref_data[i] = refToVar;`
`243`		`- }`
	`238`	`+ std::transform(indices,indices+n,ref_data,boost::lambda::bind(&SimVars::getRealVar,this,boost::lambda::_1));`
`244`	`239`	`}`
`245`	`240`
`246`	`241`	`void SimVars::initRealAliasArray(std::vector<int> indices, double* ref_data[])`