nonlinear homotopy solver with lapack

 - added -nls-LS option with options [lapack|totalpivot]
 - make lapack as default, since it's faster
 - in the next step klu will added

Author: Willi Braun

SimulationRuntime/c/simulation/simulation_runtime.cpp

@@ -301,6 +301,28 @@ static int getNewtonStrategy()
   return NEWTON_NONE;
 }
 
+static int getNlsLSSolver()
+{
+  int i;
+  const char *cflags = omc_flagValue[FLAG_NLS_LS];
+  const string *method = cflags ? new string(cflags) : NULL;
+
+  if(!method)
+    return NLS_LS_LAPACK; /* default method */
+
+  for(i=1; i<NLS_LS_MAX; ++i)
+    if(*method == NLS_LS_METHOD[i])
+      return i;
+
+  warningStreamPrint(LOG_STDOUT, 1, "unrecognized option -nls=%s, current options are:", method->c_str());
+  for(i=1; i<NLS_LS_MAX; ++i)
+    warningStreamPrint(LOG_STDOUT, 0, "%-18s [%s]", NLS_LS_METHOD[i], NLS_LS_METHOD_DESC[i]);
+  messageClose(LOG_STDOUT);
+  throwStreamPrint(NULL,"see last warning");
+
+  return NLS_LS_UNKNOWN;
+}
+
 static double getFlagReal(enum _FLAG flag, double res)
 {
   const char *flagStr = omc_flagValue[flag];
@@ -792,6 +814,7 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data, threadData_t *thr
   data->simulationInfo->lssMethod = getlinearSparseSolverMethod();
   data->simulationInfo->newtonStrategy = getNewtonStrategy();
   data->simulationInfo->nlsCsvInfomation = omc_flag[FLAG_NLS_INFO];
+  data->simulationInfo->nlsLinearSolver = getNlsLSSolver();
 
   if(omc_flag[FLAG_LSS_MAX_DENSITY]) {
     linearSparseSolverMaxDensity = atof(omc_flagValue[FLAG_LSS_MAX_DENSITY]);

SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c

@@ -49,6 +49,16 @@
 #include "nonlinearSolverHomotopy.h"
 #include "nonlinearSolverHybrd.h"
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
+
+#ifdef __cplusplus
+}
+#endif
+
 /*! \typedef DATA_HOMOTOPY
  * define memory structure for nonlinear system solver
  *  \author bbachmann
@@ -953,6 +963,7 @@ static int wrapper_fvec_homotopy_fixpoint_der(DATA_HOMOTOPY* solverData, double*
   return 0;
 }
 
+
 /*! \fn getIndicesOfPivotElement for calculating pivot element
  *
  *  \author bbachmann
@@ -976,6 +987,7 @@ static int wrapper_fvec_homotopy_fixpoint_der(DATA_HOMOTOPY* solverData, double*
   }
 }
 
+
 /*! \fn solveSystemWithTotalPivotSearch for solution of overdetermined linear system
  *  used for the homotopy solver, for calculating the direction
  *  used for the newton solver, for calculating the Newton step
@@ -994,6 +1006,7 @@ int solveSystemWithTotalPivotSearch(int n, double* x, double* A, int* indRow, in
    double *res;
 
    debugMatrixDouble(LOG_NLS_JAC,"Linear System Matrix [Jac res]:",A, n, m);
+   debugVectorDouble(LOG_NLS_JAC,"vector b:", A+n*n, n);
 
    /* assume full rank of matrix [n x (n+1)] */
    *rank = n;
@@ -1077,20 +1090,85 @@ int solveSystemWithTotalPivotSearch(int n, double* x, double* A, int* indRow, in
     *pos=indCol[n];
   }
 
+  return 0;
+}
+
+
+/*! \fn linearSolverWrapper
+ */
+int linearSolverWrapper(int n, double* x, double* A, int* indRow, int* indCol, int *pos, int *rank, int method)
+{
+  /* First try to use lapack and if it fails then
+   * use solveSystemWithTotalPivotSearch */
+  int returnValue = -1;
+  int solverinfo;
+  int nrhs = 1;
+  int lda = n;
+
+  debugMatrixDouble(LOG_NLS_JAC,"Linear System Matrix [Jac res]:", A, n, n+1);
+  debugVectorDouble(LOG_NLS_JAC,"vector b:", x, n);
+
+  switch(method){
+    case (1): /* NLS_LS_TOTALPIVOT */
+
+      solverinfo = solveSystemWithTotalPivotSearch(n, x, A, indRow, indCol, pos, rank);
+      /* in case of failing */
+      if (solverinfo != 0)
+      {
+        /* debug information */
+        debugString(LOG_NLS_V, "Linear total pivot solver failed!!!");
+        debugString(LOG_NLS_V, "******************************************************");
+      }
+      else
+      {
+        returnValue = 0;
+      }
+      break;
+    case 2: /* NLS_LS_LAPACK */
+      /* Solve system with lapack */
+      dgesv_((int*) &n,
+          (int*) &nrhs,
+          A,
+          (int*) &lda,
+          indRow,
+          x,
+          (int*) &n,
+          &solverinfo);
+
+      debugMatrixDouble(LOG_NLS_JAC,"Linear system matrix [Jac res] after decomposition:", A, n, n+1);
+      /* in case of failing */
+      if (solverinfo != 0)
+      {
+        /* debug information */
+        debugString(LOG_NLS_V, "Linear lapack solver failed!!!");
+        debugString(LOG_NLS_V, "******************************************************");
+      }
+      else
+      {
+        vecScalarMult(n, x, -1, x);
+        returnValue = 0;
+      }
+      break;
+    default:
+      warningStreamPrint(LOG_STDOUT, 0, "Non-Linear solver try to run with a unknown linear solver.");
+  }
+
   /* Debugging error of linear system */
   if(ACTIVE_STREAM(LOG_NLS_JAC))
   {
-    res = (double*) calloc(n,sizeof(double));
-    debugVectorDouble(LOG_NLS_JAC,"solution:", x, m);
-    matVecMult(n, m, A, x, res);
+    double* res = (double*) calloc(n,sizeof(double));
+    debugVectorDouble(LOG_NLS_JAC,"solution:", x, n);
+    matVecMult(n, n, A, x, res);
     debugVectorDouble(LOG_NLS_JAC,"test solution:", res, n);
     debugDouble(LOG_NLS_JAC,"error of linear system = ", vec2Norm(n, res));
     free(res);
     messageClose(LOG_NLS_JAC);
   }
 
-  return 0;
+  return returnValue;
 }
