- add error message for division by zero to initial equations

- use c comments
- better error message vor pivot.c

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@14718 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

Author: Jens Frenkel

Compiler/BackEnd/SimCodeUtil.mo

@@ -1570,6 +1570,7 @@ algorithm
       startValueEquations = List.map(startValueEquations,addDivExpErrorMsgtoSimEqSystem);
       parameterEquations = List.map(parameterEquations,addDivExpErrorMsgtoSimEqSystem);
       removedEquations = List.map(removedEquations,addDivExpErrorMsgtoSimEqSystem);
+      initialEquations = List.map(initialEquations,addDivExpErrorMsgtoSimEqSystem);
 
       // generate jacobian or linear model matrices
       LinearMatrices = createJacobianLinearCode(symJacs, modelInfo, uniqueEqIndex);

SimulationRuntime/c/math-support/pivot.c

@@ -33,24 +33,24 @@
 /* 
 find the maximum element below (and including) line/row start
 */
-void maxsearch( double *A, modelica_integer start, modelica_integer n_rows, modelica_integer n_cols, modelica_integer *rowInd, modelica_integer *colInd, modelica_integer *maxrow, modelica_integer *maxcol, double *maxabsval)
+int maxsearch( double *A, modelica_integer start, modelica_integer n_rows, modelica_integer n_cols, modelica_integer *rowInd, modelica_integer *colInd, modelica_integer *maxrow, modelica_integer *maxcol, double *maxabsval)
 {
-  // temporary variables
+  /* temporary variables */
   modelica_integer row;
   modelica_integer col;
 
-  // Initialization
+  /* Initialization */
   modelica_integer mrow = -1;
   modelica_integer mcol = -1;
   double mabsval = 0.0;
 
-  // go through all rows and columns
+  /* go through all rows and columns */
   for(row=start; row<n_rows; row++)
   {
     for(col=start; col<n_cols; col++)
     {
       double tmp = fabs(get_pivot_matrix_elt(A,row,col));
-      // Compare element to current maximum
+      /* Compare element to current maximum */
       if (tmp > mabsval)
       {
         mrow = row;
@@ -60,11 +60,10 @@ void maxsearch( double *A, modelica_integer start, modelica_integer n_rows, mode
     }
   }
 
-  // assert that the matrix is not identical to zero
-  assert(mrow >= 0);
-  assert(mcol >= 0);
+  /* assert that the matrix is not identical to zero */
+  if ((mrow < 0) || (mcol < 0)) return -1;
 
-  // return result
+  // return result */
   *maxrow = mrow;
   *maxcol = mcol;
   *maxabsval = mabsval;
@@ -78,52 +77,54 @@ rowInd and colInd are vectors of length nrwos and n_cols respectively.
 They hold the old (and new) pivoting information, such that
   A_pivoted[i,j] = A[rowInd[i], colInd[j]]
 */
-void pivot( double *A, modelica_integer n_rows, modelica_integer n_cols, modelica_integer *rowInd, modelica_integer *colInd )
+int pivot( double *A, modelica_integer n_rows, modelica_integer n_cols, modelica_integer *rowInd, modelica_integer *colInd )
 {
-  // parameter, determines how much larger an element should be before rows and columns are interchanged
+  /* parameter, determines how much larger an element should be before rows and columns are interchanged */
   const double fac = 1.125; // approved by dymola ;)
 
-  // temporary variables
+  /* temporary variables */
   modelica_integer row;
   modelica_integer i,j;
   modelica_integer maxrow;
   modelica_integer maxcol;
   double maxabsval;
   double pivot;
 
-  // go over all pivot elements
+  /* go over all pivot elements */
   for(row=0; row<min(n_rows,n_cols); row++)
   {
-    // get current pivot
+    /* get current pivot */
     pivot = fabs(get_pivot_matrix_elt(A,row,row));
 
-    // find the maximum element in matrix
-    // result is stored in maxrow, maxcol and maxabsval
-    maxsearch(A, row, n_rows, n_cols, rowInd, colInd, &maxrow, &maxcol, &maxabsval);
+    /* find the maximum element in matrix
+       result is stored in maxrow, maxcol and maxabsval */
+    if (maxsearch(A, row, n_rows, n_cols, rowInd, colInd, &maxrow, &maxcol, &maxabsval) != 0) return -1;
 
-    // compare max element and pivot (scaled by fac)
+
+    /* compare max element and pivot (scaled by fac) */
     if (maxabsval > (fac*pivot))
     {
-      // row interchange
+      /* row interchange */
       swap(rowInd[row], rowInd[maxrow]);
-      // column interchange
+      /* column interchange */
       swap(colInd[row], colInd[maxcol]);
     }
 
-    // get pivot (without abs, may have changed because of row/column interchange
+    /* get pivot (without abs, may have changed because of row/column interchange */
     pivot = get_pivot_matrix_elt(A,row,row);
+    /* internal error, pivot element should never be zero if maxsearch succeeded */
     assert(pivot != 0);
 
-    // perform one step of Gaussian Elimination
+    /* perform one step of Gaussian Elimination */
     for(i=row+1;i<n_rows;i++)
     {
       double leader = get_pivot_matrix_elt(A,i,row);
       if (leader != 0.0)
       {
         double scale = -leader/pivot;
-        // set leader to zero
+        /* set leader to zero */
         set_pivot_matrix_elt(A,i,row, 0.0);
-        // subtract scaled equation from pivot row from current row
+        /* subtract scaled equation from pivot row from current row */
         for(j=row+1;j<n_cols;j++)
         {
           double t1 = get_pivot_matrix_elt(A,i,j);
@@ -134,5 +135,7 @@ void pivot( double *A, modelica_integer n_rows, modelica_integer n_cols, modelic
       }
     }
   }
+  /* all fine */
+  return 0;
 }
 

SimulationRuntime/c/simulation/solver/stateset.c

@@ -31,6 +31,7 @@
 
 #include "stateset.h"
 #include "matrix.h"
+#include "omc_error.h"
 
 #include <memory.h>
 
@@ -212,6 +213,7 @@ int comparePivot(modelica_integer *oldPivot, modelica_integer *newPivot, modelic
 int stateSelection(DATA *data)
 {
   long i=0;
+  long j=0;
   int res=0;
   /* go troug all state sets*/
   for (i=0;i<data->modelData.nStateSets;i++)
@@ -222,7 +224,21 @@ int stateSelection(DATA *data)
     getAnalyticalJacobianSet(data,i);
     /* call pivoting function to select the states */
     memcpy(oldColPivot,set->colPivot,set->nCandidates*sizeof(modelica_integer));
-    pivot(set->J,set->nDummyStates,set->nCandidates,set->rowPivot,set->colPivot);
+    if (pivot(set->J,set->nDummyStates,set->nCandidates,set->rowPivot,set->colPivot) != 0)
+    {
+      /* error, report the matrix and the time */
+      for(i=0;  i < data->simulationInfo.analyticJacobians[set->jacobianIndex].sizeRows;i++)
+      {
+        for(j=0;  j < data->simulationInfo.analyticJacobians[set->jacobianIndex].sizeCols;j++)
+          printf("% .5e ",set->J[i*data->simulationInfo.analyticJacobians[set->jacobianIndex].sizeCols+j]);
+        printf("\n");
+      }
+      for (i=0;i<set->nCandidates;i++)
+      {
+        printf("%s\n",set->statescandidates[i]->name);
+      }
+      THROW1("Error, singular Jacobian for dynamic state selection at time %f\nUse -lv LOG_DSSJAC to get the Jacobian\n",data->localData[0]->timeValue);
+    }
     /* if we have a new set throw event for reinitialization 
        and set the A matrix for set.x=A*(states) */
     res = comparePivot(oldColPivot,set->colPivot,set->nCandidates,set->nDummyStates,set->nStates,set->A,set->states,set->statescandidates,data);

SimulationRuntime/c/simulation/solver/stateset.h

@@ -41,7 +41,7 @@ void initializeStateSetJacobians(DATA *data);
 int stateSelection(DATA *data);
 
 /* do pivoting */
-extern void pivot( double *A, modelica_integer n_rows, modelica_integer n_cols, modelica_integer *rowInd, modelica_integer *colInd );
+extern int pivot( double *A, modelica_integer n_rows, modelica_integer n_cols, modelica_integer *rowInd, modelica_integer *colInd );
 
 
 #endif /* OMC_STATESET_H */