Skip to content

Commit

Permalink
add jacobian timings to LOG_STATS
Browse files Browse the repository at this point in the history
  • Loading branch information
@OpenModelica-Hudson
wibraun authored and OpenModelica-Hudson committed Mar 1, 2017
1 parent afb17e2 commit 0db7959
Show file tree
Hide file tree
Showing 9 changed files with 112 additions and 236 deletions.
4 changes: 2 additions & 2 deletions SimulationRuntime/c/simulation/modelinfo.c
View file
Expand Up @@ -381,7 +381,7 @@ int printModelInfo(DATA *data, threadData_t *threadData, const char *filename, c
indent(fout, 2); fprintf(fout, "<initTime>%f</initTime>\n", rt_accumulated(SIM_TIMER_INIT));
indent(fout, 2); fprintf(fout, "<eventTime>%f</eventTime>\n", rt_accumulated(SIM_TIMER_EVENT));
indent(fout, 2); fprintf(fout, "<outputTime>%f</outputTime>\n", rt_accumulated(SIM_TIMER_OUTPUT));
indent(fout, 2); fprintf(fout, "<linearizeTime>%f</linearizeTime>\n", rt_accumulated(SIM_TIMER_LINEARIZE));
indent(fout, 2); fprintf(fout, "<jacobianTime>%f</jacobianTime>\n", rt_accumulated(SIM_TIMER_JACOBIAN));
indent(fout, 2); fprintf(fout, "<totalTime>%f</totalTime>\n", rt_accumulated(SIM_TIMER_TOTAL));
indent(fout, 2); fprintf(fout, "<totalStepsTime>%f</totalStepsTime>\n", rt_total(SIM_TIMER_STEP));
indent(fout, 2); fprintf(fout, "<numStep>%d</numStep>\n", (int) rt_ncall_total(SIM_TIMER_STEP));
Expand Down Expand Up @@ -595,7 +595,7 @@ int printModelInfoJSON(DATA *data, threadData_t *threadData, const char *filenam
fprintf(fout, ",\n\"initTime\":%g",rt_accumulated(SIM_TIMER_INIT));
fprintf(fout, ",\n\"eventTime\":%g",rt_accumulated(SIM_TIMER_EVENT));
fprintf(fout, ",\n\"outputTime\":%g",rt_accumulated(SIM_TIMER_OUTPUT));
fprintf(fout, ",\n\"linearizeTime\":%g",rt_accumulated(SIM_TIMER_LINEARIZE));
fprintf(fout, ",\n\"jacobianTime\":%g",rt_accumulated(SIM_TIMER_JACOBIAN));
fprintf(fout, ",\n\"totalTime\":%g",rt_accumulated(SIM_TIMER_TOTAL));
fprintf(fout, ",\n\"totalStepsTime\":%g",rt_accumulated(SIM_TIMER_STEP));
fprintf(fout, ",\n\"totalTimeProfileBlocks\":%g",totalTimeEqs); /* The overhead the profiling is huge if small equations are profiled */
Expand Down
4 changes: 2 additions & 2 deletions SimulationRuntime/c/simulation/simulation_runtime.cpp
View file
Expand Up @@ -588,9 +588,9 @@ int startNonInteractiveSimulation(int argc, char**argv, DATA* data, threadData_t
}

if(0 == retVal && create_linearmodel) {
rt_tick(SIM_TIMER_LINEARIZE);
rt_tick(SIM_TIMER_JACOBIAN);
retVal = linearize(data, threadData);
rt_accumulate(SIM_TIMER_LINEARIZE);
rt_accumulate(SIM_TIMER_JACOBIAN);
infoStreamPrint(LOG_STDOUT, 0, "Linear model is created!");
}

Expand Down
230 changes: 69 additions & 161 deletions SimulationRuntime/c/simulation/solver/dassl.c
View file
Expand Up @@ -73,14 +73,16 @@ dummy_zeroCrossing(int *neqm, double *t, double *y, double *yp,
return 0;
}

static int JacobianSymbolic(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int JacobianSymbolicColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int JacobianOwnNum(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int JacobianOwnNumColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int callJacobian(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int jacA_num(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int jacA_numColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int jacA_sym(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
static int jacA_symColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);

void DDASKR(
int (*res) (double *t, double *y, double *yprime, double* cj, double *delta, int *ires, double *rpar, int* ipar),
Expand Down Expand Up @@ -345,17 +347,17 @@ int dassl_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo,
switch (dasslData->dasslJacobian){
case COLOREDNUMJAC:
data->simulationInfo->jacobianEvals = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern.maxColors;
dasslData->jacobianFunction = JacobianOwnNumColored;
dasslData->jacobianFunction = jacA_numColored;
break;
case COLOREDSYMJAC:
data->simulationInfo->jacobianEvals = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern.maxColors;
dasslData->jacobianFunction = JacobianSymbolicColored;
dasslData->jacobianFunction = jacA_symColored;
break;
case SYMJAC:
dasslData->jacobianFunction = JacobianSymbolic;
dasslData->jacobianFunction = jacA_sym;
break;
case NUMJAC:
dasslData->jacobianFunction = JacobianOwnNum;
dasslData->jacobianFunction = jacA_num;
break;
case INTERNALNUMJAC:
dasslData->jacobianFunction = dummy_Jacobian;
Expand Down Expand Up @@ -603,7 +605,7 @@ int dassl_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
&solverInfo->currentTime, states, stateDer, &tout,
dasslData->info, dasslData->rtol, dasslData->atol, &dasslData->idid,
dasslData->rwork, &dasslData->lrw, dasslData->iwork, &dasslData->liw,
(double*) (void*) dasslData->rpar, dasslData->ipar, dasslData->jacobianFunction, dummy_precondition,
(double*) (void*) dasslData->rpar, dasslData->ipar, callJacobian, dummy_precondition,
dasslData->zeroCrossingFunction, (int*) &dasslData->ng, dasslData->jroot);

/* closing new step message */
Expand Down Expand Up @@ -969,10 +971,19 @@ int function_ZeroCrossingsDASSL(int *neqm, double *t, double *y, double *yp,
return 0;
}


int functionJacAColored(DATA* data, threadData_t *threadData, double* jac)
/* \fn jacA_symColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
*
*
* This function calculates symbolically the jacobian matrix and exploiting the coloring.
*/
int jacA_symColored(double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, int *ipar)
{
TRACE_PUSH
DATA* data = (DATA*)(void*)((double**)rpar)[0];
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];

const int index = data->callback->INDEX_JAC_A;
unsigned int i,j,l,k,ii;

Expand All @@ -993,7 +1004,7 @@ int functionJacAColored(DATA* data, threadData_t *threadData, double* jac)
{
l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l;
jac[k] = data->simulationInfo->analyticJacobians[index].resultVars[l];
matrixA[k] = data->simulationInfo->analyticJacobians[index].resultVars[l];
ii++;
};
}
Expand All @@ -1007,10 +1018,18 @@ int functionJacAColored(DATA* data, threadData_t *threadData, double* jac)
return 0;
}


int functionJacASym(DATA* data, threadData_t *threadData, double* jac)
/* \fn jacA_sym(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
*
*
* This function calculates symbolically the jacobian matrix.
*/
int jacA_sym(double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, int *ipar)
{
TRACE_PUSH
DATA* data = (DATA*)(void*)((double**)rpar)[0];
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
const int index = data->callback->INDEX_JAC_A;
unsigned int i,j,k;

Expand All @@ -1023,7 +1042,7 @@ int functionJacASym(DATA* data, threadData_t *threadData, double* jac)

for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeRows; j++)
{
jac[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j];
matrixA[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j];
}

data->simulationInfo->analyticJacobians[index].seedVars[i] = 0.0;
Expand All @@ -1033,105 +1052,19 @@ int functionJacASym(DATA* data, threadData_t *threadData, double* jac)
return 0;
}

/*
* provides a analytical Jacobian to be used with DASSL
*/

static int JacobianSymbolicColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
{
TRACE_PUSH
DATA* data = (DATA*)(void*)((double**)rpar)[0];
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
double* backupStates;
double timeBackup;
int i;
int j;

setContext(data, t, CONTEXT_JACOBIAN);

backupStates = data->localData[0]->realVars;
timeBackup = data->localData[0]->timeValue;

data->localData[0]->timeValue = *t;
data->localData[0]->realVars = y;
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data, threadData);
/* eval ode*/
data->callback->functionODE(data, threadData);
functionJacAColored(data, threadData, pd);

/* add cj to the diagonal elements of the matrix */
j = 0;
for(i = 0; i < data->modelData->nStates; i++)
{
pd[j] -= (double) *cj;
j += data->modelData->nStates + 1;
}
data->localData[0]->realVars = backupStates;
data->localData[0]->timeValue = timeBackup;

unsetContext(data);

TRACE_POP
return 0;
}

/*
* provides a analytical Jacobian to be used with DASSL
/* \fn jacA_num(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
*
*
* This function calculates a jacobian matrix by
* numerical with forward finite differences.
*/
static int JacobianSymbolic(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
int jacA_num(double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, int *ipar)
{
TRACE_PUSH
DATA* data = (DATA*)(void*)((double**)rpar)[0];
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];
double* backupStates;
double timeBackup;
int i;
int j;

setContext(data, t, CONTEXT_JACOBIAN);

backupStates = data->localData[0]->realVars;
timeBackup = data->localData[0]->timeValue;

data->localData[0]->timeValue = *t;
data->localData[0]->realVars = y;
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data, threadData);
/* eval ode*/
data->callback->functionODE(data, threadData);
functionJacASym(data, threadData, pd);

/* add cj to the diagonal elements of the matrix */
j = 0;
for(i = 0; i < data->modelData->nStates; i++)
{
pd[j] -= (double) *cj;
j += data->modelData->nStates + 1;
}
data->localData[0]->realVars = backupStates;
data->localData[0]->timeValue = timeBackup;

unsetContext(data);

TRACE_POP
return 0;
}

/*
* function calculates a jacobian matrix by
* numerical method finite differences
*/
int jacA_num(DATA* data, double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, int *ipar)
{
TRACE_PUSH
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];

double delta_h = numericalDifferentiationDeltaXsolver;
double delta_hh,delta_hhh, deltaInv;
Expand Down Expand Up @@ -1179,59 +1112,22 @@ int jacA_num(DATA* data, double *t, double *y, double *yprime, double *delta, do
return 0;
}

/*
* provides a numerical Jacobian to be used with DASSL
/* \fn jacA_numColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
*
*
* This function calculates a jacobian matrix by
* numerical with forward finite differences and exploiting the coloring.
*/
static int JacobianOwnNum(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
int jacA_numColored(double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, int *ipar)
{
TRACE_PUSH

DATA* data = (DATA*)(void*)((double**)rpar)[0];
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];

setContext(data, t, CONTEXT_JACOBIAN);

if(jacA_num(data, t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
{
throwStreamPrint(threadData, "Error, can not get Matrix A ");
TRACE_POP
return 1;
}

/* debug */
if (ACTIVE_STREAM(LOG_JAC)){
_omc_matrix* dumpJac = _omc_createMatrix(dasslData->N, dasslData->N, pd);
_omc_printMatrix(dumpJac, "DASSL-Solver: Matrix A", LOG_JAC);
_omc_destroyMatrix(dumpJac);
}


/* add cj to diagonal elements and store in pd */
if (!dasslData->daeMode){
int i, j = 0;
for(i = 0; i < dasslData->N; i++)
{
pd[j] -= (double) *cj;
j += dasslData->N + 1;
}
}
unsetContext(data);

TRACE_POP
return 0;
}


/*
* function calculates a jacobian matrix by
* numerical method finite differences
*/
int jacA_numColored(DATA* data, double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, int *ipar)
{
TRACE_PUSH
const int index = data->callback->INDEX_JAC_A;
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
double delta_h = numericalDifferentiationDeltaXsolver;
double delta_hhh;
int ires;
Expand Down Expand Up @@ -1294,26 +1190,37 @@ int jacA_numColored(DATA* data, double *t, double *y, double *yprime, double *de
return 0;
}

/*
* provides a numerical Jacobian to be used with DASSL
/* \fn callJacobian(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
*
*
* This function is called by dassl to calculate the jacobian matrix.
*
*/
static int JacobianOwnNumColored(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
static int callJacobian(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar)
{
TRACE_PUSH
DATA* data = (DATA*)(void*)((double**)rpar)[0];
DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1];
threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2];

/* set context for the start values extrapolation of non-linear algebraic loops */
setContext(data, t, CONTEXT_JACOBIAN);

if(jacA_numColored(data, t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
/* profiling */
rt_tick(SIM_TIMER_JACOBIAN);

if(dasslData->jacobianFunction(t, y, yprime, deltaD, pd, cj, h, wt, rpar, ipar))
{
throwStreamPrint(threadData, "Error, can not get Matrix A ");
TRACE_POP
return 1;
}

/* profiling */
rt_accumulate(SIM_TIMER_JACOBIAN);

/* debug */
if (ACTIVE_STREAM(LOG_JAC)){
_omc_matrix* dumpJac = _omc_createMatrix(dasslData->N, dasslData->N, pd);
Expand All @@ -1331,6 +1238,7 @@ static int JacobianOwnNumColored(double *t, double *y, double *yprime, double *d
j += dasslData->N + 1;
}
}
/* set context for the start values extrapolation of non-linear algebraic loops */
unsetContext(data);

TRACE_POP
Expand Down
5 changes: 3 additions & 2 deletions SimulationRuntime/c/simulation/solver/dassl.h
View file
Expand Up @@ -74,9 +74,10 @@ typedef struct DASSL_DATA{
double *stateDer;
double *states;

/* function pointer of provied functions */
/* function pointer of provided functions */
int (*residualFunction)(double *t, double *x, double *xprime, double *cj, double *delta, int *ires, double *rpar, int* ipar);
void* jacobianFunction;
int (*jacobianFunction)(double *t, double *y, double *yprime, double *deltaD, double *pd, double *cj, double *h, double *wt,
double *rpar, int* ipar);
void* zeroCrossingFunction;
} DASSL_DATA;

Expand Down

0 comments on commit 0db7959

Please sign in to comment.