+
+
 /*! \fn solve system with damped Newton-Raphson
  *
  *  \author bbachmann
@@ -1114,6 +1192,7 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
   int assert = 1;
   threadData_t *threadData = solverData->threadData;
   NONLINEAR_SYSTEM_DATA* nonlinsys = &(solverData->data->simulationInfo->nonlinearSystemData[solverData->data->simulationInfo->currentNonlinearSystemIndex]);
+  int linearSolverMethod = solverData->data->simulationInfo->nlsLinearSolver;
 
   /* debug information */
   debugString(LOG_NLS_V, "******************************************************");
@@ -1135,7 +1214,7 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
     debugInt(LOG_NLS_V, "Iteration:", numberOfIterations);
 
     /* solve jacobian and function value (both stored in hJac, last column is fvec), side effects: jacobian matrix is changed */
-    if ((numberOfIterations>1) && (solveSystemWithTotalPivotSearch(solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank) != 0))
+    if ((numberOfIterations>1) && linearSolverWrapper(solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank, linearSolverMethod) != 0)
     {
       /* report solver abortion */
       solverData->info=-1;
@@ -1392,6 +1471,7 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
     matVecMultAbsBB(solverData->n, solverData->fJac, solverData->ones, solverData->resScaling);
     debugVectorDouble(LOG_NLS_JAC, "residuum scaling:", solverData->resScaling, solverData->n);
     scaleMatrixRows(solverData->n, solverData->m, solverData->fJac);
+    vecCopy(n, solverData->fJac + n*n, solverData->dy0);
   }
   return 0;
 }

SimulationRuntime/c/simulation_data.h

@@ -547,6 +547,7 @@ typedef struct SIMULATION_INFO
   int nlsMethod;                       /* nonlinear solver */
   int newtonStrategy;                  /* newton damping strategy solver */
   int nlsCsvInfomation;                /* = 1 csv files with detailed nonlinear solver process are generated */
+  int nlsLinearSolver;                 /* nls linear solver setting =1 totalpivot, =2 lapack, =3=klu */
 
   /* current context evaluation, set by dassl and used for extrapolation
    * of next non-linear guess */

SimulationRuntime/c/util/simulation_options.c

