Skip to content

Commit

Permalink
cleaned up LinearSolver and DgesvSolver. Also fixed minor errors.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OpenModelica-Hudson
qichenghua authored and OpenModelica-Hudson committed Jan 12, 2017
1 parent a558ce2 commit 40f96b8
Show file tree
Hide file tree
Showing 5 changed files with 59 additions and 148 deletions.
12 changes: 6 additions & 6 deletions SimulationRuntime/cpp/CMakeLists.txt
View file
Expand Up @@ -850,12 +850,6 @@ IF(USE_SUNDIALS)
GET_TARGET_PROPERTY(libKinsol ${KinsolName} LOCATION)
GET_FILENAME_COMPONENT(libKinsolName ${libKinsol} NAME)

GET_TARGET_PROPERTY(libLinearSolver ${LinearSolverName} LOCATION)
GET_FILENAME_COMPONENT(libLinearSolverName ${libLinearSolver} NAME)

GET_TARGET_PROPERTY(libDgesvSolver ${DgesvSolverName} LOCATION)
GET_FILENAME_COMPONENT(libDgesvSolverName ${libDgesvSolver} NAME)

#GET_TARGET_PROPERTY(libIdas ${IdasName} LOCATION)
#GET_FILENAME_COMPONENT(libIdasName ${libIdas} NAME)
#GET_TARGET_PROPERTY(libIda ${IdaName} LOCATION)
Expand All @@ -878,6 +872,12 @@ IF(USE_TRILINOS)
SET(NOX_LIB ${libNoxName})
ENDIF(USE_TRILINOS)

GET_TARGET_PROPERTY(libLinearSolver ${LinearSolverName} LOCATION)
GET_FILENAME_COMPONENT(libLinearSolverName ${libLinearSolver} NAME)

GET_TARGET_PROPERTY(libDgesvSolver ${DgesvSolverName} LOCATION)
GET_FILENAME_COMPONENT(libDgesvSolverName ${libDgesvSolver} NAME)

GET_TARGET_PROPERTY(libPeer ${PeerName} LOCATION)
GET_FILENAME_COMPONENT(libPeerName ${libPeer} NAME)

Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Core/Utils/Modelica/ModelicaSimulationError.h
View file
Expand Up @@ -8,7 +8,7 @@
//Enum for error types that can occur
enum SIMULATION_ERROR {
SOLVER, //all errors occur in solver (Euler,CVode)
ALGLOOP_SOLVER, //all errors occur in non-,lin-solver (Kinsol,Newton,Hybrj)
ALGLOOP_SOLVER, //all errors occur in non-,lin-solver (Nox,Kinsol,Newton,Hybrj)
MODEL_EQ_SYSTEM, //all errors occur in model system class during simulation
ALGLOOP_EQ_SYSTEM,//all errors occur in algloop system class during simulation
MODEL_FACTORY, //all errors occur model system factory classes
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/cpp/Include/Solver/LinearSolver/LinearSolver.h
View file
Expand Up @@ -57,8 +57,8 @@ class LinearSolver : public IAlgLoopSolver
*_zeroVec, ///zero vector
*_scale; //scaling parameter to prevent overflow in singular systems
bool _sparse;
bool _generateoutput; //prints nothing, if set to false. Prints Matrix, right hand side, and solution of the linear system, if set to true.
#if defined(klu)
int _dim;
klu_symbolic* _kluSymbolic ;
klu_numeric* _kluNumeric ;
klu_common* _kluCommon ;
Expand Down
9 changes: 5 additions & 4 deletions SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.cpp
View file
Expand Up @@ -106,11 +106,11 @@ void DgesvSolver::solve()
_iterationStatus = CONTINUE;

//use lapack
long int dimRHS = 1; // Dimension of right hand side of linear system (=_b)
long int irtrn = 0; // Return-flag of Fortran code
long int dimRHS = 1; // Dimension of right hand side of linear system (=_b)
long int irtrn = 0; // Return-flag of Fortran code

if(_algLoop->isLinearTearing())
_algLoop->setReal(_zeroVec); //if the system is linear Tearing it means that the system is of the form Ax-b=0, so plugging in x=0 yields -b for the left hand side
_algLoop->setReal(_zeroVec); //if the system is linear Tearing it means that the system is of the form Ax-b=0, so plugging in x=0 yields -b for the left hand side

_algLoop->evaluate();
_algLoop->getb(_b);
Expand All @@ -132,6 +132,7 @@ void DgesvSolver::solve()
_iterationStatus = DONE;


//we need to revert the sign of y, because the sign of b was changed before.
if(_algLoop->isLinearTearing()){
for(int i=0; i<_dimSys; i++)
_y[i]=-_b[i];
Expand All @@ -140,7 +141,7 @@ void DgesvSolver::solve()
}

_algLoop->setReal(_y);
if(_algLoop->isLinearTearing()) _algLoop->evaluate();//warum nur in diesem Fall??
if(_algLoop->isLinearTearing()) _algLoop->evaluate();//resets the right hand side to zero in the case of linear tearing. Otherwise, the b vector on the right hand side needs no update.
}

