fix ticket:4213, added scaling of linear systems back again

Author: qichenghua

SimulationRuntime/cpp/Include/Solver/Dgesv/DgesvSolver.h

@@ -49,6 +49,7 @@ class DgesvSolver : public IAlgLoopSolver
       *_y_new,				//stores new solution
 	  *_b,                  ///< right hand side
 	  *_A,				///coefficients of linear system
-	  *_zeroVec;			///zero vector
+	  *_zeroVec,			///zero vector
+	  *_fNominal;
 };
 /** @} */ // end of solverLinearSolver

SimulationRuntime/cpp/Include/Solver/LinearSolver/LinearSolver.h

@@ -66,6 +66,8 @@ class LinearSolver : public IAlgLoopSolver
   int* _Ap;
   double* _Ax;
   int _nonzeros;
+#else
+	double *_fNominal;// klu scales the matrix entries already
 #endif
 
 };

SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.cpp

@@ -31,6 +31,7 @@ DgesvSolver::DgesvSolver(ILinearAlgLoop* algLoop, ILinSolverSettings* settings)
 	, _zeroVec            (NULL)
 	, _iterationStatus    (CONTINUE)
 	, _firstCall          (true)
+    , _fNominal (NULL)
 {
 }
 
@@ -44,6 +45,7 @@ DgesvSolver::~DgesvSolver()
 	if(_A)              delete []  _A;
 	if(_ihelpArray)       delete []  _ihelpArray;
 	if(_zeroVec)          delete []  _zeroVec;
+    if (_fNominal)		  delete []    _fNominal;
 }
 
 void DgesvSolver::initialize()
@@ -70,6 +72,7 @@ void DgesvSolver::initialize()
 			if(_A)             delete []  _A;
 			if(_ihelpArray)      delete []  _ihelpArray;
 			if(_zeroVec)         delete []  _zeroVec;
+			if (_fNominal)		 delete []    _fNominal;
 
 			_y                = new double[_dimSys];
 			_y0               = new double[_dimSys];
@@ -79,6 +82,7 @@ void DgesvSolver::initialize()
 			_A              = new double[_dimSys*_dimSys];
 			_ihelpArray       = new long int[_dimSys];
 			_zeroVec          = new double[_dimSys];
+			_fNominal          = new double[_dimSys];
 
 			_algLoop->getReal(_y);
 			_algLoop->getReal(_y0);
@@ -120,6 +124,18 @@ void DgesvSolver::solve()
 
 	memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double));
 
+
+    for (int j = 0, idx = 0; j < _dimSys; j++)
+		for (int i = 0; i < _dimSys; i++, idx++)
+			_fNominal[i] = std::max(std::abs(Atemp[idx]), _fNominal[i]);
+
+	for (int j = 0, idx = 0; j < _dimSys; j++)
+		for (int i = 0; i < _dimSys; i++, idx++)
+			_A[idx] /= _fNominal[i];
+
+	for (int i = 0; i < _dimSys; i++)
+        _b[i] /= _fNominal[i];
+
 	dgesv_(&_dimSys,&dimRHS,_A,&_dimSys,_ihelpArray,_b,&_dimSys,&irtrn);
 
 	if  (irtrn != 0){

SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp

@@ -47,6 +47,7 @@ LinearSolver::LinearSolver(ILinearAlgLoop* algLoop, ILinSolverSettings* settings
 	, _firstCall          (true)
 	, _scale			  (NULL)
 	, _generateoutput (false)
+    , _fNominal (NULL)
 {
 	_sparse = _algLoop->getUseSparseFormat();
 }
@@ -63,6 +64,7 @@ LinearSolver::~LinearSolver()
 	if (_jhelpArray)       delete[]  _jhelpArray;
 	if(_zeroVec)          delete []  _zeroVec;
 	if (_scale)            delete[]  _scale;
+    if (_fNominal)		  delete []    _fNominal;
 
 	#if defined(klu)
 		if(_sparse == true)
@@ -109,6 +111,7 @@ void LinearSolver::initialize()
 			if (_jhelpArray)       delete[]  _jhelpArray;
 			if(_zeroVec)         delete []  _zeroVec;
 			if (_scale)			 delete[]  _scale;
+			if (_fNominal)		 delete []    _fNominal;
 
 			_y                = new double[_dimSys];
 			_y0               = new double[_dimSys];
@@ -120,6 +123,7 @@ void LinearSolver::initialize()
 			_jhelpArray		  = new long int[_dimSys];
 			_zeroVec          = new double[_dimSys];
 			_scale			  = new double[_dimSys];
+			_fNominal          = new double[_dimSys];
 
 			_algLoop->getReal(_y);
 			_algLoop->getReal(_y0);
@@ -204,6 +208,28 @@ void LinearSolver::solve()
 
 		memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double));
 
+		for (int j = 0, idx = 0; j < _dimSys; j++){
+			for (int i = 0; i < _dimSys; i++, idx++){
+				_fNominal[i] = std::max(std::abs(Atemp[idx]), _fNominal[i]);
+			}
+		}
+
+		for (int i=0;i<_dimSys;i++){
+			if (_fNominal[i]==0.0){
+				_fNominal[i]==1.0;// if the row contains only zeros, there is no need to scale that row.
+			}
+		}
+
+
+		for (int j = 0, idx = 0; j < _dimSys; j++)
+			for (int i = 0; i < _dimSys; i++, idx++)
+				_A[idx] /= _fNominal[i];
+
+		for (int i = 0; i < _dimSys; i++)
+			_b[i] /= _fNominal[i];
+
+
+
 		if(_generateoutput){
 			std::cout << std::endl;
 			std::cout << "We solve a linear system with coefficient matrix" << std::endl;