Skip to content

Commit

Permalink
extends solver stats for all solver
Browse files Browse the repository at this point in the history
  • Loading branch information
Willi Braun committed Apr 6, 2016
1 parent e8fbdf6 commit df3f1c9
Show file tree
Hide file tree
Showing 7 changed files with 85 additions and 44 deletions.
8 changes: 1 addition & 7 deletions SimulationRuntime/c/simulation/solver/dassl.c
View file
Expand Up @@ -159,10 +159,6 @@ int dassl_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo,
dasslData->rtol = (double*) malloc(data->modelData->nStates*sizeof(double));
dasslData->info = (int*) calloc(infoLength, sizeof(int));
assertStreamPrint(threadData, 0 != dasslData->info,"out of memory");
dasslData->dasslStatistics = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));
assertStreamPrint(threadData, 0 != dasslData->dasslStatistics,"out of memory");
dasslData->dasslStatisticsTmp = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));
assertStreamPrint(threadData, 0 != dasslData->dasslStatisticsTmp,"out of memory");

dasslData->idid = 0;

Expand Down Expand Up @@ -407,8 +403,6 @@ int dassl_deinitial(DASSL_DATA *dasslData)
free(dasslData->delta_hh);
free(dasslData->newdelta);
free(dasslData->stateDer);
free(dasslData->dasslStatistics);
free(dasslData->dasslStatisticsTmp);

free(dasslData);

Expand Down Expand Up @@ -637,7 +631,7 @@ int dassl_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
for(ui = 0; ui < numStatistics; ui++)
{
assert(10 + ui < dasslData->liw);
dasslData->dasslStatisticsTmp[ui] = dasslData->iwork[10 + ui];
solverInfo->solverStatsTmp[ui] = dasslData->iwork[10 + ui];
}

infoStreamPrint(LOG_DASSL, 0, "Finished DASSL step.");
Expand Down
4 changes: 0 additions & 4 deletions SimulationRuntime/c/simulation/solver/dassl.h
View file
Expand Up @@ -36,7 +36,6 @@
#define DDASKR _daskr_ddaskr_

static const unsigned int maxOrder = 5;
static const unsigned int numStatistics = 5;
static const unsigned int infoLength = 20;