IAlgLoopSolver::ITERATIONSTATUS DgesvSolver::getIterationStatus()
Expand Down
182 changes: 46 additions & 136 deletions SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp
View file
Expand Up @@ -46,6 +46,7 @@ LinearSolver::LinearSolver(ILinearAlgLoop* algLoop, ILinSolverSettings* settings
, _iterationStatus (CONTINUE)
, _firstCall (true)
, _scale (NULL)
, _generateoutput (false)
{
_sparse = _algLoop->getUseSparseFormat();
}
Expand Down Expand Up @@ -187,38 +188,37 @@ void LinearSolver::solve()
_iterationStatus = CONTINUE;

//use lapack
long int dimRHS = 1; // Dimension of right hand side of linear system (=_b)
long int irtrn = 0; // Return-flag of Fortran code
long int dimRHS = 1; // Dimension of right hand side of linear system (=_b)
long int irtrn = 0; // Return-flag of Fortran code

if(_algLoop->isLinearTearing())
_algLoop->setReal(_zeroVec); //if the system is linear Tearing it means that the system is of the form Ax-b=0, so plugging in x=0 yields -b for the left hand side
_algLoop->setReal(_zeroVec); //if the system is linear tearing it means that the system is of the form Ax-b=0, so plugging in x=0 yields -b for the left hand side

_algLoop->evaluate();
_algLoop->getb(_b);

if (_sparse == false){
//if !_sparse, we use LAPACK routines, otherwise we use KLU to solve the linear system
if (!_sparse){
const matrix_t& A = _algLoop->getAMatrix();
const double* Atemp = A.data().begin();

memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double));

/*
//output routine
std::cout << std::endl;
std::cout << "We solve a linear system with coefficient matrix" << std::endl;
for (int i=0;i<_dimSys;i++){
for (int j=0;j<_dimSys;j++){
std::cout << Atemp[i+j*_dimSys] << " ";
if(_generateoutput){
std::cout << std::endl;
std::cout << "We solve a linear system with coefficient matrix" << std::endl;
for (int i=0;i<_dimSys;i++){
for (int j=0;j<_dimSys;j++){
std::cout << Atemp[i+j*_dimSys] << " ";
}
std::cout << std::endl;
}
std::cout << "and right hand side" << std::endl;
for (int i=0;i<_dimSys;i++){
std::cout << _b[i] << " ";
}
std::cout << std::endl;
}
std::cout << "and right hand side" << std::endl;
for (int i=0;i<_dimSys;i++){
std::cout << _b[i] << " ";
}
std::cout << std::endl;
*/


dgesv_(&_dimSys,&dimRHS,_A,&_dimSys,_ihelpArray,_b,&_dimSys,&irtrn);

Expand All @@ -232,122 +232,41 @@ void LinearSolver::solve()
else
_iterationStatus = DONE;
}else{

#if defined(klu)


//testing sparse with dense
/*this version is a test. it extracts the dense format out of the sparse format and uses the dense lapack solver to sove the dense problem.
//writing entries of A
sparsematrix_t& A = _algLoop->getSparseAMatrix();
_Ax= boost::numeric::bindings::begin_value (A);

if (_generateoutput){
double* a = new double[_dimSys*_dimSys];
memset(a, 0, _dimSys*_dimSys*sizeof(double));

double** a = new double*[_dimSys];
double* asdf = new double[_dimSys*_dimSys];
for(int i=0;i<_dimSys;i++){
a[i]=new double[_dimSys];
}
for(int i=0;i<_dimSys;i++){
for(int j=0;j<_dimSys;j++){
a[i][j]=0;
for(int k=_Ap[j];k<_Ap[j+1];k++)
if(i==_Ai[k])
a[i][j]=_Ax[k];
}
}
for(int i=0;i<_dimSys;i++){
for(int j=0;j<_dimSys;j++){
asdf[i+j*_dimSys]=a[i][j];
for(int i=0;i<_dimSys;i++){
for(int j=0;j<_dimSys;j++){
for(int k=_Ap[j];k<_Ap[j+1];k++)
if(i==_Ai[k])
a[i+j*_dimSys]=_Ax[k];
}
}
}

//output routine
std::cout << std::endl;
std::cout << "We solve a linear system with coefficient matrix" << std::endl;
for (int i=0;i<_dimSys;i++){
for (int j=0;j<_dimSys;j++){
std::cout << asdf[i+j*_dimSys] << " ";
}
std::cout << std::endl;
}
std::cout << "and right hand side" << std::endl;
for (int i=0;i<_dimSys;i++){
std::cout << _b[i] << " ";
}
std::cout << std::endl;
dgesv_(&_dimSys,&dimRHS,asdf,&_dimSys,_ihelpArray,_b,&_dimSys,&irtrn);
_iterationStatus = DONE;
for(int i=0;i<_dimSys;i++){
delete [] a[i];
}
delete [] a;
delete [] asdf;
/*test for ccs-dense-conversion
// testSparse2.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
int _tmain(int argc, _TCHAR* argv[])
{
int _dimSys = 3;
int _Ap[4] = {0,3,5,6};
int _Ai[6] = { 0, 1, 2, 1, 2, 2 };
int _Ax[6] = { 5, 4, 3, 2, 1, 8 };
double** a = new double*[_dimSys];
for (int i = 0; i<_dimSys; i++){
a[i] = new double[_dimSys];
}
for (int i = 0; i<_dimSys; i++){
for (int j = 0; j<_dimSys; j++){
a[i][j] = 0;
for (int k = _Ap[j]; k<_Ap[j + 1]; k++)
if (i == _Ai[k])
a[i][j] += _Ax[k];
std::cout << "We solve a linear system with coefficient matrix" << std::endl;
for (int i=0;i<_dimSys;i++){
for (int j=0;j<_dimSys;j++){
std::cout << a[i+j*_dimSys] << " ";
}
std::cout << std::endl;
}

for (int i = 0; i < _dimSys; i++){
for (int j = 0; j < _dimSys; j++){
printf("a(%i,%i)=%e\n", i, j, a[i][j]);
}
}
delete [] a;


for (int i = 0; i<_dimSys; i++){
delete[] a[i];
std::cout << "and right hand side" << std::endl;
for (int i=0;i<_dimSys;i++){
std::cout << _b[i] << " ";
}
delete[] a;
return 0;
std::cout << std::endl;
}
*/


//version with klu

/*
//testing, whether Ax is modified
double *Ax=NULL;
Ax = new double[_nonzeros];
for(int i=0;i<_nonzeros;i++) Ax[i]=_Ax[i];
*/


//writing entries of A
sparsematrix_t& A = _algLoop->getSparseAMatrix();
_Ax= boost::numeric::bindings::begin_value (A);

int ok = klu_refactor (_Ap, _Ai, _Ax, _kluSymbolic, _kluNumeric, _kluCommon) ;

Expand All @@ -363,38 +282,29 @@ void LinearSolver::solve()
ok=klu_solve (_kluSymbolic, _kluNumeric, _dimSys, 1, _b, _kluCommon) ;
if (ok!=1) throw ModelicaSimulationError(ALGLOOP_SOLVER,"error solving Sparse Solver KLU");

/*
//testing, whether Ax is modified
for(int i=0;i<_nonzeros;i++) if(Ax[i]!=_Ax[i]) std::cout << "Ax was modified" << std::endl;
if(Ax)
delete [] Ax;
*/

#else
throw ModelicaSimulationError(ALGLOOP_SOLVER,"error solving linear system with klu not implemented");
#endif
}

//we need to revert the sign of y, because the sign of b was changed before.
if(_algLoop->isLinearTearing()){
for(int i=0; i<_dimSys; i++)
_y[i]=-_b[i];
}else{
memcpy(_y,_b,_dimSys*sizeof(double));
}


/*
//output routine
std::cout << "The solution of the linear system is given by" << std::endl;
for (int i=0;i<_dimSys;i++){
std::cout << _y[i] << " ";
if (_generateoutput){
std::cout << "The solution of the linear system is given by" << std::endl;
for (int i=0;i<_dimSys;i++){
std::cout << _y[i] << " ";
}
std::cout << std::endl;
}
std::cout << std::endl;
*/


_algLoop->setReal(_y);
if(_algLoop->isLinearTearing()) _algLoop->evaluate();//warum nur in diesem Fall??
if(_algLoop->isLinearTearing()) _algLoop->evaluate();//resets the right hand side to zero in the case of linear tearing. Otherwise, the b vector on the right hand side needs no update.
}

IAlgLoopSolver::ITERATIONSTATUS LinearSolver::getIterationStatus()
Expand Down

0 comments on commit 40f96b8

Please sign in to comment.