Skip to content

Commit

Permalink
Linear solvers: Check residuals and fixes for Lis
Browse files Browse the repository at this point in the history 
  - Fixed wrong index in analytic Jacobian lead index for Lis solver
  - Switched column and rows in getAnalytical Jacobian for Lis
  - Added check for all linear solver using jacobian matrices for residual vector equals null
  - Added missing linear solver for linear solver tests
  • Loading branch information
@AnHeuermann @lochel
AnHeuermann authored and lochel committed Jul 1, 2019
1 parent 4dc44f9 commit ea1e2c1
Show file tree
Hide file tree
Showing 10 changed files with 189 additions and 28 deletions.
16 changes: 15 additions & 1 deletion OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c
View file
Expand Up @@ -41,6 +41,7 @@
#include "simulation_data.h"
#include "simulation/simulation_info_json.h"
#include "util/omc_error.h"
#include "omc_math.h"
#include "util/varinfo.h"
#include "model_help.h"

Expand Down Expand Up @@ -182,6 +183,7 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
void *dataAndThreadData[2] = {data, threadData};
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
DATA_KLU* solverData = (DATA_KLU*)systemData->solverData[0];
_omc_scalar residualNorm = 0;

int i, j, status = 0, success = 0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
double tmpJacEvalTime;
Expand Down Expand Up @@ -299,14 +301,26 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)

/* update inner equations */
residual_wrapper(aux_x, solverData->work, dataAndThreadData, sysNumber);
residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);

if ((isnan(residualNorm)) || (residualNorm>1e-4)){
warningStreamPrint(LOG_LS, 0,
"Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
(int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
success = 0;
}
} else {
/* the solution is automatically in x */
memcpy(aux_x, systemData->b, sizeof(double)*systemData->size);
}

if (ACTIVE_STREAM(LOG_LS_V))
{
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
if (1 == systemData->method) {
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
} else {
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
}
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);

for(i = 0; i < systemData->size; ++i)
Expand Down
6 changes: 5 additions & 1 deletion OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c
View file
Expand Up @@ -304,7 +304,11 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux
}

if (ACTIVE_STREAM(LOG_LS_V)){
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
if (1 == systemData->method) {
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
} else {
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
}
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);