enum DASSL_JACOBIAN
Expand All @@ -59,9 +58,6 @@ typedef struct DASSL_DATA{
int dasslJacobian; /* specifices the method to calculate the jacobian matrix */
int dasslAvoidEventRestart; /* if TRUE then no restart after an event is performed */

unsigned int* dasslStatistics;
unsigned int* dasslStatisticsTmp;

int* info;

int idid;
Expand Down
8 changes: 4 additions & 4 deletions SimulationRuntime/c/simulation/solver/perform_simulation.c
View file
Expand Up @@ -292,14 +292,14 @@ static void retrySimulationStep(DATA* data, threadData_t *threadData, SOLVER_INF
solverInfo->didEventStep = 1;
}

static void saveDasslStats(SOLVER_INFO* solverInfo)
static void saveIntegratorStats(SOLVER_INFO* solverInfo)
{
int ui;
if (solverInfo->didEventStep == 1 && solverInfo->solverMethod == S_DASSL)
if (solverInfo->didEventStep == 1)
{
for(ui=0; ui<numStatistics; ui++)
{
((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[ui] += ((DASSL_DATA*)solverInfo->solverData)->dasslStatisticsTmp[ui];
solverInfo->solverStats[ui] += solverInfo->solverStatsTmp[ui];
}
}
}
Expand Down Expand Up @@ -407,7 +407,7 @@ int prefixedName_performSimulation(DATA* data, threadData_t *threadData, SOLVER_
retry = 0; /* reset retry */

fmtEmitStep(data, threadData, &fmt, solverInfo->didEventStep);
saveDasslStats(solverInfo);
saveIntegratorStats(solverInfo);
checkSimulationTerminated(data, solverInfo);

/* terminate for some cases:
Expand Down
46 changes: 40 additions & 6 deletions SimulationRuntime/c/simulation/solver/radau.c
View file
Expand Up @@ -652,31 +652,37 @@ static int lobatto6Res(N_Vector x, N_Vector f, void* user_data)
return 0;
}

int kinsolOde(void* ode)
int kinsolOde(SOLVER_INFO* solverInfo)
{
KINODE *kinOde = (KINODE*) ode;
KINODE *kinOde = (KINODE*) solverInfo->solverData;
KDATAODE *kData = kinOde->kData;
int i;
int i, flag, dense=0;
long int tmp;
initKinsol(kinOde);
for(i = 0;i <3; ++i)
for(i = 0; i < 3; ++i)
{

kData->error_code = KINSol(kData->kmem, /* KINSol memory block */
kData->x, /* initial guess on input; solution vector */
kData->glstr, /* global stragegy choice */
kData->glstr, /* global strategy choice */
kData->sVars, /* scaling vector, for the variable cc */
kData->sEqns );

if(kData->error_code>=0)
return 0;
{
break;
}

if(i == 0)
{
KINDense(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates);
dense=1;
infoStreamPrint(LOG_SOLVER,0,"Restart Kinsol: change linear solver to KINDense.");
}
else if(i == 1)
{
KINSptfqmr(kinOde->kData->kmem, kinOde->N*kinOde->nlp->nStates);
dense=0;
infoStreamPrint(LOG_SOLVER,0,"Restart Kinsol: change linear solver to KINSptfqmr.");
}
else if(i == 2)
Expand All @@ -685,6 +691,34 @@ int kinsolOde(void* ode)
infoStreamPrint(LOG_SOLVER,0,"Restart Kinsol: change linear solver to KINSpbcg.");
}
}
/* save stats */
/* steps */
solverInfo->solverStatsTmp[0] += 1;
/* functionODE evaluations */
tmp = 0;
flag = KINGetNumFuncEvals(kData->kmem, &tmp);
if (flag == KIN_SUCCESS)
{
solverInfo->solverStatsTmp[1] += tmp;
}

/* Jacobians evaluations */
tmp = 0;
flag = KINDlsGetNumJacEvals(kData->kmem, &tmp);
if (flag == KIN_SUCCESS)
{
solverInfo->solverStatsTmp[2] += tmp;
}

/* beta-condition failures evaluations */
tmp = 0;
flag = KINGetNumBetaCondFails(kData->kmem, &tmp);
if (flag == KIN_SUCCESS)
{
solverInfo->solverStatsTmp[4] += tmp;
}


return (kData->error_code<0) ? -1 : 0;
}
#else
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/c/simulation/solver/radau.h
View file
Expand Up @@ -100,7 +100,7 @@
#endif /* SUNDIALS */
int allocateKinOde(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, int flag, int N);
int freeKinOde(DATA* data, SOLVER_INFO* solverInfo, int N);
int kinsolOde(void* ode);
int kinsolOde(SOLVER_INFO* solverInfo);
#ifdef __cplusplus
};
#endif
Expand Down
56 changes: 34 additions & 22 deletions SimulationRuntime/c/simulation/solver/solver_main.c
View file
Expand Up @@ -185,6 +185,8 @@ int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solv
solverInfo->didEventStep = 0;
solverInfo->stateEvents = 0;
solverInfo->sampleEvents = 0;
solverInfo->solverStats = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));
solverInfo->solverStatsTmp = (unsigned int*) calloc(numStatistics, sizeof(unsigned int));

