[C] Unify Jacobian evaluation (wip) (#13709)

TODO:
- sparse memory format
- unify remaining functions
- numeric jacobian
- scaling (x, residual?)

Author: phannebohm

OMCompiler/Compiler/Template/CodegenC.tpl

@@ -1489,7 +1489,7 @@ matrixData getJacobianMatrixF(DATA * data, threadData_t * threadData, ofstream &
 {
   // initialize the jacobian call
   const int index = data->callback->INDEX_JAC_F;
-  ANALYTIC_JACOBIAN *jacobian = &(data->simulationInfo->analyticJacobians[index]);
+  JACOBIAN *jacobian = &(data->simulationInfo->analyticJacobians[index]);
   data->callback->initialAnalyticJacobianF(data, threadData, jacobian);
   int cols = jacobian->sizeCols;
   int rows = jacobian->sizeRows;
@@ -1533,7 +1533,7 @@ matrixData getJacobianMatrixH(DATA * data, threadData_t * threadData, ofstream &
 {
   // initialize the jacobian call
   const int index = data->callback->INDEX_JAC_H;
-  ANALYTIC_JACOBIAN *jacobian = &(data->simulationInfo->analyticJacobians[index]);
+  JACOBIAN *jacobian = &(data->simulationInfo->analyticJacobians[index]);
   data->callback->initialAnalyticJacobianH(data, threadData, jacobian);
   int cols = jacobian->sizeCols;
   int rows = jacobian->sizeRows;

OMCompiler/SimulationRuntime/c/dataReconciliation/dataReconciliation.cpp

@@ -332,7 +332,7 @@ int functionJacBD_num(DATA* data, threadData_t *threadData, double *matrixB, dou
 int functionJacA(DATA* data, threadData_t *threadData, double* jac){
 
   const int index = data->callback->INDEX_JAC_A;
-  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
+  JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
   unsigned int i,j,k;
   k = 0;
   if (jacobian->constantEqns != NULL) {
@@ -378,7 +378,7 @@ int functionJacA(DATA* data, threadData_t *threadData, double* jac){
 int functionJacB(DATA* data, threadData_t *threadData, double* jac){
 
   const int index = data->callback->INDEX_JAC_B;
-  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
+  JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
 
   unsigned int i,j,k;
   k = 0;
@@ -424,7 +424,7 @@ int functionJacB(DATA* data, threadData_t *threadData, double* jac){
 int functionJacC(DATA* data, threadData_t *threadData, double* jac){
 
   const int index = data->callback->INDEX_JAC_C;
-  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
+  JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
   unsigned int i,j,k;
   k = 0;
   if (jacobian->constantEqns != NULL) {
@@ -467,7 +467,7 @@ int functionJacC(DATA* data, threadData_t *threadData, double* jac){
 int functionJacD(DATA* data, threadData_t *threadData, double* jac){
 
   const int index = data->callback->INDEX_JAC_D;
-  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
+  JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
   unsigned int i,j,k;
   k = 0;
   if (jacobian->constantEqns != NULL) {
@@ -572,7 +572,7 @@ int linearize(DATA* data, threadData_t *threadData)
     if (data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sizeTmpVars > 0){
         /* Retrieve symbolic Jacobian */
         /* Determine Matrix A */
-        ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
+        JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
         if(!data->callback->initialAnalyticJacobianA(data, threadData, jacobian)){
             assertStreamPrint(threadData,0==functionJacA(data, threadData, matrixA),"Error, can not get Matrix A ");
         }

OMCompiler/SimulationRuntime/c/linearization/linearize.cpp

@@ -274,12 +274,12 @@ struct OpenModelicaGeneratedFunctionCallbacks {
   /*
   * These functions calculate specific jacobian column.
   */
-  analyticalJacobianColumn_func_ptr functionJacA_column;
-  analyticalJacobianColumn_func_ptr functionJacB_column;
-  analyticalJacobianColumn_func_ptr functionJacC_column;
-  analyticalJacobianColumn_func_ptr functionJacD_column;
-  analyticalJacobianColumn_func_ptr functionJacF_column;
-  analyticalJacobianColumn_func_ptr functionJacH_column;
+  jacobianColumn_func_ptr functionJacA_column;
+  jacobianColumn_func_ptr functionJacB_column;
+  jacobianColumn_func_ptr functionJacC_column;
+  jacobianColumn_func_ptr functionJacD_column;
+  jacobianColumn_func_ptr functionJacF_column;
+  jacobianColumn_func_ptr functionJacH_column;
   /*#endif*/
 
   const char *(*linear_model_frame)(void); /* printf format-string with holes for 6 strings */
@@ -373,14 +373,14 @@ struct OpenModelicaGeneratedFunctionCallbacks {
   * FMU continuous partial derivative functions
   */
   initialAnalyticalJacobian_func_ptr initialPartialFMIDER;
-  analyticalJacobianColumn_func_ptr functionJacFMIDER_column;
+  jacobianColumn_func_ptr functionJacFMIDER_column;
   const int INDEX_JAC_FMIDER;
 
   /*
   * FMU partial derivative functions for initilization DAE
   */
   initialAnalyticalJacobian_func_ptr initialPartialFMIDERINIT;
-  analyticalJacobianColumn_func_ptr functionJacFMIDERINIT_column;
+  jacobianColumn_func_ptr functionJacFMIDERINIT_column;
   const int INDEX_JAC_FMIDERINIT;
 };
 

OMCompiler/SimulationRuntime/c/openmodelica_func.h

@@ -811,7 +811,7 @@ void diffSynColoredOptimizerSystem(OptData *optData, modelica_real **J, const in
   int i,j,l,ii, ll;
 
   const int h_index = optData->s.indexABCD[index];
-  ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]);
+  JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]);
   const long double * scaldt = optData->bounds.scaldt[m];
   const unsigned int * const cC = jacobian->sparsePattern->colorCols;
   const unsigned int * const lindex  = jacobian->sparsePattern->leadindex;
@@ -877,7 +877,7 @@ void diffSynColoredOptimizerSystemF(OptData *optData, modelica_real **J){
     int i,j,l,ii, ll;
     const int index = 4;
     const int h_index = optData->s.indexABCD[index];
-    ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]);
+    JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]);
     const unsigned int * const cC = jacobian->sparsePattern->colorCols;
     const unsigned int * const lindex  = jacobian->sparsePattern->leadindex;
     const int nx = jacobian->sizeCols;

OMCompiler/SimulationRuntime/c/optimization/DataManagement/MoveData.c

@@ -48,13 +48,15 @@
  *                                  NULL if not available.
  * @param sparsePattern             Pointer to sparsity pattern of Jacobian.
  */
-void initAnalyticJacobian(ANALYTIC_JACOBIAN* jacobian, unsigned int sizeCols, unsigned int sizeRows, unsigned int sizeTmpVars, int (*constantEqns)(void* data, threadData_t *threadData, void* thisJacobian, void* parentJacobian), SPARSE_PATTERN* sparsePattern) {
+void initJacobian(JACOBIAN* jacobian, unsigned int sizeCols, unsigned int sizeRows, unsigned int sizeTmpVars, jacobianColumn_func_ptr evalColumn, jacobianColumn_func_ptr constantEqns, SPARSE_PATTERN* sparsePattern)
+{
   jacobian->sizeCols = sizeCols;
   jacobian->sizeRows = sizeRows;
   jacobian->sizeTmpVars = sizeTmpVars;
   jacobian->seedVars = (modelica_real*) calloc(sizeCols, sizeof(modelica_real));
   jacobian->resultVars = (modelica_real*) calloc(sizeRows, sizeof(modelica_real));
   jacobian->tmpVars = (modelica_real*) calloc(sizeTmpVars, sizeof(modelica_real));
+  jacobian->evalColumn = evalColumn;
   jacobian->constantEqns = constantEqns;
   jacobian->sparsePattern = sparsePattern;
   jacobian->availability = JACOBIAN_UNKNOWN;
@@ -67,16 +69,18 @@ void initAnalyticJacobian(ANALYTIC_JACOBIAN* jacobian, unsigned int sizeCols, un
  * Sparsity pattern is not copied, only the pointer to it.
  *
  * @param source                  Jacobian that should be copied.
- * @return ANALYTIC_JACOBIAN*     Copy of source.
+ * @return JACOBIAN*              Copy of source.
  */
-ANALYTIC_JACOBIAN* copyAnalyticJacobian(ANALYTIC_JACOBIAN* source) {
-  ANALYTIC_JACOBIAN* jacobian = (ANALYTIC_JACOBIAN*) malloc(sizeof(ANALYTIC_JACOBIAN));
-  initAnalyticJacobian(jacobian,
-                       source->sizeCols,
-                       source->sizeRows,
-                       source->sizeTmpVars,
-                       source->constantEqns,
-                       source->sparsePattern);
+JACOBIAN* copyJacobian(JACOBIAN* source)
+{
+  JACOBIAN* jacobian = (JACOBIAN*) malloc(sizeof(JACOBIAN));
+  initJacobian(jacobian,
+    source->sizeCols,
+    source->sizeRows,
+    source->sizeTmpVars,
+    source->evalColumn,
+    source->constantEqns,
+    source->sparsePattern);
 
   return jacobian;
 }
@@ -88,7 +92,8 @@ ANALYTIC_JACOBIAN* copyAnalyticJacobian(ANALYTIC_JACOBIAN* source) {
  *
  * @param jac   Pointer to Jacobian.
  */
-void freeAnalyticJacobian(ANALYTIC_JACOBIAN *jac) {
+void freeJacobian(JACOBIAN *jac)
+{
   if (jac == NULL) {
     return;
   }
@@ -99,6 +104,53 @@ void freeAnalyticJacobian(ANALYTIC_JACOBIAN *jac) {
   free(jac->sparsePattern); jac->sparsePattern = NULL;
 }
 
+/*! \fn evalJacobian
+ *
+ *  compute entries of Jacobian
+ *  uses coloring (from sparsity pattern) if available
+ *
+ *  \param [ref] [data]
+ *  \param [ref] [threadData]
+ *  \param [ref] [jacobian]   Pointer to Jacobian
+ *  \param [out] [jac]        Contains jacobian values on exit, column-major
+ */
+void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac)
+{
+  int i,j,k,l,ii;
+  const SPARSE_PATTERN* sp = jacobian->sparsePattern;
+
+  memset(jac, 0.0, (jacobian->sizeRows) * (jacobian->sizeCols) * sizeof(modelica_real));
+
+  /* evaluate constant equations of Jacobian */
+  if (jacobian->constantEqns != NULL) {
+    jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
+  }
+
+  /* evaluate Jacobian */
+  for (i = 0; i < sp->maxColors; i++) {
+    /* activate seed variable for the corresponding color */
+    for (j = 0; j < jacobian->sizeCols; j++)
+      if (sp->colorCols[j]-1 == i)
+        jacobian->seedVars[j] = 1.0;
+
+    /* evaluate Jacobian column */
+    jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
+
+    for (j = 0; j < jacobian->sizeCols; j++) {
+      if (sp->colorCols[j]-1 == i) {
+        for (ii = sp->leadindex[j]; ii < sp->leadindex[j+1]; ii++) {
+          l = sp->index[ii];
+          k = j*jacobian->sizeRows + l;
+          jac[k] = jacobian->resultVars[l]; //* solverData->xScaling[j];
+        }
+        /* de-activate seed variable for the corresponding color */
+        jacobian->seedVars[j] = 0.0;
+      }
+    }
+  }
+}
+
+
 /**
  * @brief Allocate memory for sparsity pattern.
  *
@@ -108,7 +160,8 @@ void freeAnalyticJacobian(ANALYTIC_JACOBIAN *jac) {
  * @param maxColors           Maximum number of colors of Matrix.
  * @return SPARSE_PATTERN*    Pointer ot allocated sparsity pattern of Matrix.
  */
-SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int numberOfNonZeros, unsigned int maxColors) {
+SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int numberOfNonZeros, unsigned int maxColors)
+{
   SPARSE_PATTERN* sparsePattern = (SPARSE_PATTERN*) malloc(sizeof(SPARSE_PATTERN));
   sparsePattern->leadindex = (unsigned int*) malloc((n_leadIndex+1)*sizeof(unsigned int));
   sparsePattern->index = (unsigned int*) malloc(numberOfNonZeros*sizeof(unsigned int));
@@ -125,7 +178,8 @@ SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int number
  *
  * @param spp   Pointer to sparsity pattern
  */
-void freeSparsePattern(SPARSE_PATTERN *spp) {
+void freeSparsePattern(SPARSE_PATTERN *spp)
+{
   if (spp != NULL) {
     free(spp->index); spp->index = NULL;
     free(spp->colorCols); spp->colorCols = NULL;
@@ -141,7 +195,8 @@ void freeSparsePattern(SPARSE_PATTERN *spp) {
  * @param filename    String for the filename.
  * @return FILE*      Pointer to sparsity pattern stream.
  */
-FILE * openSparsePatternFile(DATA* data, threadData_t *threadData, const char* filename) {
+FILE * openSparsePatternFile(DATA* data, threadData_t *threadData, const char* filename)
+{
   FILE* pFile;
   const char* fullPath = NULL;
 
@@ -168,7 +223,8 @@ FILE * openSparsePatternFile(DATA* data, threadData_t *threadData, const char* f
  * @param color         Current color index.
  * @param length        Number of columns in color `color`.
  */
-void readSparsePatternColor(threadData_t* threadData, FILE * pFile, unsigned int* colorCols, unsigned int color, unsigned int length, unsigned int maxIndex) {
+void readSparsePatternColor(threadData_t* threadData, FILE * pFile, unsigned int* colorCols, unsigned int color, unsigned int length, unsigned int maxIndex)
+{
   unsigned int i, index;
   size_t count;
 
@@ -192,7 +248,8 @@ void readSparsePatternColor(threadData_t* threadData, FILE * pFile, unsigned int
  * @param flagValue               Flag value of FLAG_JACOBIAN. Can be NULL.
  * @return JACOBIAN_METHOD   Returns jacobian method that is availble.
  */
-JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILITY availability, const char* flagValue){
+JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILITY availability, const char* flagValue)
+{
   JACOBIAN_METHOD jacobianMethod = JAC_UNKNOWN;
   assertStreamPrint(threadData, availability != JACOBIAN_UNKNOWN, "Jacobian availablity status is unknown.");
 
@@ -205,7 +262,7 @@ JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILIT
       }
     }
     // Error case
-    if(jacobianMethod == JAC_UNKNOWN){
+    if (jacobianMethod == JAC_UNKNOWN) {
       errorStreamPrint(OMC_LOG_STDOUT, 0, "Unknown value `%s` for flag `-jacobian`", flagValue);
       infoStreamPrint(OMC_LOG_STDOUT, 1, "Available options are");
       for (int method=1; method < JAC_MAX; method++) {
@@ -271,7 +328,8 @@ JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILIT
   return jacobianMethod;
 }
 
-void freeNonlinearPattern(NONLINEAR_PATTERN *nlp) {
+void freeNonlinearPattern(NONLINEAR_PATTERN *nlp)
+{
   if (nlp != NULL) {
     free(nlp->indexVar); nlp->indexVar = NULL;
     free(nlp->indexEqn); nlp->indexEqn = NULL;
@@ -280,19 +338,20 @@ void freeNonlinearPattern(NONLINEAR_PATTERN *nlp) {
   }
 }
 
-unsigned int* getNonlinearPatternCol(NONLINEAR_PATTERN *nlp, int var_idx){
+unsigned int* getNonlinearPatternCol(NONLINEAR_PATTERN *nlp, int var_idx)
+{
   unsigned int idx_start = nlp->indexVar[var_idx];
   unsigned int idx_stop;
-  if (var_idx == nlp->numberOfVars){
+  if (var_idx == nlp->numberOfVars) {
     idx_stop = nlp->numberOfNonlinear;
-  }else{
+  } else {
     idx_stop = nlp->indexVar[var_idx + 1];
   }
 
   unsigned int* col = (unsigned int*) malloc((idx_stop - idx_start + 1)*sizeof(unsigned int));
 
   int index = 0;
-  for(int i = idx_start; i < idx_stop + 1; i++){
+  for (int i = idx_start; i < idx_stop + 1; i++) {
     col[index] = nlp->columns[i];
     index++;
   }
@@ -305,12 +364,13 @@ unsigned int* getNonlinearPatternCol(NONLINEAR_PATTERN *nlp, int var_idx){
   return col;
 }
 
-unsigned int* getNonlinearPatternRow(NONLINEAR_PATTERN *nlp, int eqn_idx){
+unsigned int* getNonlinearPatternRow(NONLINEAR_PATTERN *nlp, int eqn_idx)
+{
   unsigned int idx_start = nlp->indexEqn[eqn_idx];
   unsigned int idx_stop;
-  if (eqn_idx == nlp->numberOfEqns){
+  if (eqn_idx == nlp->numberOfEqns) {
     idx_stop = nlp->numberOfNonlinear;
-  }else{
+  } else {
     idx_stop = nlp->indexEqn[eqn_idx + 1];
   }
   //printf("   eqn_idx   = %d\n", eqn_idx);
@@ -319,7 +379,7 @@ unsigned int* getNonlinearPatternRow(NONLINEAR_PATTERN *nlp, int eqn_idx){
   unsigned int* row = (unsigned int*) malloc((idx_stop - idx_start + 1)*sizeof(unsigned int));
 
   int index = 0;
-  for(int i = idx_start; i < idx_stop + 1; i++){
+  for (int i = idx_start; i < idx_stop + 1; i++) {
     row[index] = nlp->rows[i];
     //printf("      row[index] = row[%d] = %d\n", index, row[index]);
     index++;

OMCompiler/SimulationRuntime/c/simulation/jacobian_util.c

@@ -40,9 +40,11 @@
 extern "C" {
 #endif
 
-void initAnalyticJacobian(ANALYTIC_JACOBIAN* jacobian, unsigned int sizeCols, unsigned int sizeRows, unsigned int sizeTmpVars, int (*constantEqns)(void* data, threadData_t *threadData, void* thisJacobian, void* parentJacobian), SPARSE_PATTERN* sparsePattern);
-ANALYTIC_JACOBIAN* copyAnalyticJacobian(ANALYTIC_JACOBIAN* source);
-void freeAnalyticJacobian(ANALYTIC_JACOBIAN* jac);
+void initJacobian(JACOBIAN* jacobian, unsigned int sizeCols, unsigned int sizeRows, unsigned int sizeTmpVars, jacobianColumn_func_ptr evalColumn, jacobianColumn_func_ptr constantEqns, SPARSE_PATTERN* sparsePattern);
+JACOBIAN* copyJacobian(JACOBIAN* source);
+void freeJacobian(JACOBIAN* jac);
+
+void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac);
 
 SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int numberOfNonZeros, unsigned int maxColors);
 void freeSparsePattern(SPARSE_PATTERN *spp);

OMCompiler/SimulationRuntime/c/simulation/jacobian_util.h

@@ -526,7 +526,7 @@ int cvode_solver_initial(DATA *data, threadData_t *threadData, SOLVER_INFO *solv
   int flag;
   int i;
   double *abstol_tmp;
-  ANALYTIC_JACOBIAN *jacobian;
+  JACOBIAN *jacobian;
 
   /* Log cvode_initial */
   infoStreamPrint(OMC_LOG_SOLVER_V, 0, "### Start initialize of CVODE solver ###");