for(i = 0; i < systemData->size; ++i) {
Expand Down
55 changes: 36 additions & 19 deletions OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.c
View file
Expand Up @@ -38,6 +38,7 @@
#include "simulation_data.h"
#include "simulation/simulation_info_json.h"
#include "util/omc_error.h"
#include "omc_math.h"
#include "util/varinfo.h"
#include "model_help.h"

Expand Down Expand Up @@ -91,8 +92,7 @@ allocateLisData(int n_row, int n_col, int nz, void** voiddata)
/*! \fn free memory for linear system solver Lis
*
*/
int
freeLisData(void **voiddata)
int freeLisData(void **voiddata)
{
DATA_LIS* data = (DATA_LIS*) *voiddata;

Expand All @@ -107,18 +107,26 @@ freeLisData(void **voiddata)
}


/*!
* Print LIS_MATRIX provided in column sparse row format
*
* \param [in] A Matrix A of linear problem A*x=b in CSR format
* \param [in] n Dimension of A
*/
void printLisMatrixCSR(LIS_MATRIX A, int n)
{
int i, j;
/* A matrix */
infoStreamPrint(LOG_LS_V, 1, "A matrix [%dx%d] nnz = %d ", n, n, A->nnz);
infoStreamPrint(LOG_LS_V, 1, "A matrix [%dx%d] nnz = %d", n, n, A->nnz);
infoStreamPrint(LOG_LS_V, 0, "Column Sparse Row format. Print tuple (index,value) for each row:");
for(i=0; i<n; i++)
{
char *buffer = (char*)malloc(sizeof(char)*A->ptr[i+1]*50);
buffer[0] = 0;
for(j=A->ptr[i]; j<A->ptr[i+1]; j++)
sprintf(buffer, "column %d: ", i);
for(j = A->ptr[i]; j < A->ptr[i+1]; j++)
{
sprintf(buffer, "%s(%d,%d,%g) ", buffer, i, A->index[j], A->value[j]);
sprintf(buffer, "%s(%d,%g) ", buffer, A->index[j], A->value[j]);
}
infoStreamPrint(LOG_LS_V, 0, "%s", buffer);
free(buffer);
Expand Down Expand Up @@ -158,12 +166,12 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber
{
if(jacobian->seedVars[j] == 1)
{
ii = jacobian->sparsePattern.leadindex[j-1];
ii = jacobian->sparsePattern.leadindex[j];
while(ii < jacobian->sparsePattern.leadindex[j+1])
{
l = jacobian->sparsePattern.index[ii];
/*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -jacobian->resultVars[l]); */
systemData->setAElement(i, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
/* infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", l, i, nth, -jacobian->resultVars[l]); */
systemData->setAElement(l, i, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
nth++;
ii++;
};
Expand All @@ -175,7 +183,7 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber
return 0;
}

/*! \fn wrapper_fvec_umfpack for the residual function
/*! \fn wrapper_fvec_lis for the residual function
*
*/
static int wrapper_fvec_lis(double* x, double* f, void** data, int sysNumber)
Expand All @@ -193,16 +201,15 @@ static int wrapper_fvec_lis(double* x, double* f, void** data, int sysNumber)
* [sysNumber] index of the corresponding linear system
*
*/

int
solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
{
void *dataAndThreadData[2] = {data, threadData};
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
DATA_LIS* solverData = (DATA_LIS*)systemData->solverData[0];
int i, ret, success = 1, ni, iflag = 1, n = systemData->size, eqSystemNumber = systemData->equationIndex;
char *lis_returncode[] = {"LIS_SUCCESS", "LIS_ILL_OPTION", "LIS_BREAKDOWN", "LIS_OUT_OF_MEMORY", "LIS_MAXITER", "LIS_NOT_IMPLEMENTED", "LIS_ERR_FILE_IO"};
LIS_INT err;
_omc_scalar residualNorm = 0;

int indexes[2] = {1,eqSystemNumber};
double tmpJacEvalTime;
Expand All @@ -216,10 +223,10 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
}

rt_ext_tp_tick(&(solverData->timeClock));

lis_matrix_set_size(solverData->A, solverData->n_row, 0);
if (0 == systemData->method)
{

lis_matrix_set_size(solverData->A, solverData->n_row, 0);
/* set A matrix */
systemData->setA(data, threadData, systemData);
lis_matrix_assemble(solverData->A);
Expand All @@ -228,8 +235,6 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
systemData->setb(data, threadData, systemData);

} else {

lis_matrix_set_size(solverData->A, solverData->n_row, 0);
/* calculate jacobian -> matrix A*/
if(systemData->jacobianIndex != -1){
getAnalyticalJacobianLis(data, threadData, sysNumber);
Expand Down Expand Up @@ -280,25 +285,37 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
free(buffer);
}

/* print solution */
/* Log solution */
if (1 == success){

if (1 == systemData->method){
if (1 == systemData->method){ /* Case calculate jacobian -> matrix A*/
/* take the solution */
lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x);
for(i = 0; i < solverData->n_row; ++i)
aux_x[i] += solverData->work[i];

/* update inner equations */
wrapper_fvec_lis(aux_x, solverData->work, dataAndThreadData, sysNumber);
residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);

if ((isnan(residualNorm)) || (residualNorm>1e-4)){
warningStreamPrint(LOG_LS, 0,
"Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
(int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
success = 0;
}
} else {
/* write solution */
lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x);
}

if (ACTIVE_STREAM(LOG_LS_V))
{
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
if (1 == systemData->method) {
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
} else {
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
}
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);

for(i = 0; i < systemData->size; ++i)
Expand Down
1 change: 0 additions & 1 deletion OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.h
View file
Expand Up @@ -57,6 +57,5 @@ typedef struct DATA_LIS
int allocateLisData(int n_row, int n_col, int nz, void **data);
int freeLisData(void **data);
int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x);

#endif

12 changes: 9 additions & 3 deletions OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c
View file
Expand Up @@ -38,6 +38,7 @@
#include "../../simulation_data.h"
#include "../simulation_info_json.h"
#include "../../util/omc_error.h"
#include "omc_math.h"
#include "../../util/varinfo.h"
#include "model_help.h"

Expand Down Expand Up @@ -320,8 +321,8 @@ int freeTotalPivotData(void** voiddata)
*/
int getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, double* jac, int sysNumber)
{
int i,j,k,l,ii,currentSys = sysNumber;
LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[currentSys]);
int i,j,k,l,ii;
LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]);

const int index = systemData->jacobianIndex;
ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
Expand Down Expand Up @@ -394,6 +395,7 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double*
int eqSystemNumber = systemData->equationIndex;
int indexes[2] = {1,eqSystemNumber};
int rank;
_omc_scalar residualNorm = 0;

/* We are given the number of the linear system.
* We want to look it up among all equations. */
Expand Down Expand Up @@ -459,7 +461,11 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double*
}

if (ACTIVE_STREAM(LOG_LS_V)){
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
if (1 == systemData->method) {
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
} else {
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
}
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
for(i=0; i<systemData->size; ++i)
{
Expand Down
16 changes: 15 additions & 1 deletion OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c
View file
Expand Up @@ -41,6 +41,7 @@
#include "simulation_data.h"
#include "simulation/simulation_info_json.h"
#include "util/omc_error.h"
#include "omc_math.h"
#include "util/varinfo.h"
#include "model_help.h"

Expand Down Expand Up @@ -196,6 +197,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
void *dataAndThreadData[2] = {data, threadData};
LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
DATA_UMFPACK* solverData = (DATA_UMFPACK*)systemData->solverData[0];
_omc_scalar residualNorm = 0;

int i, j, status = UMFPACK_OK, success = 0, ni=0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
int casualTearingSet = systemData->strictTearingFunctionCall != NULL;
Expand Down Expand Up @@ -305,13 +307,25 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)

/* update inner equations */
wrapper_fvec_umfpack(aux_x, solverData->work, dataAndThreadData, sysNumber);
residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);

if ((isnan(residualNorm)) || (residualNorm>1e-4)){
warningStreamPrint(LOG_LS, 0,
"Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
(int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
success = 0;
}
} else {
/* the solution is automatically in x */
}

if (ACTIVE_STREAM(LOG_LS_V))
{
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
if (1 == systemData->method) {
infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
} else {
infoStreamPrint(LOG_LS_V, 1, "Solution x:");
}
infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);

for(i = 0; i < systemData->size; ++i)
Expand Down
19 changes: 19 additions & 0 deletions OMCompiler/SimulationRuntime/c/simulation/solver/omc_math.c
View file
Expand Up @@ -804,6 +804,25 @@ _omc_scalar _omc_euclideanVectorNorm(const _omc_vector* vec)
return sqrt(result);
}

/*! \fn _omc_scalar _omc_gen_euclideanVectorNorm(_omc_vector* vec)
*
* calculates the euclidean vector norm
*
* \param [in] [vec] !TODO: DESCRIBE ME!
*/
_omc_scalar _omc_gen_euclideanVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size)
{
_omc_size i;
_omc_scalar result = 0;
assertStreamPrint(NULL, vec_size > 0, "Vector size is greater than zero");
assertStreamPrint(NULL, NULL != vec_data, "Vector data is NULL pointer");
for (i = 0; i < vec_size; ++i) {
result += pow(fabs(vec_data[i]),2.0);
}

return sqrt(result);
}

/*! \fn _omc_scalar _omc_maximumVectorNorm(_omc_vector* vec)
*
* calculates the maximum vector norm
Expand Down
1 change: 1 addition & 0 deletions OMCompiler/SimulationRuntime/c/simulation/solver/omc_math.h
View file
Expand Up @@ -119,6 +119,7 @@ void _omc_printMatrix(_omc_matrix* mat, const char* name, const int logLevel);

/* norm functions */
_omc_scalar _omc_euclideanVectorNorm(const _omc_vector* vec);
_omc_scalar _omc_gen_euclideanVectorNorm(const _omc_scalar* vec_data, const _omc_size vec_size);
_omc_scalar _omc_maximumVectorNorm(const _omc_vector* vec);

#endif

0 comments on commit ea1e2c1

Please sign in to comment.