/* if FLAG_NOEQUIDISTANT_GRID is set, choose dassl step method */
if (omc_flag[FLAG_NOEQUIDISTANT_GRID])
Expand Down Expand Up @@ -547,31 +549,22 @@ int finishSimulation(DATA* data, threadData_t *threadData, SOLVER_INFO* solverIn
infoStreamPrint(LOG_STATS, 0, "%5ld time events", solverInfo->sampleEvents);
messageClose(LOG_STATS);

#if defined(WITH_DASSL)
if(S_DASSL == solverInfo->solverMethod)
if(S_OPTIMIZATION == solverInfo->solverMethod || /* skip solver statistics for optimization */
S_QSS == solverInfo->solverMethod) /* skip also for qss, since not available*/
{
/* save dassl stats before print */
for(ui=0; ui<numStatistics; ui++)
((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[ui] += ((DASSL_DATA*)solverInfo->solverData)->dasslStatisticsTmp[ui];

infoStreamPrint(LOG_STATS, 1, "solver: DASSL");
infoStreamPrint(LOG_STATS, 0, "%5d steps taken", ((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[0]);
infoStreamPrint(LOG_STATS, 0, "%5d calls of functionODE", ((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[1]);
infoStreamPrint(LOG_STATS, 0, "%5d evaluations of jacobian", ((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[2]);
infoStreamPrint(LOG_STATS, 0, "%5d error test failures", ((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[3]);
infoStreamPrint(LOG_STATS, 0, "%5d convergence test failures", ((DASSL_DATA*)solverInfo->solverData)->dasslStatistics[4]);
messageClose(LOG_STATS);
}
else
#endif
if(S_OPTIMIZATION == solverInfo->solverMethod)
{
/* skip solver statistics for optimization */
}
else
{
infoStreamPrint(LOG_STATS, 1, "solver");
infoStreamPrint(LOG_STATS, 0, "sorry - no solver statistics available. [not yet implemented]");
/* save stats before print */
for(ui=0; ui<numStatistics; ui++)
solverInfo->solverStats[ui] += solverInfo->solverStatsTmp[ui];

infoStreamPrint(LOG_STATS, 1, "solver: %s", SOLVER_METHOD_NAME[solverInfo->solverMethod]);
infoStreamPrint(LOG_STATS, 0, "%5d steps taken", solverInfo->solverStats[0]);
infoStreamPrint(LOG_STATS, 0, "%5d calls of functionODE", solverInfo->solverStats[1]);
infoStreamPrint(LOG_STATS, 0, "%5d evaluations of jacobian", solverInfo->solverStats[2]);
infoStreamPrint(LOG_STATS, 0, "%5d error test failures", solverInfo->solverStats[3]);
infoStreamPrint(LOG_STATS, 0, "%5d convergence test failures", solverInfo->solverStats[4]);
messageClose(LOG_STATS);
}

Expand Down Expand Up @@ -749,6 +742,12 @@ static int euler_ex_step(DATA* data, SOLVER_INFO* solverInfo)
}
sData->timeValue = solverInfo->currentTime;

/* save stats */
/* steps */
solverInfo->solverStatsTmp[0] += 1;
/* function ODE evaluation is done directly after this function */
solverInfo->solverStatsTmp[1] += 1;

return 0;
}

Expand All @@ -772,6 +771,12 @@ static int sym_euler_im_step(DATA* data, threadData_t *threadData, SOLVER_INFO*
/* update der(x)*/
for(i=0, j=data->modelData->nStates; i<data->modelData->nStates; ++i, ++j)
data->localData[0]->realVars[j] = (data->localData[0]->realVars[i]-data->localData[1]->realVars[i])/solverInfo->currentStepSize;

/* save stats */
/* steps */
solverInfo->solverStatsTmp[0] += 1;
/* function ODE evaluation is done directly after this */
solverInfo->solverStatsTmp[1] += 1;
return retVal;
}

Expand Down Expand Up @@ -824,6 +829,13 @@ static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* so
sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * sum;
}
sData->timeValue = solverInfo->currentTime;

/* save stats */
/* steps */
solverInfo->solverStatsTmp[0] += 1;
/* function ODE evaluation is done directly after this */
solverInfo->solverStatsTmp[1] += rk->work_states_ndims+1;

return 0;
}

Expand All @@ -845,7 +857,7 @@ static int ipopt_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverI
/*************************************** Radau/Lobatto ***********************************/
int radau_lobatto_step(DATA* data, SOLVER_INFO* solverInfo)
{
if(kinsolOde(solverInfo->solverData) == 0)
if(kinsolOde(solverInfo) == 0)
{
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
Expand Down
5 changes: 5 additions & 0 deletions SimulationRuntime/c/simulation/solver/solver_main.h
View file
Expand Up @@ -42,6 +42,8 @@
#include "util/list.h"
#include "util/simulation_options.h"

static const unsigned int numStatistics = 5;

typedef struct SOLVER_INFO
{
double currentTime;
Expand All @@ -66,6 +68,9 @@ typedef struct SOLVER_INFO
/* stats */
unsigned long stateEvents;
unsigned long sampleEvents;
/* integrator stats */
unsigned int* solverStats;
unsigned int* solverStatsTmp;

/* further options */
int integratorSteps;
Expand Down

0 comments on commit df3f1c9

Please sign in to comment.