Skip to content

Commit 9dde075

Browse files
linuslangenkamplinuslangenkamp
authored
[C] unify Jacobian evaluation (#13971)
1 parent 5f7cb2c commit 9dde075

File tree

12 files changed

+79
-58
lines changed

12 files changed

+79
-58
lines changed

‎OMCompiler/SimulationRuntime/c/simulation/jacobian_util.c‎

Lines changed: 31 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -106,51 +106,62 @@ void freeJacobian(JACOBIAN *jac)
106106

107107
/*! \fn evalJacobian
108108
*
109-
* compute entries of Jacobian
110-
* uses coloring (from sparsity pattern) if available
109+
* compute entries of Jacobian in sparse CSC or dense format
110+
* uses coloring (sparsePattern non NULL)
111111
*
112112
* \param [ref] [data]
113113
* \param [ref] [threadData]
114-
* \param [ref] [jacobian] Pointer to Jacobian
115-
* \param [out] [jac] Contains jacobian values on exit, column-major
114+
* \param [ref] [jacobian] Pointer to Jacobian
115+
* \param [ref] [parentJacobian] Pointer to parent Jacobian
116+
* \param [out] [jac] Output buffer, size nnz (sparse) or #rows * #cols (dense), non zero-initialized
117+
* \param [ref] [isDense] Flag to set dense / sparse output
116118
*/
117-
void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac)
119+
void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac, modelica_boolean isDense)
118120
{
119-
int i,j,k,l,ii;
121+
int color, column, row, nz;
120122
const SPARSE_PATTERN* sp = jacobian->sparsePattern;
121123

122-
memset(jac, 0.0, (jacobian->sizeRows) * (jacobian->sizeCols) * sizeof(modelica_real));
123-
124124
/* evaluate constant equations of Jacobian */
125125
if (jacobian->constantEqns != NULL) {
126126
jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
127127
}
128128

129+
if (isDense) {
130+
/* memset to zero for dense, since solvers might destroy "hard zeros"
131+
* does not apply for sparse, since the values are overwritten */
132+
memset(jac, 0.0, jacobian->sizeRows * jacobian->sizeCols * sizeof(modelica_real));
133+
}
134+
129135
/* evaluate Jacobian */
130-
for (i = 0; i < sp->maxColors; i++) {
136+
for (color = 0; color < sp->maxColors; color++) {
131137
/* activate seed variable for the corresponding color */
132-
for (j = 0; j < jacobian->sizeCols; j++)
133-
if (sp->colorCols[j]-1 == i)
134-
jacobian->seedVars[j] = 1.0;
138+
for (column = 0; column < jacobian->sizeCols; column++)
139+
if (sp->colorCols[column]-1 == color)
140+
jacobian->seedVars[column] = 1.0;
135141

136142
/* evaluate Jacobian column */
137143
jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
138144

139-
for (j = 0; j < jacobian->sizeCols; j++) {
140-
if (sp->colorCols[j]-1 == i) {
141-
for (ii = sp->leadindex[j]; ii < sp->leadindex[j+1]; ii++) {
142-
l = sp->index[ii];
143-
k = j*jacobian->sizeRows + l;
144-
jac[k] = jacobian->resultVars[l]; //* solverData->xScaling[j];
145+
for (column = 0; column < jacobian->sizeCols; column++) {
146+
if (sp->colorCols[column]-1 == color) {
147+
for (nz = sp->leadindex[column]; nz < sp->leadindex[column+1]; nz++) {
148+
row = sp->index[nz];
149+
if (!isDense) {
150+
/* sparse case */
151+
jac[nz] = jacobian->resultVars[row]; //* solverData->xScaling[j];
152+
}
153+
else {
154+
/* dense case */
155+
jac[column * jacobian->sizeRows + row] = jacobian->resultVars[row]; //* solverData->xScaling[j];
156+
}
145157
}
146158
/* de-activate seed variable for the corresponding color */
147-
jacobian->seedVars[j] = 0.0;
159+
jacobian->seedVars[column] = 0.0;
148160
}
149161
}
150162
}
151163
}
152164

153-
154165
/**
155166
* @brief Allocate memory for sparsity pattern.
156167
*

‎OMCompiler/SimulationRuntime/c/simulation/jacobian_util.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ void initJacobian(JACOBIAN* jacobian, unsigned int sizeCols, unsigned int sizeRo
4444
JACOBIAN* copyJacobian(JACOBIAN* source);
4545
void freeJacobian(JACOBIAN* jac);
4646

47-
void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac);
47+
void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac, modelica_boolean isDense);
4848

4949
SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int numberOfNonZeros, unsigned int maxColors);
5050
void freeSparsePattern(SPARSE_PATTERN *spp);

‎OMCompiler/SimulationRuntime/c/simulation/solver/kinsolSolver.c‎

Lines changed: 9 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -600,7 +600,7 @@ int nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
600600
const SPARSE_PATTERN* sp = jacobian->sparsePattern;
601601
assertStreamPrint(threadData, NULL != sp, "sp is NULL");
602602
double *xScaling = NV_DATA_S(kinsolData->xScale);
603-
long int i, j, ii, l;
603+
long int column, nz;
604604

605605
if (SUNMatGetID(Jac) != SUNMATRIX_SPARSE || SM_SPARSETYPE_S(Jac) == CSR_MAT) {
606606
errorStreamPrint(OMC_LOG_STDOUT, 0,
@@ -611,36 +611,15 @@ int nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SUNMatrix Jac,
611611
/* performance measurement */
612612
rt_ext_tp_tick(&nlsData->jacobianTimeClock);
613613

614-
/* reset matrix */
615-
SUNMatZero(Jac);
616-
617-
/* Evaluate constant equations of Jacobian */
618-
if (jacobian->constantEqns != NULL) {
619-
jacobian->constantEqns(data, threadData, jacobian, NULL);
620-
}
621-
622-
/* Evaluate Jacobian */
623-
for (i = 0; i < sp->maxColors; i++) {
624-
/* Set seed variables */
625-
for (j = 0; j < jacobian->sizeCols; j++)
626-
if (sp->colorCols[j] - 1 == i)
627-
jacobian->seedVars[j] = 1.0;
628-
629-
/* Evaluate Jacobian column */
630-
jacobian->evalColumn(data, threadData, jacobian, NULL);
614+
/* call generic sparse Jacobian with CSC buffer "SM_DATA_S(Jac)" */
615+
evalJacobian(data, threadData, jacobian, NULL, SM_DATA_S(Jac), FALSE);
616+
setSundialsSparsePattern(jacobian, Jac);
631617

632-
/* Save column in Jac and unset seed variables */
633-
for (j = 0; j < jacobian->sizeCols; j++) {
634-
if (sp->colorCols[j] - 1 == i) {
635-
for (ii = sp->leadindex[j]; ii < sp->leadindex[j + 1]; ii++) {
636-
l = sp->index[ii];
637-
if (kinsolData->nominalJac) {
638-
setJacElementSundialsSparse(l, j, ii, jacobian->resultVars[l] / xScaling[j], Jac, SM_CONTENT_S(Jac)->M);
639-
} else {
640-
setJacElementSundialsSparse(l, j, ii, jacobian->resultVars[l], Jac, SM_CONTENT_S(Jac)->M);
641-
}
642-
}
643-
jacobian->seedVars[j] = 0.0;
618+
/* scaling */
619+
if (kinsolData->nominalJac) {
620+
for (column = 0; column < jacobian->sizeCols; column++) {
621+
for (nz = sp->leadindex[column]; nz < sp->leadindex[column + 1]; nz++) {
622+
SM_DATA_S(Jac)[nz] / xScaling[column];
644623
}
645624
}
646625
}

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -113,7 +113,8 @@ void getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, LINEAR_SY
113113
JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
114114
JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
115115

116-
evalJacobian(data, threadData, jacobian, parentJacobian, jac);
116+
/* call generic dense Jacobian */
117+
evalJacobian(data, threadData, jacobian, parentJacobian, jac, TRUE);
117118

118119
for (k = 0; k < (jacobian->sizeRows) * (jacobian->sizeCols); k++)
119120
jac[k] = -jac[k];

‎OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -341,7 +341,8 @@ void getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, LINEA
341341
JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
342342
JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
343343

344-
evalJacobian(data, threadData, jacobian, parentJacobian, jac);
344+
/* call generic dense Jacobian */
345+
evalJacobian(data, threadData, jacobian, parentJacobian, jac, TRUE);
345346
}
346347

347348
/*! \fn wrapper_fvec_hybrd for the residual Function

‎OMCompiler/SimulationRuntime/c/simulation/solver/newton_diagnostics.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,8 @@ double** getJacobian( DATA* data, threadData_t *threadData, NONLINEAR_SYSTEM_DAT
135135
jac = (modelica_real*) calloc(jacobian->sizeRows * jacobian->sizeCols, sizeof(modelica_real*));
136136
assertStreamPrint(threadData, NULL != jac, "out of memory");
137137

138-
evalJacobian(data, threadData, jacobian, NULL, jac);
138+
/* call generic dense Jacobian */
139+
evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
139140

140141
/* copy jacobian from column-major to row-major */
141142
for (i = 0; i < jacobian->sizeRows; ++i)

‎OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -832,7 +832,8 @@ int getAnalyticalJacobianHomotopy(DATA_HOMOTOPY* solverData, double* jac)
832832
JACOBIAN* jacobian = solverData->userData->analyticJacobian;
833833
const SPARSE_PATTERN* sp = jacobian->sparsePattern;
834834

835-
evalJacobian(data, threadData, jacobian, NULL, jac);
835+
/* call generic dense Jacobian */
836+
evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
836837

837838
/* apply scaling to each column */
838839
for (j = 0; j < jacobian->sizeCols; j++) {

‎OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -265,7 +265,8 @@ static int getAnalyticalJacobian(NLS_USERDATA* hybrdUserData, double* jac)
265265
DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData);
266266
JACOBIAN* jacobian = hybrdUserData->analyticJacobian;
267267

268-
evalJacobian(data, threadData, jacobian, NULL, jac);
268+
/* call generic dense Jacobian */
269+
evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
269270

270271
memcpy(solverData->fjacobian, jac, (jacobian->sizeRows) * (jacobian->sizeCols) * sizeof(modelica_real));
271272

‎OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,8 @@ int wrapper_fvec_newton(int n, double* x, double* fvec, NLS_USERDATA* userData,
9696
rt_ext_tp_tick(&nlsData->jacobianTimeClock);
9797

9898
if(nlsData->jacobianIndex != -1 && jacobian != NULL ) {
99-
evalJacobian(data, threadData, jacobian, NULL, solverData->fjac);
99+
/* call generic dense Jacobian */
100+
evalJacobian(data, threadData, jacobian, NULL, solverData->fjac, TRUE);
100101
} else {
101102
double delta_h = sqrt(solverData->epsfcn);
102103
double delta_hh;

‎OMCompiler/SimulationRuntime/c/simulation/solver/stateset.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -188,7 +188,8 @@ static void getAnalyticalJacobianSet(DATA* data, threadData_t *threadData, unsig
188188

189189
modelica_real* jac = data->simulationInfo->stateSetData[index].J;
190190

191-
evalJacobian(data, threadData, jacobian, NULL, jac);
191+
/* call generic dense Jacobian */
192+
evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
192193

193194
if(OMC_ACTIVE_STREAM(OMC_LOG_DSS_JAC))
194195
{

0 commit comments

Comments
 (0)