@@ -83,6 +83,7 @@ const char *FLAG_NAME[FLAG_MAX+1] = {
   /* FLAG_NEWTON_STRATEGY */       "newton",
   /* FLAG_NLS */                   "nls",
   /* FLAG_NLS_INFO */              "nlsInfo",
+  /* FLAG_NLS_LS */                "nlsLS",
   /* FLAG_NOEMIT */                "noemit",
   /* FLAG_NOEQUIDISTANT_GRID */    "noEquidistantTimeGrid",
   /* FLAG_NOEQUIDISTANT_OUT_FREQ*/ "noEquidistantOutputFrequency",
@@ -161,6 +162,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = {
   /* FLAG_NEWTON_STRATEGY */       "value specifies the damping strategy for the newton solver",
   /* FLAG_NLS */                   "value specifies the nonlinear solver",
   /* FLAG_NLS_INFO */              "outputs detailed information about solving process of non-linear systems into csv files.",
+  /* FLAG_NLS_LS */                "value specifies the linear solver used by the non-linear solver\m nlsLS=[totalpivot|lapack|klu]",
   /* FLAG_NOEMIT */                "do not emit any results to the result file",
   /* FLAG_NOEQUIDISTANT_GRID */    "stores results not in equidistant time grid as given by stepSize or numberOfIntervals, instead the variable step size of dassl is used.",
   /* FLAG_NOEQUIDISTANT_OUT_FREQ*/ "value controls the output frequency in noEquidistantTimeGrid mode",
@@ -328,6 +330,8 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
   "  * mixed",
   /* FLAG_NLS_INFO */
   "  Outputs detailed information about solving process of non-linear systems into csv files.",
+  /* FLAG_NLS_LS */
+  "value specifies the linear solver used by the non-linear solver\m nlsLS=[totalpivot|lapack|klu]",
   /* FLAG_NOEMIT */
   "  Do not emit any results to the result file.",
   /* FLAG_NOEQUIDISTANT_GRID */
@@ -443,6 +447,7 @@ const int FLAG_TYPE[FLAG_MAX] = {
   /* FLAG_NEWTON_STRATEGY */       FLAG_TYPE_OPTION,
   /* FLAG_NLS */                   FLAG_TYPE_OPTION,
   /* FLAG_NLS_INFO */              FLAG_TYPE_FLAG,
+  /* FLAG_NLS_LS */                FLAG_TYPE_OPTION,
   /* FLAG_NOEMIT */                FLAG_TYPE_FLAG,
   /* FLAG_NOEQUIDISTANT_GRID*/     FLAG_TYPE_FLAG,
   /* FLAG_NOEQUIDISTANT_OUT_FREQ*/ FLAG_TYPE_OPTION,
@@ -676,3 +681,22 @@ const char *IDA_LS_METHOD_DESC[IDA_LS_MAX+1] = {
 
   "IDA_LS_MAX"
 };
+
+const char *NLS_LS_METHOD[NLS_LS_MAX+1] = {
+  "unknown",
+
+  "totalpivot",
+  "lapack",
+
+  "NLS_LS_MAX"
+};
+
+const char *NLS_LS_METHOD_DESC[NLS_LS_MAX+1] = {
+  "unknown",
+
+  "internal total pivot implementation. Solve in some case even under-determined systems.",
+  "use external lapack implementation.",
+
+  "NLS_LS_MAX"
+};
+

SimulationRuntime/c/util/simulation_options.h

@@ -91,6 +91,7 @@ enum _FLAG
   FLAG_NEWTON_STRATEGY,
   FLAG_NLS,
   FLAG_NLS_INFO,
+  FLAG_NLS_LS,
   FLAG_NOEMIT,
   FLAG_NOEQUIDISTANT_GRID,
   FLAG_NOEQUIDISTANT_OUT_FREQ,
@@ -271,8 +272,18 @@ enum IDA_LS
 extern const char *IDA_LS_METHOD[IDA_LS_MAX+1];
 extern const char *IDA_LS_METHOD_DESC[IDA_LS_MAX+1];
 
+enum NLS_LS
+{
+  NLS_LS_UNKNOWN = 0,
+
+  NLS_LS_TOTALPIVOT,
+  NLS_LS_LAPACK,
 
+  NLS_LS_MAX
+};
 
+extern const char *NLS_LS_METHOD[NLS_LS_MAX+1];
+extern const char *NLS_LS_METHOD_DESC[NLS_LS_MAX+1];
 
 
 #if defined(__cplusplus)