diff --git a/OMCompiler/Compiler/Template/CodegenC.tpl b/OMCompiler/Compiler/Template/CodegenC.tpl index 7b24a79f99c..1c2353af597 100644 --- a/OMCompiler/Compiler/Template/CodegenC.tpl +++ b/OMCompiler/Compiler/Template/CodegenC.tpl @@ -117,6 +117,7 @@ end translateModel; #include "simulation/simulation_info_json.h" #include "simulation/simulation_runtime.h" #include "util/omc_error.h" + #include "util/parallel_helper.h" #include "simulation/solver/model_help.h" #include "simulation/solver/delay.h" #include "simulation/solver/linearSystem.h" @@ -1957,7 +1958,6 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].setA = setLinearMatrixA<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].setb = setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; - linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; >> case SOME(__) then let size = listLength(ls.vars) @@ -1976,10 +1976,9 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac%>; linearSystemData[<%ls.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac%>; linearSystemData[<%ls.indexLinearSystem%>].jacobianIndex = <%jacIndex%>; - linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%ls.index%>; - linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].setA = NULL; //setLinearMatrixA<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; - linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; >> else error(sourceInfo(), ' No jacobian create for linear system <%ls.index%>.') @@ -2006,7 +2005,6 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].setA = setLinearMatrixA<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].setb = setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; - linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; assertStreamPrint(NULL, nLinearSystems > <%at.indexLinearSystem%>, "Internal Error: nLinearSystems mismatch!"); linearSystemData[<%at.indexLinearSystem%>].equationIndex = <%at.index%>; @@ -2017,7 +2015,6 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%at.indexLinearSystem%>].setA = setLinearMatrixA<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].setb = setLinearVectorb<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%at.index%>; - linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; linearSystemData[<%at.indexLinearSystem%>].checkConstraints = checkConstraints<%at.index%>; >> case SOME(__) then @@ -2043,9 +2040,8 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac%>; linearSystemData[<%ls.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac%>; linearSystemData[<%ls.indexLinearSystem%>].jacobianIndex = <%jacIndex%>; - linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; - linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%ls.index%>; - linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].setA = NULL; //setLinearMatrixA<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; assertStreamPrint(NULL, nLinearSystems > <%at.indexLinearSystem%>, "Internal Error: indexlinearSystem mismatch!"); @@ -2058,9 +2054,8 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%at.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac2%>; linearSystemData[<%at.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac2%>; linearSystemData[<%at.indexLinearSystem%>].jacobianIndex = <%jacIndex2%>; - linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; - linearSystemData[<%at.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%at.index%>; - linearSystemData[<%at.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%at.index%>; + linearSystemData[<%at.indexLinearSystem%>].setA = NULL; //setLinearMatrixA<%at.index%>; + linearSystemData[<%at.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].checkConstraints = checkConstraints<%at.index%>; >> @@ -2131,7 +2126,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]); const int equationIndexes[2] = {1,<%ls.index%>}; <% if ls.partOfJac then - 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;' + 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian;' %> ANALYTIC_JACOBIAN* jacobian = NULL; <%varDeclsRes%> @@ -2182,7 +2177,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%> <%varDecls%> <%MatrixA%> } @@ -2192,7 +2187,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%> <%varDecls2%> <%vectorb%> } @@ -2270,7 +2265,9 @@ template functionSetupLinearSystemsTemp(list linearSystems, String DATA *data = (DATA*) ((void**)dataIn[0]); threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]); const int equationIndexes[2] = {1,<%ls.index%>}; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%> + <% if ls.partOfJac then + 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian;' + %> ANALYTIC_JACOBIAN* jacobian = NULL; <%varDeclsRes%> <% if profileAll() then 'SIM_PROF_TICK_EQ(<%ls.index%>);' %> @@ -2299,7 +2296,9 @@ template functionSetupLinearSystemsTemp(list linearSystems, String DATA *data = (DATA*) ((void**)dataIn[0]); threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]); const int equationIndexes[2] = {1,<%at.index%>}; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%> + <% if ls.partOfJac then + 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian;' + %> ANALYTIC_JACOBIAN* jacobian = NULL; <%varDeclsRes2%> <% if profileAll() then 'SIM_PROF_TICK_EQ(<%at.index%>);' %> @@ -2373,7 +2372,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%> <%varDecls%> <%MatrixA%> } @@ -2383,7 +2382,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%> <%varDecls2%> <%vectorb%> } @@ -2403,7 +2402,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%at.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%> <%varDecls3%> <%MatrixA2%> } @@ -2413,7 +2412,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%at.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; - <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parDynamicData[omc_get_thread_num()].parentJacobian;'%> <%varDecls4%> <%vectorb2%> } @@ -5585,7 +5584,7 @@ case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = at) th } <% if profileSome() then 'SIM_PROF_TICK_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex);' %> <% if ls.partOfJac then - 'data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = jacobian;' + 'data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parDynamicData[omc_get_thread_num()].parentJacobian = jacobian;' %> retValue = solve_linear_system(data, threadData, <%ls.indexLinearSystem%>, &aux_x[0]); @@ -5599,6 +5598,7 @@ case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = at) th /* write solution */ <%ls.vars |> SIMVAR(__) hasindex i0 => '<%contextCref(name, context, auxFunctions)%> = aux_x[<%i0%>];' ;separator="\n"%> <% if profileSome() then 'SIM_PROF_ACC_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex);' %> + <%returnval%> >> end equationLinear; diff --git a/OMCompiler/SimulationRuntime/c/Makefile.common b/OMCompiler/SimulationRuntime/c/Makefile.common index 8c9793cbe34..eb956bebf61 100644 --- a/OMCompiler/SimulationRuntime/c/Makefile.common +++ b/OMCompiler/SimulationRuntime/c/Makefile.common @@ -103,6 +103,7 @@ RUNTIMEUTIL_HEADERS = \ ./util/generic_array.h \ ./util/index_spec.h \ ./util/integer_array.h \ +./util/jacobian_util.h \ ./util/java_interface.h \ ./util/modelica.h \ ./util/modelica_string.h \ @@ -111,6 +112,7 @@ RUNTIMEUTIL_HEADERS = \ ./util/omc_mmap.h \ ./util/omc_msvc.h \ ./util/omc_spinlock.h \ +./util/parallel_helper.h \ ./util/read_matlab4.c \ ./util/read_matlab4.h \ ./util/read_csv.c \ diff --git a/OMCompiler/SimulationRuntime/c/Makefile.objs b/OMCompiler/SimulationRuntime/c/Makefile.objs index feb5d8d544f..98b5cc6eaa7 100644 --- a/OMCompiler/SimulationRuntime/c/Makefile.objs +++ b/OMCompiler/SimulationRuntime/c/Makefile.objs @@ -28,7 +28,7 @@ else UTIL_OBJS_NO_FMI= endif -UTIL_OBJS_MINIMAL=base_array$(OBJ_EXT) boolean_array$(OBJ_EXT) omc_error$(OBJ_EXT) omc_file$(OBJ_EXT) division$(OBJ_EXT) generic_array$(OBJ_EXT) index_spec$(OBJ_EXT) integer_array$(OBJ_EXT) list$(OBJ_EXT) modelica_string$(OBJ_EXT) real_array$(OBJ_EXT) ringbuffer$(OBJ_EXT) string_array$(OBJ_EXT) utility$(OBJ_EXT) varinfo$(OBJ_EXT) ModelicaUtilities$(OBJ_EXT) omc_msvc$(OBJ_EXT) simulation_options$(OBJ_EXT) rational$(OBJ_EXT) modelica_string_lit$(OBJ_EXT) omc_init$(OBJ_EXT) omc_mmap$(OBJ_EXT) $(UTIL_OBJS_NO_FMI) +UTIL_OBJS_MINIMAL=base_array$(OBJ_EXT) boolean_array$(OBJ_EXT) omc_error$(OBJ_EXT) omc_file$(OBJ_EXT) division$(OBJ_EXT) generic_array$(OBJ_EXT) index_spec$(OBJ_EXT) integer_array$(OBJ_EXT) list$(OBJ_EXT) modelica_string$(OBJ_EXT) real_array$(OBJ_EXT) ringbuffer$(OBJ_EXT) string_array$(OBJ_EXT) utility$(OBJ_EXT) varinfo$(OBJ_EXT) ModelicaUtilities$(OBJ_EXT) omc_msvc$(OBJ_EXT) parallel_helper$(OBJ_EXT) simulation_options$(OBJ_EXT) rational$(OBJ_EXT) modelica_string_lit$(OBJ_EXT) omc_init$(OBJ_EXT) omc_mmap$(OBJ_EXT) jacobian_util$(OBJ_EXT) $(UTIL_OBJS_NO_FMI) UTIL_HFILES_MINIMAL=base_array.h boolean_array.h division.h generic_array.h omc_error.h omc_file.h index_spec.h integer_array.h list.h modelica.h modelica_string.h read_write.h real_array.h ringbuffer.h rtclock.h string_array.h utility.h varinfo.h simulation_options.h omc_mmap.h modelica_string_lit.h omc_init.h ifeq ($(OMC_MINIMAL_RUNTIME),) diff --git a/OMCompiler/SimulationRuntime/c/simulation/simulation_runtime.cpp b/OMCompiler/SimulationRuntime/c/simulation/simulation_runtime.cpp index f078d3853ad..3aa1073ec08 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/simulation_runtime.cpp +++ b/OMCompiler/SimulationRuntime/c/simulation/simulation_runtime.cpp @@ -86,6 +86,7 @@ #include "simulation/solver/initialization/initialization.h" #include "simulation/solver/dae_mode.h" #include "dataReconciliation/dataReconciliation.h" +#include "util/parallel_helper.h" #ifdef _OMC_QSS_LIB #include "solver_qss/solver_qss.h" @@ -940,6 +941,36 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data, threadData_t *thr data->simulationInfo->minStepSize = 4.0 * DBL_EPSILON * fmax(fabs(data->simulationInfo->startTime),fabs(data->simulationInfo->stopTime)); rt_accumulate(SIM_TIMER_INIT_XML); + /* Set the maximum number of threads prior to any allocation w.r.t. + * linear systems and Jacobians in order to avoid memory leaks. + */ +#ifdef USE_PARJAC + int num_threads = omc_get_max_threads(); + if (omc_flag[FLAG_JACOBIAN_THREADS]) { + int num_threads_tmp = atoi(omc_flagValue[FLAG_JACOBIAN_THREADS]); + infoStreamPrint(LOG_STDOUT, 0, + "Number of threads passed via -jacobianThreads: %d", + num_threads_tmp); + if (0 >= num_threads_tmp) { + warningStreamPrint(LOG_STDOUT, 0, + "Number of desired OpenMP threads for parallel Jacobian evaluation is <= 0."); + warningStreamPrint(LOG_STDOUT, 0, "Use omp_get_max_threads()."); + } else { + num_threads = num_threads_tmp; + } + } + omp_set_num_threads(num_threads); + + infoStreamPrint(LOG_STDOUT, 0, + "Number of OpenMP threads for parallel Jacobian evaluation: %d", + omc_get_max_threads()); +#else + if (omc_flag[FLAG_JACOBIAN_THREADS]) { + warningStreamPrint(LOG_STDOUT, 0, + "Simulation flag jacobianThreads not available. Make sure you have configured omc with \"--enable-parjac\" and build with a compiler supporting OpenMP."); + } +#endif + /* initialize static data of mixed/linear/non-linear system solvers */ initializeMixedSystems(data, threadData); initializeLinearSystems(data, threadData); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.c b/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.c index 1ee3815f8fc..4a0da41c4b8 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.c @@ -27,6 +27,11 @@ * CONDITIONS OF OSMC-PL. * */ +#ifdef USE_PARJAC + #include + #define GC_THREADS + #include +#endif #include #include @@ -36,6 +41,7 @@ #include "simulation_data.h" #include "util/omc_error.h" +#include "util/parallel_helper.h" #include "gc/omc_gc.h" #include "simulation/options.h" @@ -193,11 +199,11 @@ int dassl_initial(DATA* data, threadData_t *threadData, dasslData->idid = 0; dasslData->ysave = (double*) malloc(N*sizeof(double)); - dasslData->ypsave = (double*) malloc(N*sizeof(double)); dasslData->delta_hh = (double*) malloc(N*sizeof(double)); dasslData->newdelta = (double*) malloc(N*sizeof(double)); dasslData->stateDer = (double*) calloc(N, sizeof(double)); dasslData->states = (double*) malloc(N*sizeof(double)); + dasslData->allocatedParMem = 0; /* false */ data->simulationInfo->currentContext = CONTEXT_ALGEBRAIC; @@ -368,9 +374,17 @@ int dassl_initial(DATA* data, threadData_t *threadData, case COLOREDSYMJAC: data->simulationInfo->jacobianEvals = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern->maxColors; dasslData->jacobianFunction = jacA_symColored; +#ifdef USE_PARJAC + allocateThreadLocalJacobians(data, &(dasslData->jacColumns)); + dasslData->allocatedParMem = 1; /* true */ +#endif break; case SYMJAC: dasslData->jacobianFunction = jacA_sym; +#ifdef USE_PARJAC + allocateThreadLocalJacobians(data, &(dasslData->jacColumns)); + dasslData->allocatedParMem = 1; /* true */ +#endif break; case NUMJAC: dasslData->jacobianFunction = jacA_num; @@ -432,17 +446,25 @@ int dassl_deinitial(DASSL_DATA *dasslData) /* free work arrays for DASSL */ free(dasslData->rwork); free(dasslData->iwork); + free(dasslData->jroot); free(dasslData->rpar); free(dasslData->ipar); free(dasslData->atol); free(dasslData->rtol); free(dasslData->info); - free(dasslData->jroot); free(dasslData->ysave); free(dasslData->delta_hh); free(dasslData->newdelta); - free(dasslData->states); free(dasslData->stateDer); + free(dasslData->states); + + +#ifdef USE_PARJAC + if (dasslData->allocatedParMem) { + freeAnalyticalJacobian(&(dasslData->jacColumns)); + dasslData->allocatedParMem = 0; + } +#endif free(dasslData); @@ -920,8 +942,13 @@ int jacA_symColored(double *t, double *y, double *yprime, double *delta, threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2]; DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1]; +#ifdef USE_PARJAC + ANALYTIC_JACOBIAN* jac = (dasslData->jacColumns); +#else const int index = data->callback->INDEX_JAC_A; ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]); +#endif + unsigned int columns = jac->sizeCols; unsigned int rows = jac->sizeRows; unsigned int sizeTmpVars = jac->sizeTmpVars; @@ -939,6 +966,7 @@ int jacA_symColored(double *t, double *y, double *yprime, double *delta, * * * This function calculates symbolically the jacobian matrix. + * Can calculate the jacobian in parallel. */ int jacA_sym(double *t, double *y, double *yprime, double *delta, double *matrixA, double *cj, double *h, double *wt, double *rpar, @@ -951,13 +979,35 @@ int jacA_sym(double *t, double *y, double *yprime, double *delta, threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2]; const int index = data->callback->INDEX_JAC_A; - ANALYTIC_JACOBIAN* t_jac = &(data->simulationInfo->analyticJacobians[index]); - unsigned int columns = t_jac->sizeCols; - unsigned int rows = t_jac->sizeRows; - unsigned int sizeTmpVars = t_jac->sizeTmpVars; + ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]); + unsigned int columns = jac->sizeCols; + unsigned int rows = jac->sizeRows; + unsigned int sizeTmpVars = jac->sizeTmpVars; + unsigned int i; - unsigned int i,j; +#ifdef USE_PARJAC + GC_allow_register_threads(); +#endif + +#pragma omp parallel default(none) firstprivate(columns, rows, sizeTmpVars) shared(i, matrixA, data, threadData, dasslData) +{ +#ifdef USE_PARJAC + /* Register omp-thread in GC */ + if(!GC_thread_is_registered()) { + struct GC_stack_base sb; + memset (&sb, 0, sizeof(sb)); + GC_get_stack_base(&sb); + GC_register_my_thread (&sb); + } + // Use thread local analytic Jacobians + ANALYTIC_JACOBIAN* t_jac = &(dasslData->jacColumns[omc_get_thread_num()]); + //printf("index= %d, t_jac->sizeCols= %d, t_jac->sizeRows = %d, t_jac->sizeTmpVars = %d \n",index, t_jac->sizeCols , t_jac->sizeRows, t_jac->sizeTmpVars); +#else + ANALYTIC_JACOBIAN* t_jac = jac; +#endif + unsigned int j; +#pragma omp for schedule(runtime) for(i=0; i < columns; i++) { t_jac->seedVars[i] = 1.0; @@ -967,7 +1017,8 @@ int jacA_sym(double *t, double *y, double *yprime, double *delta, matrixA[i*columns+j] = t_jac->resultVars[j]; t_jac->seedVars[i] = 0.0; - } + } // for loop +} // omp parallel TRACE_POP return 0; @@ -992,7 +1043,6 @@ int jacA_num(double *t, double *y, double *yprime, double *delta, double delta_h = numericalDifferentiationDeltaXsolver; double delta_hh,delta_hhh, deltaInv; double ysave; - double ypsave; int ires; int i,j; @@ -1055,7 +1105,6 @@ int jacA_numColored(double *t, double *y, double *yprime, double *delta, int ires; double* delta_hh = dasslData->delta_hh; double* ysave = dasslData->ysave; - double* ypsave = dasslData->ypsave; unsigned int i,j,l,k,ii; diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.h b/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.h index 37b07fa2f47..daf3ef1c5e9 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/dassl.h @@ -67,7 +67,6 @@ typedef struct DASSL_DATA{ /* variables used in jacobian calculation */ double *ysave; - double *ypsave; double *delta_hh; double *newdelta; double *stateDer; @@ -81,6 +80,10 @@ typedef struct DASSL_DATA{ double *rpar, int* ipar); void* zeroCrossingFunction; +#ifdef USE_PARJAC + ANALYTIC_JACOBIAN* jacColumns; /* thread local analytic jacobians */ +#endif + int allocatedParMem; /* indicated if parallel memory was allocated, 0=false, 1=true*/ } DASSL_DATA; /* main dassl function to make a step */ diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.c b/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.c index aa67d61ef6a..685809f659d 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.c @@ -31,6 +31,12 @@ /*! \file ida_solver.c */ +#ifdef USE_PARJAC + #include + #define GC_THREADS + #include +#endif + #include #include @@ -40,6 +46,7 @@ #include "simulation_data.h" #include "util/omc_error.h" +#include "util/parallel_helper.h" #include "gc/omc_gc.h" #include "simulation/options.h" @@ -79,14 +86,15 @@ #include #include +/* Function prototypes */ static int callDenseJacobian(long int Neq, double tt, double cj, - N_Vector yy, N_Vector yp, N_Vector rr, - DlsMat Jac, void *user_data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + N_Vector yy, N_Vector yp, N_Vector rr, + DlsMat Jac, void *user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); -static int callSparseJacobian(realtype tt, realtype cj, - N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, void *user_data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); +static int callSparseJacobian(realtype currentTime, realtype cj, + N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, void *user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); static int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData); int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* userData); @@ -98,9 +106,12 @@ static int idaReScaleData(IDA_SOLVER *idaData); static int idaScaleVector(N_Vector vec, double* factors, unsigned int size); static int idaReScaleVector(N_Vector vec, double* factors, unsigned int size); +int ida_event_update(DATA* data, threadData_t *threadData); + +/* Static variables */ static IDA_SOLVER *idaDataGlobal; static int initializedSolver = 0; -int ida_event_update(DATA* data, threadData_t *threadData); + int checkIDAflag(int flag) { @@ -121,7 +132,7 @@ int checkIDAflag(int flag) } void errOutputIDA(int error_code, const char *module, const char *function, - char *msg, void *userData) + char *msg, void *userData) { TRACE_PUSH DATA* data = (DATA*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->data); @@ -132,10 +143,10 @@ void errOutputIDA(int error_code, const char *module, const char *function, TRACE_POP } -/* initial main ida data */ -int -ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, IDA_SOLVER *idaData){ - +/* initialize main ida data */ +int ida_solver_initial(DATA* data, threadData_t *threadData, + SOLVER_INFO* solverInfo, IDA_SOLVER *idaData) +{ TRACE_PUSH int flag; @@ -269,11 +280,11 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo /* daeMode: set nominal values for algebraic variables */ if (idaData->daeMode) { - getAlgebraicDAEVarNominals(data, idaData->yScale + data->modelData->nStates); - for(i=data->modelData->nStates; i < idaData->N; ++i) - { - idaData->ypScale[i] = 1.0; - } + getAlgebraicDAEVarNominals(data, idaData->yScale + data->modelData->nStates); + for(i=data->modelData->nStates; i < idaData->N; ++i) + { + idaData->ypScale[i] = 1.0; + } } infoStreamPrint(LOG_SOLVER_V, 1, "The scale factors for all ida states: "); for(i=0; i < idaData->N; ++i) @@ -497,12 +508,18 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo throwStreamPrint(threadData, "##IDA## Setup of linear solver KLU failed!"); } + idaData->allocatedParMem = 0; /* false */ + switch (idaData->jacobianMethod){ case SYMJAC: case NUMJAC: case COLOREDSYMJAC: case COLOREDNUMJAC: flag = IDASlsSetSparseJacFn(idaData->ida_mem, callSparseJacobian); +#ifdef USE_PARJAC + allocateThreadLocalJacobians(data, &(idaData->jacColumns)); + idaData->allocatedParMem = 1; /* true */ +#endif break; default: throwStreamPrint(threadData,"For the klu solver jacobian calculation method has to be %s or %s", JACOBIAN_METHOD[COLOREDSYMJAC], JACOBIAN_METHOD[COLOREDNUMJAC]); @@ -518,6 +535,10 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo case NUMJAC: /* set jacobian function */ flag = IDADlsSetDenseJacFn(idaData->ida_mem, callDenseJacobian); +#ifdef USE_PARJAC + allocateThreadLocalJacobians(data, &(idaData->jacColumns)); + idaData->allocatedParMem = 1; /* true */ +#endif break; case INTERNALNUMJAC: break; @@ -677,9 +698,9 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo return 0; } -/* deinitial ida data */ -int -ida_solver_deinitial(IDA_SOLVER *idaData){ +/* deinitialize ida data */ +int ida_solver_deinitial(IDA_SOLVER *idaData) +{ TRACE_PUSH free(idaData->simData); @@ -706,14 +727,20 @@ ida_solver_deinitial(IDA_SOLVER *idaData){ N_VDestroy_Serial(idaData->errwgt); N_VDestroy_Serial(idaData->newdelta); +#ifdef USE_PARJAC + if (idaData->allocatedParMem) { + freeAnalyticalJacobian(&(idaData->jacColumns)); + idaData->allocatedParMem = 0; + } +#endif + IDAFree(&idaData->ida_mem); TRACE_POP return 0; } -int -ida_event_update(DATA* data, threadData_t *threadData) +int ida_event_update(DATA* data, threadData_t *threadData) { IDA_SOLVER *idaData = idaDataGlobal; int flag; @@ -798,8 +825,7 @@ ida_event_update(DATA* data, threadData_t *threadData) } /* main ida function to make a step */ -int -ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo) +int ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo) { TRACE_PUSH double tout = 0; @@ -1353,7 +1379,8 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* * into a dense DlsMat matrix */ static -int jacColoredNumericalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat Jac, double cj, void *userData) +int jacColoredNumericalDense(double currentTime, N_Vector yy, N_Vector yp, + N_Vector rr, DlsMat Jac, double cj, void *userData) { TRACE_PUSH IDA_SOLVER* idaData = (IDA_SOLVER*)userData; @@ -1396,7 +1423,7 @@ int jacColoredNumericalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, D sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern; } - setContext(data, &tt, CONTEXT_JACOBIAN); + setContext(data, ¤tTime, CONTEXT_JACOBIAN); for(i = 0; i < sparsePattern->maxColors; i++) { @@ -1420,7 +1447,7 @@ int jacColoredNumericalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, D } } - (*idaData->residualFunction)(tt, yy, yp, idaData->newdelta, userData); + (*idaData->residualFunction)(currentTime, yy, yp, idaData->newdelta, userData); increaseJacContext(data); @@ -1450,41 +1477,62 @@ int jacColoredNumericalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, D } /* - * function calculates the Jacobian matrix symbolical - * with considering also the coloring and pass it in a - * dense DlsMat matrix. + * function calculates the Jacobian matrix symbolically with considering also + * the coloring and pass it in a dense DlsMat matrix. */ -static -int jacColoredSymbolicalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat Jac, double cj, void *userData) +static int jacColoredSymbolicalDense(double currentTime, N_Vector yy, + N_Vector yp, N_Vector rr, DlsMat Jac, + double cj, void *userData) { TRACE_PUSH IDA_SOLVER* idaData = (IDA_SOLVER*)userData; DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data); threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)idaData->simData)->threadData); void* ida_mem = idaData->ida_mem; + long int N = idaData->N; const int index = data->callback->INDEX_JAC_A; unsigned int i,ii,j, nth; - ANALYTIC_JACOBIAN* jacData = &(data->simulationInfo->analyticJacobians[index]); SPARSE_PATTERN* sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern; - long int N = idaData->N; /* prepare variables */ double *states = N_VGetArrayPointer(yy); double *yprime = N_VGetArrayPointer(yp); - setContext(data, &tt, CONTEXT_SYM_JACOBIAN); + setContext(data, ¤tTime, CONTEXT_SYM_JACOBIAN); /* Reuse jacobian matrix in KLU solver */ +#ifdef USE_PARJAC + GC_allow_register_threads(); +#endif + +#pragma omp parallel default(none) firstprivate(N) shared(i, sparsePattern, idaData, data, threadData, Jac) private(ii, j, nth) +{ +#ifdef USE_PARJAC + /* Register omp-thread in GC */ + if(!GC_thread_is_registered()) { + struct GC_stack_base sb; + memset (&sb, 0, sizeof(sb)); + GC_get_stack_base(&sb); + GC_register_my_thread (&sb); + } + // ToDo Use always a thread local analyticJacobians (replace simulationInfo->analyticaJacobians) + // These are not the Jacobians of the linear systems! (SimulationInfo->linearSystemData[idx].jacobian) + ANALYTIC_JACOBIAN* t_jac = &(idaData->jacColumns[omc_get_thread_num()]); +#else + ANALYTIC_JACOBIAN* t_jac = &(data->simulationInfo->analyticJacobians[index]); +#endif + +#pragma omp for for(i = 0; i < sparsePattern->maxColors; i++) { for(ii=0; ii < N; ii++) { if(sparsePattern->colorCols[ii]-1 == i) { - jacData->seedVars[ii] = 1; + t_jac->seedVars[ii] = 1; } } - data->callback->functionJacA_column(data, threadData, jacData, NULL); + data->callback->functionJacA_column(data, threadData, t_jac, NULL); increaseJacContext(data); for(ii = 0; ii < N; ii++) @@ -1495,8 +1543,8 @@ int jacColoredSymbolicalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, while(nth < sparsePattern->leadindex[ii+1]) { j = sparsePattern->index[nth]; - infoStreamPrint(LOG_JAC, 0, "### symbolical jacobian at [%d,%d] = %f ###", j, ii, jacData->resultVars[j]); - DENSE_ELEM(Jac, j, ii) = jacData->resultVars[j]; + infoStreamPrint(LOG_JAC, 0, "### symbolical jacobian at [%d,%d] = %f ###", j, ii, t_jac->resultVars[j]); + DENSE_ELEM(Jac, j, ii) = t_jac->resultVars[j]; nth++; }; } @@ -1504,9 +1552,11 @@ int jacColoredSymbolicalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, for(ii=0; ii < idaData->N; ii++) { - jacData->seedVars[ii] = 0; + t_jac->seedVars[ii] = 0; } - } + } // for column +} // omp parallel + unsetContext(data); TRACE_POP @@ -1598,8 +1648,7 @@ static void finishSparseColPtr(SlsMat mat, int nnz) * numerical method finite differences with coloring * into a sparse SlsMat matrix */ -static -int jacoColoredNumericalSparse(double tt, N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, double cj, void *userData) +static int jacoColoredNumericalSparse(double currentTime, N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, double cj, void *userData) { TRACE_PUSH IDA_SOLVER* idaData = (IDA_SOLVER*)userData; @@ -1650,7 +1699,7 @@ int jacoColoredNumericalSparse(double tt, N_Vector yy, N_Vector yp, N_Vector rr, /* it's needed to clear the matrix */ SlsSetToZero(Jac); - setContext(data, &tt, CONTEXT_JACOBIAN); + setContext(data, ¤tTime, CONTEXT_JACOBIAN); /* rescale idaData->y and idaData->yp * the evaluation of the residual function @@ -1684,7 +1733,7 @@ int jacoColoredNumericalSparse(double tt, N_Vector yy, N_Vector yp, N_Vector rr, } } idaData->disableScaling = 1; - (*idaData->residualFunction)(tt, yy, yp, idaData->newdelta, userData); + (*idaData->residualFunction)(currentTime, yy, yp, idaData->newdelta, userData); idaData->disableScaling = disBackup; increaseJacContext(data); @@ -1733,7 +1782,7 @@ int jacoColoredNumericalSparse(double tt, N_Vector yy, N_Vector yp, N_Vector rr, /* * This function calculates the jacobian matrix symbolically while exploiting coloring. */ -int jacColoredSymbolicalSparse(double tt, N_Vector yy, N_Vector yp, +int jacColoredSymbolicalSparse(double currentTime, N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, double cj, void *userData) { @@ -1747,15 +1796,20 @@ int jacColoredSymbolicalSparse(double tt, N_Vector yy, N_Vector yp, double *states = N_VGetArrayPointer(yy); double *yprime = N_VGetArrayPointer(yp); +#ifdef USE_PARJAC + ANALYTIC_JACOBIAN* jacData = (idaData->jacColumns); +#else ANALYTIC_JACOBIAN* jacData = &(data->simulationInfo->analyticJacobians[index]); - SPARSE_PATTERN* sparsePattern = jacData->sparsePattern; +#endif unsigned int columns = jacData->sizeCols; unsigned int rows = jacData->sizeRows; + SPARSE_PATTERN* sparsePattern = jacData->sparsePattern; + int maxColors = sparsePattern->maxColors; /* it's needed to clear the matrix */ SlsSetToZero(Jac); - setContext(data, &tt, CONTEXT_SYM_JACOBIAN); + setContext(data, ¤tTime, CONTEXT_SYM_JACOBIAN); /* Reuse jacobian matrix in KLU solver */ genericColoredSymbolicJacobianEvaluation(rows, columns, sparsePattern, Jac, jacData, data, threadData, &setJacElementKluSparse); @@ -1770,10 +1824,10 @@ int jacColoredSymbolicalSparse(double tt, N_Vector yy, N_Vector yp, /* * Wrapper function to call numerical or symbolical jacobian matrix */ -static int callSparseJacobian(double tt, double cj, - N_Vector yy, N_Vector yp, N_Vector rr, - SlsMat Jac, void *user_data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) +static int callSparseJacobian(double currentTime, double cj, + N_Vector yy, N_Vector yp, N_Vector rr, + SlsMat Jac, void *user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { TRACE_PUSH @@ -1787,11 +1841,11 @@ static int callSparseJacobian(double tt, double cj, if (idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC) { - jacColoredSymbolicalSparse(tt, yy, yp, rr, Jac, cj, user_data); + jacColoredSymbolicalSparse(currentTime, yy, yp, rr, Jac, cj, user_data); } else if (idaData->jacobianMethod == COLOREDNUMJAC || idaData->jacobianMethod == NUMJAC) { - jacoColoredNumericalSparse(tt, yy, yp, rr, Jac, cj, user_data); + jacoColoredNumericalSparse(currentTime, yy, yp, rr, Jac, cj, user_data); } /* debug */ @@ -1821,8 +1875,7 @@ static int callSparseJacobian(double tt, double cj, return 0; } -static -int getScalingFactors(DATA* data, IDA_SOLVER *idaData, SlsMat inScaleMatrix) +static int getScalingFactors(DATA* data, IDA_SOLVER *idaData, SlsMat inScaleMatrix) { int i; @@ -1892,8 +1945,7 @@ int getScalingFactors(DATA* data, IDA_SOLVER *idaData, SlsMat inScaleMatrix) return 0; } -static -int idaScaleVector(N_Vector vec, double* factors, unsigned int size) +static int idaScaleVector(N_Vector vec, double* factors, unsigned int size) { int i; double *data = N_VGetArrayPointer(vec); @@ -1906,8 +1958,7 @@ int idaScaleVector(N_Vector vec, double* factors, unsigned int size) return 0; } -static -int idaReScaleVector(N_Vector vec, double* factors, unsigned int size) +static int idaReScaleVector(N_Vector vec, double* factors, unsigned int size) { int i; double *data = N_VGetArrayPointer(vec); @@ -1921,8 +1972,7 @@ int idaReScaleVector(N_Vector vec, double* factors, unsigned int size) return 0; } -static -int idaScaleData(IDA_SOLVER *idaData) +static int idaScaleData(IDA_SOLVER *idaData) { infoStreamPrint(LOG_SOLVER_V, 1, "Scale y"); idaScaleVector(idaData->y, idaData->yScale, idaData->N); @@ -1934,8 +1984,7 @@ int idaScaleData(IDA_SOLVER *idaData) return 0; } -static -int idaReScaleData(IDA_SOLVER *idaData) +static int idaReScaleData(IDA_SOLVER *idaData) { infoStreamPrint(LOG_SOLVER_V, 1, "Re-Scale y"); idaReScaleVector(idaData->y, idaData->yScale, idaData->N); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.h b/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.h index 29f72e0bffb..7122b65f4fd 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/ida_solver.h @@ -107,24 +107,26 @@ typedef struct IDA_SOLVER N_Vector* ySp; N_Vector* ySResult; +#ifdef USE_PARJAC + ANALYTIC_JACOBIAN* jacColumns; +#endif + int allocatedParMem; /* indicated if parallel memory was allocated, 0=false, 1=true*/ + }IDA_SOLVER; -/* initial main ida Data */ -int -ida_solver_initial(DATA* simData, threadData_t *threadData, SOLVER_INFO* solverInfo, IDA_SOLVER *idaData); +/* initialize main ida Data */ +int ida_solver_initial(DATA* data, threadData_t *threadData, + SOLVER_INFO* solverInfo, IDA_SOLVER *idaData); -/* deinitial main ida Data */ -int -ida_solver_deinitial(IDA_SOLVER *idaData); +/* deinitialize main ida Data */ +int ida_solver_deinitial(IDA_SOLVER *idaData); /* main ida function to make a step */ -int -ida_solver_step(DATA* simData, threadData_t *threadData, SOLVER_INFO* solverInfo); +int ida_solver_step(DATA* simData, threadData_t *threadData, SOLVER_INFO* solverInfo); /* event handing reinitialization function */ -int -ida_event_update(DATA* data, threadData_t *threadData); +int ida_event_update(DATA* data, threadData_t *threadData); -#endif +#endif /* #ifdef WITH_SUNDIALS */ -#endif +#endif /* #ifndef OMC_IDA_SOLVER_H*/ diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.c b/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.c index ed67a931a88..8e07296a62d 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.c @@ -31,8 +31,62 @@ /*! \file jacobian_symbolical.c */ +#ifdef USE_PARJAC + #define GC_THREADS + #include +#endif + #include "simulation/solver/jacobianSymbolical.h" +#ifdef USE_PARJAC +/** Allocate thread local Jacobians in case of OpenMP-parallel Jacobian computation. + * + * (symbolical only), used in IDA and Dassl. + */ +// ToDo AHEu: Make this usable without OpenMP and use it as default! +void allocateThreadLocalJacobians(DATA* data, ANALYTIC_JACOBIAN** jacColumns) +{ + int maxTh = omc_get_max_threads(); + *jacColumns = (ANALYTIC_JACOBIAN*) malloc(maxTh*sizeof(ANALYTIC_JACOBIAN)); + const int index = data->callback->INDEX_JAC_A; + ANALYTIC_JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[index]); + SPARSE_PATTERN* sparsePattern = data->simulationInfo->analyticJacobians[index].sparsePattern; + + unsigned int columns = jac->sizeCols; + unsigned int rows = jac->sizeRows; + unsigned int sizeTmpVars = jac->sizeTmpVars; + + unsigned int i; + +#ifdef USE_PARJAC + GC_allow_register_threads(); +#endif + +#pragma omp parallel default(none) firstprivate(maxTh, columns, rows, sizeTmpVars, index) shared(sparsePattern, jacColumns, i) + /* Benchmarks indicate that it is beneficial to initialize and malloc the jacColumns using a parallel for loop. */ + { + /* Register omp-thread in GC */ + if(!GC_thread_is_registered()) { + struct GC_stack_base sb; + memset (&sb, 0, sizeof(sb)); + GC_get_stack_base(&sb); + GC_register_my_thread (&sb); + } + +#pragma omp for schedule(runtime) + for (i = 0; i < maxTh; ++i) { + (*jacColumns)[i].sizeCols = columns; + (*jacColumns)[i].sizeRows = rows; + (*jacColumns)[i].sizeTmpVars = sizeTmpVars; + (*jacColumns)[i].tmpVars = (double*) calloc(sizeTmpVars, sizeof(double)); + (*jacColumns)[i].resultVars = (double*) calloc(rows, sizeof(double)); + (*jacColumns)[i].seedVars = (double*) calloc(columns, sizeof(double)); + (*jacColumns)[i].sparsePattern = sparsePattern; + } + } +} +#endif + /** * \brief Generic parallel computation of the colored Jacobian. @@ -55,9 +109,28 @@ void genericColoredSymbolicJacobianEvaluation(int rows, int columns, SPARSE_PATT threadData_t* threadData, void (*setJacElement)(int, int, int, double, void*, int)) { - ANALYTIC_JACOBIAN* t_jac = jacColumns; +#ifdef USE_PARJAC + GC_allow_register_threads(); +#endif + +#pragma omp parallel default(none) firstprivate(columns, rows) \ + shared(spp, matrixA, jacColumns, data, threadData, setJacElement) +{ +#ifdef USE_PARJAC + /* Register omp-thread in GC */ + if(!GC_thread_is_registered()) { + struct GC_stack_base sb; + memset (&sb, 0, sizeof(sb)); + GC_get_stack_base(&sb); + GC_register_my_thread (&sb); + } + // printf("My id = %d of max threads= %d\n", omc_get_thread_num(), omp_get_num_threads()); +#endif + ANALYTIC_JACOBIAN* t_jac = &(jacColumns[omc_get_thread_num()]); unsigned int i, j, currentIndex, nth; + +#pragma omp for for (i=0; i < spp->maxColors; i++) { /* Set seed vector for current color */ for (j=0; j < columns; j++) { @@ -85,4 +158,22 @@ void genericColoredSymbolicJacobianEvaluation(int rows, int columns, SPARSE_PATT t_jac->seedVars[j] = 0; } } +} // omp parallel +} + +#ifdef USE_PARJAC +/** Free ANALYTIC_JACOBIAN struct */ +void freeAnalyticalJacobian(ANALYTIC_JACOBIAN** jacColumns) +{ + int maxTh = omc_get_max_threads(); + unsigned int i; + + for (i = 0; i < maxTh; ++i) { + free((*jacColumns)[i].tmpVars); + free((*jacColumns)[i].resultVars); + free((*jacColumns)[i].seedVars); + } + + free(*jacColumns); } +#endif diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.h b/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.h index eb279b2ae51..8ca1e640b64 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/jacobianSymbolical.h @@ -35,6 +35,9 @@ #define OMC_JACOBIAN_SYMBOLICAL_H #include "../../simulation_data.h" +#include "util/parallel_helper.h" + +void allocateThreadLocalJacobians(DATA* data, ANALYTIC_JACOBIAN** jacColumns); void genericColoredSymbolicJacobianEvaluation(int rows, int columns, SPARSE_PATTERN* spp, void* matrixA, ANALYTIC_JACOBIAN* jacColumns, @@ -42,4 +45,6 @@ void genericColoredSymbolicJacobianEvaluation(int rows, int columns, SPARSE_PATT threadData_t* threadData, void (*setJacElement)(int, int, int, double, void*, int)); +void freeAnalyticalJacobian(ANALYTIC_JACOBIAN** jacColumns); + #endif diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c index 796b63fc859..2907a614038 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c @@ -41,6 +41,7 @@ #include "simulation_data.h" #include "simulation/simulation_info_json.h" #include "util/omc_error.h" +#include "util/parallel_helper.h" #include "omc_math.h" #include "util/varinfo.h" #include "model_help.h" @@ -48,7 +49,6 @@ #include "linearSystem.h" #include "linearSolverKlu.h" - static void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n); static void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n); @@ -123,8 +123,8 @@ static int getAnalyticalJacobian(DATA* data, threadData_t *threadData, LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]); const int index = systemData->jacobianIndex; - ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); - ANALYTIC_JACOBIAN* parentJacobian = systemData->parentJacobian; + ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian; + ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian; int nth = 0; int nnz = jacobian->sparsePattern->numberOfNoneZeros; @@ -179,7 +179,7 @@ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); - DATA_KLU* solverData = (DATA_KLU*)systemData->solverData[0]; + DATA_KLU* solverData = (DATA_KLU*)systemData->parDynamicData[omc_get_thread_num()].solverData[0]; _omc_scalar residualNorm = 0; int i, j, status = 0, success = 0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber}; @@ -217,7 +217,8 @@ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) /* calculate vector b (rhs) */ memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row); - residual_wrapper(solverData->work, systemData->b, dataAndThreadData, sysNumber); + + residual_wrapper(solverData->work, systemData->parDynamicData[omc_get_thread_num()].b, dataAndThreadData, sysNumber); } tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock)); systemData->jacobianTime += tmpJacEvalTime; @@ -240,11 +241,13 @@ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) messageClose(LOG_LS_V); for (i=0; in_row; i++) - infoStreamPrint(LOG_LS_V, 0, "b[%d] = %e", i, systemData->b[i]); + { + // ToDo Rework stream prints like this one to work in parallel regions + infoStreamPrint(LOG_LS_V, 0, "b[%d] = %e", i, systemData->parDynamicData[omc_get_thread_num()].b[i]); + } } rt_ext_tp_tick(&(solverData->timeClock)); - /* symbolic pre-ordering of A to reduce fill-in of L and U */ if (0 == solverData->numberSolving) { @@ -276,11 +279,11 @@ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) if (0 == solverData->common.status){ if (1 == systemData->method){ - if (klu_solve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->b, &solverData->common)){ + if (klu_solve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->parDynamicData[omc_get_thread_num()].b, &solverData->common)){ success = 1; } } else { - if (klu_tsolve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->b, &solverData->common)){ + if (klu_tsolve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->parDynamicData[omc_get_thread_num()].b, &solverData->common)){ success = 1; } } @@ -294,7 +297,7 @@ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) if (1 == systemData->method){ /* take the solution */ for(i = 0; i < solverData->n_row; ++i) - aux_x[i] += systemData->b[i]; + aux_x[i] += systemData->parDynamicData[omc_get_thread_num()].b[i]; /* update inner equations */ residual_wrapper(aux_x, solverData->work, dataAndThreadData, sysNumber); @@ -308,7 +311,7 @@ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) } } else { /* the solution is automatically in x */ - memcpy(aux_x, systemData->b, sizeof(double)*systemData->size); + memcpy(aux_x, systemData->parDynamicData[omc_get_thread_num()].b, sizeof(double)*systemData->size); } if (ACTIVE_STREAM(LOG_LS_V)) diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c index a73283f63ae..7f19f95de40 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c @@ -38,6 +38,7 @@ #include "../../simulation_data.h" #include "../simulation_info_json.h" #include "../../util/omc_error.h" +#include "../../util/parallel_helper.h" #include "omc_math.h" #include "../../util/varinfo.h" #include "model_help.h" @@ -45,6 +46,9 @@ #include "linearSystem.h" #include "linearSolverLapack.h" +#ifdef USE_PARJAC + #include +#endif extern int dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info); @@ -87,7 +91,7 @@ int freeLapackData(void **voiddata) _omc_destroyMatrix(data->A); free(data); - voiddata[0] = 0; + voiddata[0] = NULL; return 0; } @@ -108,13 +112,12 @@ int getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, double* ja LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[currentSys]); const int index = systemData->jacobianIndex; - ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); - ANALYTIC_JACOBIAN* parentJacobian = systemData->parentJacobian; + ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian; + ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian; memset(jac, 0, (systemData->size)*(systemData->size)*sizeof(double)); - for(i=0; i < jacobian->sparsePattern->maxColors; i++) - { + for(i=0; i < jacobian->sparsePattern->maxColors; i++) { /* activate seed variable for the corresponding color */ for(ii=0; ii < jacobian->sizeCols; ii++) if(jacobian->sparsePattern->colorCols[ii]-1 == i) @@ -122,13 +125,10 @@ int getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, double* ja ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian); - for(j = 0; j < jacobian->sizeCols; j++) - { - if(jacobian->seedVars[j] == 1) - { + for(j = 0; j < jacobian->sizeCols; j++) { + if(jacobian->seedVars[j] == 1) { ii = jacobian->sparsePattern->leadindex[j]; - while(ii < jacobian->sparsePattern->leadindex[j+1]) - { + while(ii < jacobian->sparsePattern->leadindex[j+1]) { l = jacobian->sparsePattern->index[ii]; k = j*jacobian->sizeRows + l; jac[k] = -jacobian->resultVars[l]; @@ -167,8 +167,8 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux void *dataAndThreadData[2] = {data, threadData}; int i, iflag = 1; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); - DATA_LAPACK* solverData = (DATA_LAPACK*)systemData->solverData[0]; + DATA_LAPACK* solverData = (DATA_LAPACK*) systemData->parDynamicData[omc_get_thread_num()].solverData[0]; int success = 1; /* We are given the number of the linear system. @@ -185,15 +185,16 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux /* set data */ _omc_setVectorData(solverData->x, aux_x); - _omc_setVectorData(solverData->b, systemData->b); - _omc_setMatrixData(solverData->A, systemData->A); + _omc_setVectorData(solverData->b, systemData->parDynamicData[omc_get_thread_num()].b); + _omc_setMatrixData(solverData->A, systemData->parDynamicData[omc_get_thread_num()].A); + // ToDo: Make time variables thread safe as this can be called in a parallel region rt_ext_tp_tick(&(solverData->timeClock)); if (0 == systemData->method) { - if (!reuseMatrixJac){ + if (!reuseMatrixJac) { /* reset matrix A */ - memset(systemData->A, 0, (systemData->size)*(systemData->size)*sizeof(double)); + memset(systemData->parDynamicData[omc_get_thread_num()].A, 0, (systemData->size)*(systemData->size)*sizeof(double)); /* update matrix A */ systemData->setA(data, threadData, systemData); } @@ -201,9 +202,9 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux /* update vector b (rhs) */ systemData->setb(data, threadData, systemData); } else { - if (!reuseMatrixJac){ + if (!reuseMatrixJac) { /* calculate jacobian -> matrix A*/ - if(systemData->jacobianIndex != -1){ + if(systemData->jacobianIndex != -1) { getAnalyticalJacobianLapack(data, threadData, solverData->A->data, sysNumber); } else { assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" ); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.h b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.h index 17abb511fd1..9d350dceff4 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.h @@ -52,7 +52,7 @@ typedef struct DATA_LAPACK } DATA_LAPACK; -int allocateLapackData(int size, void **data); +int allocateLapackData(int size, void** voiddata); int freeLapackData(void **data); int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.c b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.c index 37f41a4edf0..82d0668d7b2 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLis.c @@ -38,6 +38,7 @@ #include "simulation_data.h" #include "simulation/simulation_info_json.h" #include "util/omc_error.h" +#include "util/parallel_helper.h" #include "omc_math.h" #include "util/varinfo.h" #include "model_help.h" @@ -151,25 +152,20 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]); const int index = systemData->jacobianIndex; - ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); - ANALYTIC_JACOBIAN* parentJacobian = systemData->parentJacobian; + ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian; + ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian; int nth = 0; int nnz = jacobian->sparsePattern->numberOfNoneZeros; - for(i=0; i < jacobian->sizeRows; i++) - { + for(i=0; i < jacobian->sizeRows; i++) { jacobian->seedVars[i] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian); - for(j = 0; j < jacobian->sizeCols; j++) - { - if(jacobian->seedVars[j] == 1) - { + for(j = 0; j < jacobian->sizeCols; j++) { + if(jacobian->seedVars[j] == 1) { ii = jacobian->sparsePattern->leadindex[j]; - while(ii < jacobian->sparsePattern->leadindex[j+1]) - { + while(ii < jacobian->sparsePattern->leadindex[j+1]) { l = jacobian->sparsePattern->index[ii]; /*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", l, i, nth, -jacobian->resultVars[l]); */ systemData->setAElement(l, i, -jacobian->resultVars[l], nth, (void*) systemData, threadData); @@ -206,7 +202,8 @@ int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); - DATA_LIS* solverData = (DATA_LIS*)systemData->solverData[0]; + DATA_LIS* solverData = (DATA_LIS*)systemData->parDynamicData[omc_get_thread_num()].solverData[0]; + int i, ret, success = 1, ni, iflag = 1, n = systemData->size, eqSystemNumber = systemData->equationIndex; char *lis_returncode[] = {"LIS_SUCCESS", "LIS_ILL_OPTION", "LIS_BREAKDOWN", "LIS_OUT_OF_MEMORY", "LIS_MAXITER", "LIS_NOT_IMPLEMENTED", "LIS_ERR_FILE_IO"}; LIS_INT err; @@ -246,10 +243,11 @@ int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) /* calculate vector b (rhs) */ memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row); - wrapper_fvec_lis(solverData->work, systemData->b, dataAndThreadData, sysNumber); - /* set b vector */ - for(i=0; ib[i], solverData->b); + wrapper_fvec_lis(solverData->work, systemData->parDynamicData[omc_get_thread_num()].b, dataAndThreadData, sysNumber); + + /* set b vector */ + for(i=0; iparDynamicData[omc_get_thread_num()].b[i], solverData->b); } } tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock)); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c index 5ff392c23fa..13c0196dea3 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c @@ -38,6 +38,7 @@ #include "../../simulation_data.h" #include "../simulation_info_json.h" #include "../../util/omc_error.h" +#include "../../util/parallel_helper.h" #include "omc_math.h" #include "../../util/varinfo.h" #include "model_help.h" @@ -45,7 +46,6 @@ #include "linearSystem.h" #include "linearSolverTotalPivot.h" - void debugMatrixDoubleLS(int logName, char* matrixName, double* matrix, int n, int m) { if(ACTIVE_STREAM(logName)) @@ -304,7 +304,7 @@ int freeTotalPivotData(void** voiddata) free(data->indCol); free(voiddata[1]); - voiddata[1] = 0; + voiddata[1] = NULL; return 0; } @@ -326,8 +326,8 @@ int getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, double const int index = systemData->jacobianIndex; - ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); - ANALYTIC_JACOBIAN* parentJacobian = systemData->parentJacobian; + ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian; + ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian; memset(jac, 0, (systemData->size)*(systemData->size)*sizeof(double)); @@ -390,7 +390,8 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* void *dataAndThreadData[2] = {data, threadData}; int i, j; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); - DATA_TOTALPIVOT* solverData = (DATA_TOTALPIVOT*) systemData->solverData[1]; + DATA_TOTALPIVOT* solverData = (DATA_TOTALPIVOT*) systemData->parDynamicData[omc_get_thread_num()].solverData[1]; + int n = systemData->size, status; double fdeps = 1e-8; double xTol = 1e-8; @@ -413,19 +414,18 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* debugVectorDoubleLS(LOG_LS_V,"Old VALUES",aux_x,n); rt_ext_tp_tick(&(solverData->timeClock)); - if (0 == systemData->method){ + if (0 == systemData->method) { /* reset matrix A */ - vecConstLS(n*n, 0.0, systemData->A); + vecConstLS(n*n, 0.0, systemData->parDynamicData[omc_get_thread_num()].A); /* update matrix A -> first n columns of matrix Ab*/ systemData->setA(data, threadData, systemData); - vecCopyLS(n*n, systemData->A, solverData->Ab); + vecCopyLS(n*n, systemData->parDynamicData[omc_get_thread_num()].A, solverData->Ab); /* update vector b (rhs) -> -b is last column of matrix Ab*/ rt_ext_tp_tick(&(solverData->timeClock)); systemData->setb(data, threadData, systemData); - vecScalarMultLS(n, systemData->b, -1.0, solverData->Ab + n*n); - + vecScalarMultLS(n, systemData->parDynamicData[omc_get_thread_num()].b, -1.0, solverData->Ab + n*n); } else { /* calculate jacobian -> first n columns of matrix Ab*/ @@ -446,14 +446,18 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* status = solveSystemWithTotalPivotSearchLS(n, solverData->x, solverData->Ab, solverData->indRow, solverData->indCol, &rank); infoStreamPrint(LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock))); - if (status != 0) - { + if (status != 0) { + // ToDo Rework stream prints like this one to work in parallel regions +#ifdef USE_PARJAC + warningStreamPrint(LOG_STDOUT, 0, "Thread %u: Error solving linear system of equations (no. %d) at time %f.", omc_get_thread_num(), (int)systemData->equationIndex, data->localData[0]->timeValue); + success = 0; +#else warningStreamPrint(LOG_STDOUT, 0, "Error solving linear system of equations (no. %d) at time %f.", (int)systemData->equationIndex, data->localData[0]->timeValue); success = 0; +#endif } else { - debugVectorDoubleLS(LOG_LS_V, "SOLUTION:", solverData->x, n+1); - if (1 == systemData->method){ + if (1 == systemData->method) { /* add the solution to old solution vector*/ vecAddLS(n, aux_x, solverData->x, aux_x); wrapper_fvec_totalpivot(aux_x, solverData->b, dataAndThreadData, sysNumber); @@ -462,7 +466,7 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* vecCopyLS(n, solverData->x, aux_x); } - if (ACTIVE_STREAM(LOG_LS_V)){ + if (ACTIVE_STREAM(LOG_LS_V)) { if (1 == systemData->method) { infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm); } else { diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h index 35cce599986..d88e8a5f024 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h @@ -51,7 +51,7 @@ typedef struct DATA_TOTALPIVOT } DATA_TOTALPIVOT; -int allocateTotalPivotData(int size, void** data); +int allocateTotalPivotData(int size, void** voiddata); int freeTotalPivotData(void** data); int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c index 8684c4af489..8e7cfa8c49b 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c @@ -41,6 +41,7 @@ #include "simulation_data.h" #include "simulation/simulation_info_json.h" #include "util/omc_error.h" +#include "util/parallel_helper.h" #include "omc_math.h" #include "util/varinfo.h" #include "model_help.h" @@ -48,7 +49,6 @@ #include "linearSystem.h" #include "linearSolverUmfpack.h" - void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n); void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n); int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x); @@ -138,8 +138,8 @@ int getAnalyticalJacobianUmfPack(DATA* data, threadData_t *threadData, int sysNu LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]); const int index = systemData->jacobianIndex; - ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); - ANALYTIC_JACOBIAN* parentJacobian = systemData->parentJacobian; + ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian; + ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian; int nth = 0; int nnz = jacobian->sparsePattern->numberOfNoneZeros; @@ -197,7 +197,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); - DATA_UMFPACK* solverData = (DATA_UMFPACK*)systemData->solverData[0]; + DATA_UMFPACK* solverData = (DATA_UMFPACK*)systemData->parDynamicData[omc_get_thread_num()].solverData[0]; _omc_scalar residualNorm = 0; int i, j, status = UMFPACK_OK, success = 0, ni=0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber}; @@ -236,7 +236,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) /* calculate vector b (rhs) */ memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row); - wrapper_fvec_umfpack(solverData->work, systemData->b, dataAndThreadData, sysNumber); + wrapper_fvec_umfpack(solverData->work, systemData->parDynamicData[omc_get_thread_num()].b, dataAndThreadData, sysNumber); } tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock)); systemData->jacobianTime += tmpJacEvalTime; @@ -258,14 +258,15 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) } messageClose(LOG_LS_V); - for (i=0; in_row; i++) - infoStreamPrint(LOG_LS_V, 0, "b[%d] = %e", i, systemData->b[i]); + for (i=0; in_row; i++) { + // ToDo Rework stream prints like this one to work in parallel regions + infoStreamPrint(LOG_LS_V, 0, "b[%d] = %e", i, systemData->parDynamicData[omc_get_thread_num()].b[i]); + } } rt_ext_tp_tick(&(solverData->timeClock)); /* symbolic pre-ordering of A to reduce fill-in of L and U */ - if (0 == solverData->numberSolving) - { + if (0 == solverData->numberSolving) { status = umfpack_di_symbolic(solverData->n_col, solverData->n_row, solverData->Ap, solverData->Ai, solverData->Ax, &(solverData->symbolic), solverData->control, solverData->info); } @@ -281,9 +282,9 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) if (0 == status){ if (1 == systemData->method){ - status = umfpack_di_wsolve(UMFPACK_A, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); + status = umfpack_di_wsolve(UMFPACK_A, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->parDynamicData[omc_get_thread_num()].b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); } else { - status = umfpack_di_wsolve(UMFPACK_Aat, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); + status = umfpack_di_wsolve(UMFPACK_Aat, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->parDynamicData[omc_get_thread_num()].b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); } } @@ -371,8 +372,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) */ int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x) { - - DATA_UMFPACK* solverData = (DATA_UMFPACK*) systemData->solverData[0]; + DATA_UMFPACK* solverData = (DATA_UMFPACK*) systemData->parDynamicData[omc_get_thread_num()].solverData[0]; double *Ux, *Rs, r_ii, *b, sum, *y, *z; int *Up, *Ui, *Q, do_recip, rank = 0, current_rank, current_unz, i, j, k, l, success = 0, status, stop = 0; @@ -411,14 +411,13 @@ int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x) { for (i = 0; i < solverData->n_row; i++) { - b[i] = systemData->b[i] / Rs[i]; + b[i] = systemData->parDynamicData[omc_get_thread_num()].b[i] / Rs[i]; } } else { - for (i = 0; i < solverData->n_row; i++) - { - b[i] = systemData->b[i] * Rs[i]; + for (i = 0; i < solverData->n_row; i++) { + b[i] = systemData->parDynamicData[omc_get_thread_num()].b[i] * Rs[i]; } } diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.c b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.c index cc0e1a9787d..e2d97095fd2 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.c @@ -36,6 +36,8 @@ #include "model_help.h" #include "../../util/omc_error.h" +#include "../../util/parallel_helper.h" +#include "../../util/jacobian_util.h" #include "../../util/rtclock.h" #include "nonlinearSystem.h" #include "linearSystem.h" @@ -68,8 +70,12 @@ int initializeLinearSystems(DATA *data, threadData_t *threadData) TRACE_PUSH int i, nnz; int size; + int res; + unsigned int j, maxNumberThreads; LINEAR_SYSTEM_DATA *linsys = data->simulationInfo->linearSystemData; + maxNumberThreads = omc_get_max_threads(); + infoStreamPrint(LOG_LS, 1, "initialize linear system solvers"); infoStreamPrint(LOG_LS, 0, "%ld linear systems", data->modelData->nLinearSystems); @@ -83,23 +89,28 @@ int initializeLinearSystems(DATA *data, threadData_t *threadData) for(i=0; imodelData->nLinearSystems; ++i) { + res = allocLinSystThreadData(&(linsys[i])); + assertStreamPrint(threadData, !res ,"Out of memory"); + size = linsys[i].size; nnz = linsys[i].nnz; - linsys[i].totalTime = 0; linsys[i].failed = 0; /* allocate system data */ - linsys[i].b = (double*) malloc(size*sizeof(double)); + for (j=0; jsimulationInfo->analyticJacobians[linsys[i].jacobianIndex]); if(linsys[i].jacobianIndex != -1) { assertStreamPrint(threadData, 0 != linsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" ); } + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[linsys[i].jacobianIndex]); if(linsys[i].initialAnalyticalJacobian(data, threadData, jacobian)) { linsys[i].jacobianIndex = -1; @@ -107,6 +118,25 @@ int initializeLinearSystems(DATA *data, threadData_t *threadData) } nnz = jacobian->sparsePattern->numberOfNoneZeros; linsys[i].nnz = nnz; + +#ifdef USE_PARJAC + /* Allocate jacobian for parDynamicData */ + for (j=0; jsizeCols = jacobian->sizeCols; + linsys[i].parDynamicData[j].jacobian->sizeRows = jacobian->sizeRows; + linsys[i].parDynamicData[j].jacobian->sizeTmpVars = jacobian->sizeTmpVars; + linsys[i].parDynamicData[j].jacobian->seedVars = (modelica_real*) calloc(jacobian->sizeCols, sizeof(modelica_real)); + linsys[i].parDynamicData[j].jacobian->resultVars = (modelica_real*) calloc(jacobian->sizeRows, sizeof(modelica_real)); + linsys[i].parDynamicData[j].jacobian->tmpVars = (modelica_real*) calloc(jacobian->sizeTmpVars, sizeof(modelica_real)); + linsys[i].parDynamicData[j].jacobian->sparsePattern = jacobian->sparsePattern; + } +#else + linsys[i].parDynamicData[0].jacobian = jacobian; +#endif } if(nnz/(double)(size*size)<=linearSparseSolverMaxDensity && size>=linearSparseSolverMinSize) @@ -132,12 +162,18 @@ int initializeLinearSystems(DATA *data, threadData_t *threadData) case LSS_UMFPACK: linsys[i].setAElement = setAElementUmfpack; linsys[i].setBElement = setBElement; - allocateUmfPackData(size, size, nnz, linsys[i].solverData); + for (j=0; jsimulationInfo->lsMethod) { case LS_LAPACK: - linsys[i].A = (double*) malloc(size*size*sizeof(double)); linsys[i].setAElement = setAElement; linsys[i].setBElement = setBElement; - allocateLapackData(size, linsys[i].solverData); + for (j=0; jparDynamicData = (LINEAR_SYSTEM_THREAD_DATA*) malloc(omc_get_max_threads()*sizeof(LINEAR_SYSTEM_THREAD_DATA)); + if (!linsys->parDynamicData) + return -1; + return 0; +} + +/** + * \brief Frees memory allocated with `allocLinSystThreadData` + * + * Frees complete array of numberOfThreads many LINEAR_SYSTEM_THREAD_DATA structs. + * + * \param [in/out] linsys + */ +void freeLinSystThreadData(LINEAR_SYSTEM_DATA *linsys) +{ + free(linsys->parDynamicData); +} + /*! \fn int updateStaticDataOfLinearSystems(DATA *data) * * This function allocates memory for all linear systems. @@ -287,34 +372,67 @@ void printLinearSystemSolvingStatistics(DATA *data, int sysNumber, int logLevel) int freeLinearSystems(DATA *data, threadData_t *threadData) { TRACE_PUSH - int i; + int i,j; LINEAR_SYSTEM_DATA* linsys = data->simulationInfo->linearSystemData; infoStreamPrint(LOG_LS_V, 1, "free linear system solvers"); for(i=0; imodelData->nLinearSystems; ++i) { /* free system and solver data */ - free(linsys[i].b); + for (j=0; jsimulationInfo->analyticJacobians[linsys[i].jacobianIndex]); + freeAnalyticJacobian(jacobian); + /* Note: The Jacobian of data->simulationInfo itself will be free later. */ + +#ifdef USE_PARJAC + for (j=0; jsimulationInfo->analyticJacobians[linsys[i].jacobianIndex] + // which is free some lines above (and are invalid pointers at this point). Thus, free + // what is left. + free(linsys[i].parDynamicData[j].jacobian->seedVars); linsys[i].parDynamicData[j].jacobian->seedVars = NULL; + free(linsys[i].parDynamicData[j].jacobian->resultVars); linsys[i].parDynamicData[j].jacobian->resultVars = NULL; + free(linsys[i].parDynamicData[j].jacobian->tmpVars); linsys[i].parDynamicData[j].jacobian->tmpVars = NULL; + linsys[i].parDynamicData[j].jacobian->sparsePattern = NULL; + free(linsys[i].parDynamicData[j].jacobian); linsys[i].parDynamicData[j].jacobian = NULL; + } +#endif + } + if(linsys[i].useSparseSolver == 1) { switch(data->simulationInfo->lssMethod) { #if !defined(OMC_MINIMAL_RUNTIME) case LSS_LIS: - freeLisData(linsys[i].solverData); + for (j=0; jsimulationInfo->lssMethod); } } - - else{ - switch(data->simulationInfo->lsMethod) - { + else { // useSparseSolver == 0 + switch(data->simulationInfo->lsMethod) { case LS_LAPACK: - freeLapackData(linsys[i].solverData); - free(linsys[i].A); + for (j=0; jsimulationInfo->linearSystemData[sysNumber]); @@ -577,7 +702,11 @@ int check_linear_solution(DATA *data, int printFailingSystems, int sysNumber) TRACE_POP return 1; } +#ifdef USE_PARJAC + warningStreamPrintWithEquationIndexes(LOG_STDOUT, 1, indexes, "Thread %u: Solving linear system %d fails at time %g. For more information use -lv LOG_LS.", omc_get_thread_num(), index, data->localData[0]->timeValue); +#else warningStreamPrintWithEquationIndexes(LOG_STDOUT, 1, indexes, "Solving linear system %d fails at time %g. For more information use -lv LOG_LS.", index, data->localData[0]->timeValue); +#endif for(j=0; jmodelData->modelDataXml, (linsys[i]).equationIndex).numVar; ++j) { int done=0; @@ -629,7 +758,7 @@ int check_linear_solution(DATA *data, int printFailingSystems, int sysNumber) static void setAElement(int row, int col, double value, int nth, void *data, threadData_t *threadData) { LINEAR_SYSTEM_DATA* linsys = (LINEAR_SYSTEM_DATA*) data; - linsys->A[row + col * linsys->size] = value; + linsys->parDynamicData[omc_get_thread_num()].A[row + col * linsys->size] = value; } /*! \fn setBElement @@ -642,21 +771,21 @@ static void setAElement(int row, int col, double value, int nth, void *data, thr static void setBElement(int row, double value, void *data, threadData_t *threadData) { LINEAR_SYSTEM_DATA* linsys = (LINEAR_SYSTEM_DATA*) data; - linsys->b[row] = value; + linsys->parDynamicData[omc_get_thread_num()].b[row] = value; } #if !defined(OMC_MINIMAL_RUNTIME) static void setAElementLis(int row, int col, double value, int nth, void *data, threadData_t *threadData) { LINEAR_SYSTEM_DATA* linsys = (LINEAR_SYSTEM_DATA*) data; - DATA_LIS* sData = (DATA_LIS*) linsys->solverData[0]; + DATA_LIS* sData = (DATA_LIS*) linsys->parDynamicData[omc_get_thread_num()].solverData[0]; lis_matrix_set_value(LIS_INS_VALUE, row, col, value, sData->A); } static void setBElementLis(int row, double value, void *data, threadData_t *threadData) { LINEAR_SYSTEM_DATA* linsys = (LINEAR_SYSTEM_DATA*) data; - DATA_LIS* sData = (DATA_LIS*) linsys->solverData[0]; + DATA_LIS* sData = (DATA_LIS*) linsys->parDynamicData[omc_get_thread_num()].solverData[0]; lis_vector_set_value(LIS_INS_VALUE, row, value, sData->b); } #endif @@ -665,11 +794,11 @@ static void setBElementLis(int row, double value, void *data, threadData_t *thre static void setAElementUmfpack(int row, int col, double value, int nth, void *data, threadData_t *threadData) { LINEAR_SYSTEM_DATA* linSys = (LINEAR_SYSTEM_DATA*) data; - DATA_UMFPACK* sData = (DATA_UMFPACK*) linSys->solverData[0]; + DATA_UMFPACK* sData = (DATA_UMFPACK*) linSys->parDynamicData[omc_get_thread_num()].solverData[0]; infoStreamPrint(LOG_LS_V, 0, " set %d. -> (%d,%d) = %f", nth, row, col, value); - if (row > 0){ - if (sData->Ap[row] == 0){ + if (row > 0) { + if (sData->Ap[row] == 0) { sData->Ap[row] = nth; } } @@ -677,13 +806,14 @@ static void setAElementUmfpack(int row, int col, double value, int nth, void *da sData->Ai[nth] = col; sData->Ax[nth] = value; } + static void setAElementKlu(int row, int col, double value, int nth, void *data, threadData_t *threadData) { LINEAR_SYSTEM_DATA* linSys = (LINEAR_SYSTEM_DATA*) data; - DATA_KLU* sData = (DATA_KLU*) linSys->solverData[0]; + DATA_KLU* sData = (DATA_KLU*) linSys->parDynamicData[omc_get_thread_num()].solverData[0]; - if (row > 0){ - if (sData->Ap[row] == 0){ + if (row > 0) { + if (sData->Ap[row] == 0) { sData->Ap[row] = nth; } } @@ -692,5 +822,4 @@ static void setAElementKlu(int row, int col, double value, int nth, void *data, sData->Ax[nth] = value; } - #endif diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.h b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.h index d7c91b87f1d..36ba937fb99 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.h +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.h @@ -49,6 +49,7 @@ extern "C" { typedef void* LS_SOLVER_DATA; int initializeLinearSystems(DATA *data, threadData_t *threadData); +int allocLinSystThreadData(LINEAR_SYSTEM_DATA *linsys); int updateStaticDataOfLinearSystems(DATA *data, threadData_t *threadData); int freeLinearSystems(DATA *data, threadData_t *threadData); int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/model_help.c b/OMCompiler/SimulationRuntime/c/simulation/solver/model_help.c index f19f2f8f667..772651b810d 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/model_help.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/model_help.c @@ -52,6 +52,10 @@ #include "stateset.h" #include "../../meta/meta_modelica.h" +#ifdef USE_PARJAC + #include +#endif + int maxEventIterations = 20; double linearSparseSolverMaxDensity = 0.2; int linearSparseSolverMinSize = 201; diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c b/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c index 97a18fea9d5..1b89234a86c 100644 --- a/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c +++ b/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c @@ -43,6 +43,7 @@ #include "ida_solver.h" #include "delay.h" #include "events.h" +#include "util/parallel_helper.h" #include "util/varinfo.h" #include "radau.h" #include "model_help.h" @@ -360,6 +361,7 @@ int freeSolverData(DATA* data, SOLVER_INFO* solverInfo) int retValue = 0; int i; + freeList(solverInfo->eventLst); /* free solver statistics */ free(solverInfo->solverStats); free(solverInfo->solverStatsTmp); @@ -589,6 +591,15 @@ int finishSimulation(DATA* data, threadData_t *threadData, SOLVER_INFO* solverIn infoStreamPrint(LOG_STATS, 0, "%5d error test failures", solverInfo->solverStats[3]); infoStreamPrint(LOG_STATS, 0, "%5d convergence test failures", solverInfo->solverStats[4]); + infoStreamPrint(LOG_STATS, 0, "%gs time of jacobian evaluation", rt_accumulated(SIM_TIMER_JACOBIAN)); +#ifdef USE_PARJAC + infoStreamPrint(LOG_STATS, 0, "%i OpenMP-threads used for jacobian evaluation", omc_get_max_threads()); + int chunk_size; + omp_sched_t kind; + omp_get_schedule(&kind, &chunk_size); + infoStreamPrint(LOG_STATS, 0, "Schedule: %i Chunk Size: %i", kind, chunk_size); +#endif + messageClose(LOG_STATS); } diff --git a/OMCompiler/SimulationRuntime/c/simulation_data.h b/OMCompiler/SimulationRuntime/c/simulation_data.h index 438f04d666c..c67f8b30b06 100644 --- a/OMCompiler/SimulationRuntime/c/simulation_data.h +++ b/OMCompiler/SimulationRuntime/c/simulation_data.h @@ -120,6 +120,8 @@ typedef enum {ERROR_AT_TIME,NO_PROGRESS_START_POINT,NO_PROGRESS_FACTOR,IMPROPER_ * sizeofIndex contain number of elements in index * colorsCols contain color of colored columns * + * Use freeSparsePattern(SPARSE_PATTERM *spp) for "destruction" (see util/jacobian_util.c/h). + * */ typedef struct SPARSE_PATTERN { @@ -142,6 +144,7 @@ typedef struct SPARSE_PATTERN * resultVars contain result of one column to the corresponding jacobian * jacobian contains dense jacobian elements * + * Use freeAnalyticJacobian(ANALYTIC_JACOBIAN *jac) for "destruction" (see util/jacobian_util.c/h). */ typedef struct ANALYTIC_JACOBIAN { @@ -309,6 +312,28 @@ typedef struct NONLINEAR_SYSTEM_DATA typedef void* NONLINEAR_SYSTEM_DATA; #endif +typedef struct LINEAR_SYSTEM_THREAD_DATA +{ + void *solverData[2]; /* [1] is the totalPivot solver + [0] holds other solvers + both are used for the default solver */ + modelica_real *x; /* solution vector x */ + modelica_real *A; /* matrix A */ + modelica_real *b; /* vector b */ + + modelica_real residualError; /* not used yet*/ + + ANALYTIC_JACOBIAN* parentJacobian; /* if != NULL then it's the parent jacobian matrix */ + ANALYTIC_JACOBIAN* jacobian; /* jacobian */ + + /* Statistics for each thread */ + unsigned long numberOfCall; /* number of solving calls of this system */ + unsigned long numberOfJEval; /* number of jacobian evaluations of this system */ + double totalTime; /* save the totalTime */ + rtclock_t totalTimeClock; /* time clock for the totalTime */ + double jacobianTime; /* save the time to calculate jacobians */ +}LINEAR_SYSTEM_THREAD_DATA; + #if !defined(OMC_NUM_LINEAR_SYSTEMS) || OMC_NUM_LINEAR_SYSTEMS>0 typedef struct LINEAR_SYSTEM_DATA { @@ -322,7 +347,6 @@ typedef struct LINEAR_SYSTEM_DATA int (*analyticalJacobianColumn)(void*, threadData_t*, ANALYTIC_JACOBIAN*, ANALYTIC_JACOBIAN* parentJacobian); int (*initialAnalyticalJacobian)(void*, threadData_t*, ANALYTIC_JACOBIAN*); - modelica_integer jacobianIndex; void (*residualFunc)(void**, const double*, double*, const int*); void (*initializeStaticLSData)(void*, threadData_t *threadData, void*); @@ -337,19 +361,18 @@ typedef struct LINEAR_SYSTEM_DATA modelica_integer nnz; /* number of nonzero entries */ modelica_integer size; modelica_integer equationIndex; /* index for EQUATION_INFO */ - - void *solverData[2]; /* [1] is the totalPivot solver; [0] holds other solvers ; both are used for the default solver */ - modelica_real *x; /* solution vector x */ - modelica_real *A; /* matrix A */ - modelica_real *b; /* vector b */ + modelica_integer jacobianIndex; modelica_integer method; /* not used yet*/ - modelica_real residualError; /* not used yet*/ + modelica_boolean useSparseSolver; /* 1: use sparse solver, - else any solver */ + + LINEAR_SYSTEM_THREAD_DATA* parDynamicData; /* Array of length numMaxThreads for internal write data */ + + // ToDo: Gather information from all threads if in parallel region modelica_boolean solved; /* 1: solved in current step - else not */ modelica_boolean failed; /* 1: failed while last try with lapack - else not */ - modelica_boolean useSparseSolver; /* 1: use sparse solver, - else any solver */ - ANALYTIC_JACOBIAN* parentJacobian; /* if != NULL then it's the parent jacobian matrix */ + // ToDo: Gather information from all threads if in parallel region /* statistics */ unsigned long numberOfCall; /* number of solving calls of this system */ unsigned long numberOfJEval; /* number of jacobian evaluations of this system */ @@ -377,7 +400,7 @@ typedef struct MIXED_SYSTEM_DATA modelica_boolean** iterationPreVarsPtr; void *solverData; - modelica_integer method; /* not used yet*/ + modelica_integer method; /* not used yet */ modelica_boolean solved; /* 1: solved in current step - else not */ }MIXED_SYSTEM_DATA; #else @@ -659,7 +682,7 @@ typedef struct SIMULATION_INFO modelica_real* sensitivityMatrix; /* used by integrator for sensitivity mode */ int* sensitivityParList; /* used by integrator for sensitivity mode */ - ANALYTIC_JACOBIAN* analyticJacobians; + ANALYTIC_JACOBIAN* analyticJacobians; // ToDo Only store informations for Jacobian used by integrator here NONLINEAR_SYSTEM_DATA* nonlinearSystemData; int currentNonlinearSystemIndex; diff --git a/OMCompiler/SimulationRuntime/c/util/CMakeLists.txt b/OMCompiler/SimulationRuntime/c/util/CMakeLists.txt index db80f36ab10..f1cf3c0efa5 100644 --- a/OMCompiler/SimulationRuntime/c/util/CMakeLists.txt +++ b/OMCompiler/SimulationRuntime/c/util/CMakeLists.txt @@ -1,16 +1,70 @@ -# Quellen und Header -SET(util_sources base_array.c boolean_array.c omc_file.c omc_error.c division.c index_spec.c - integer_array.c java_interface.c libcsv.c list.c modelica_string.c - read_write.c read_matlab4.c read_csv.c real_array.c ringbuffer.c rational.c - rtclock.c simulation_options.c string_array.c utility.c varinfo.c omc_msvc.c OldModelicaTables.c omc_mmap.c - ModelicaUtilities.c modelica_string_lit.c omc_init.c write_csv.c ../gc/memory_pool.c) +# Sources and headers +SET(util_sources ../gc/memory_pool.c + base_array.c + boolean_array.c + division.c + index_spec.c + integer_array.c + jacobian_util.c + java_interface.c + libcsv.c + list.c + modelica_string_lit.c + modelica_string.c + ModelicaUtilities.c + OldModelicaTables.c + omc_error.c + omc_file.c + omc_init.c + omc_mmap.c + omc_msvc.c + parallel_helper.c + rational.c + read_csv.c + read_matlab4.c + read_write.c + real_array.c + ringbuffer.c + rtclock.c + simulation_options.c + string_array.c + utility.c + varinfo.c + write_csv.c) -SET(util_headers base_array.h boolean_array.h division.h omc_file.h omc_error.h index_spec.h integer_array.h - java_interface.h jni.h jni_md.h jni_md_solaris.h jni_md_windows.h list.h - modelica.h modelica_string.h read_write.h read_matlab4.h real_array.h rational.h - ringbuffer.h rtclock.h simulation_options.h string_array.h utility.h varinfo.h omc_mmap.h - ../ModelicaUtilities.h modelica_string_lit.h omc_init.h write_csv.h ../gc/memory_pool.h) +SET(util_headers ../gc/memory_pool.h + ../ModelicaUtilities.h + base_array.h + boolean_array.h + division.h + index_spec.h + integer_array.h + jacobian_util.h + java_interface.h + jni_md_solaris.h + jni_md_windows.h + jni_md.h + jni.h + list.h + modelica_string_lit.h + modelica_string.h + modelica.h + omc_error.h + omc_file.h + omc_init.h write_csv.h + omc_mmap.h + parallel_helper.h + rational.h + read_matlab4.h + read_write.h + real_array.h + ringbuffer.h + rtclock.h + simulation_options.h + string_array.h + utility.h + varinfo.h) if(MSVC) INCLUDE_DIRECTORIES(${OMCTRUNCHOME}/OMCompiler/3rdParty/regex-0.12) diff --git a/OMCompiler/SimulationRuntime/c/util/jacobian_util.c b/OMCompiler/SimulationRuntime/c/util/jacobian_util.c new file mode 100644 index 00000000000..4d4f1d67c90 --- /dev/null +++ b/OMCompiler/SimulationRuntime/c/util/jacobian_util.c @@ -0,0 +1,48 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2019, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE + * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: http://www.openmodelica.org or + * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica + * distribution. GNU version 3 is obtained from: + * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from: + * http://www.opensource.org/licenses/BSD-3-Clause. + * + * This program is distributed WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS + * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE + * CONDITIONS OF OSMC-PL. + * + */ + +/*! File jac_util.c + */ + +#include "jacobian_util.h" + +void freeAnalyticJacobian(ANALYTIC_JACOBIAN *jac) { + free(jac->seedVars); jac->seedVars = NULL; + free(jac->tmpVars); jac->tmpVars = NULL; + free(jac->resultVars); jac->resultVars = NULL; + freeSparsePattern(jac->sparsePattern); + free(jac->sparsePattern); jac->sparsePattern = NULL; +} + +void freeSparsePattern(SPARSE_PATTERN *spp) { + free(spp->index); spp->index = NULL; + free(spp->colorCols); spp->colorCols = NULL; + free(spp->leadindex); spp->leadindex = NULL; +} diff --git a/OMCompiler/SimulationRuntime/c/util/jacobian_util.h b/OMCompiler/SimulationRuntime/c/util/jacobian_util.h new file mode 100644 index 00000000000..f8c47434fcc --- /dev/null +++ b/OMCompiler/SimulationRuntime/c/util/jacobian_util.h @@ -0,0 +1,61 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2019, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE + * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: http://www.openmodelica.org or + * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica + * distribution. GNU version 3 is obtained from: + * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from: + * http://www.opensource.org/licenses/BSD-3-Clause. + * + * This program is distributed WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS + * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE + * CONDITIONS OF OSMC-PL. + * + */ + +/*! File jacobian_util.h + */ + +#ifndef OMC_JACOBIAN_UTIL_H +#define OMC_JACOBIAN_UTIL_H + +#include "../simulation_data.h" + +#ifdef __cplusplus +extern "C" { +#endif + +/** + * \brief Free struct ANALYTIC_JACOBIAN + * + * Frees dynamically allocated memory and sets pointers to NULL. + */ +void freeAnalyticJacobian(ANALYTIC_JACOBIAN* jac); + +/** + * \brief Free struct SPARSE_PATTERN + * + * Frees dynamically allocated memory and sets pointers to NULL. + */ +void freeSparsePattern(SPARSE_PATTERN *spp); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/OMCompiler/SimulationRuntime/c/util/parallel_helper.c b/OMCompiler/SimulationRuntime/c/util/parallel_helper.c new file mode 100644 index 00000000000..8049d66518d --- /dev/null +++ b/OMCompiler/SimulationRuntime/c/util/parallel_helper.c @@ -0,0 +1,64 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2019, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE + * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: http://www.openmodelica.org or + * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica + * distribution. GNU version 3 is obtained from: + * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from: + * http://www.opensource.org/licenses/BSD-3-Clause. + * + * This program is distributed WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS + * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE + * CONDITIONS OF OSMC-PL. + * + */ + +/*! File parallel_helper.c + */ + +#include "parallel_helper.h" + +/** + * \brief Wrapper for OpenMP function omp_get_thread_num + * + * If OpenMP is available return thread number, otherwise 0. + */ +int omc_get_thread_num(void) +{ +#ifdef USE_PARJAC + return omp_get_thread_num(); +#else + return 0; +#endif +} + +/** + * \brief Wrapper for OpenMP function omc_get_max_threads + * + * If OpenMP is available return maximum number of threads, + * otherwise 1. + */ +int omc_get_max_threads(void) +{ +#ifdef USE_PARJAC + return omp_get_max_threads(); +#else + return 1; +#endif +} + diff --git a/OMCompiler/SimulationRuntime/c/util/parallel_helper.h b/OMCompiler/SimulationRuntime/c/util/parallel_helper.h new file mode 100644 index 00000000000..86fe674c3e9 --- /dev/null +++ b/OMCompiler/SimulationRuntime/c/util/parallel_helper.h @@ -0,0 +1,53 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2014, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE + * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: http://www.openmodelica.org or + * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica + * distribution. GNU version 3 is obtained from: + * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from: + * http://www.opensource.org/licenses/BSD-3-Clause. + * + * This program is distributed WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS + * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE + * CONDITIONS OF OSMC-PL. + * + */ + +/*! File parallel_helper.h + */ + +#ifndef OMC_PARALLEL_HELPER_H +#define OMC_PARALLEL_HELPER_H + +#ifdef USE_PARJAC + #include +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +// TODO Add DLL_export to work on windows? +extern int omc_get_thread_num(void); +extern int omc_get_max_threads(void); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/OMCompiler/SimulationRuntime/c/util/simulation_options.c b/OMCompiler/SimulationRuntime/c/util/simulation_options.c index cab8f33239b..beccf7768a6 100644 --- a/OMCompiler/SimulationRuntime/c/util/simulation_options.c +++ b/OMCompiler/SimulationRuntime/c/util/simulation_options.c @@ -88,6 +88,7 @@ const char *FLAG_NAME[FLAG_MAX+1] = { /* FLAG_IPOPT_MAX_ITER */ "ipopt_max_iter", /* FLAG_IPOPT_WARM_START */ "ipopt_warm_start", /* FLAG_JACOBIAN */ "jacobian", + /* FLAG_JACOBIAN_THREADS */ "jacobianThreads", /* FLAG_L */ "l", /* FLAG_L_DATA_RECOVERY */ "l_datarec", /* FLAG_LOG_FORMAT */ "logFormat", @@ -201,6 +202,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = { /* FLAG_IPOPT_MAX_ITER */ "value specifies the max number of iteration for ipopt", /* FLAG_IPOPT_WARM_START */ "value specifies lvl for a warm start in ipopt: 1,2,3,...", /* FLAG_JACOBIAN */ "select the calculation method of the Jacobian used only by ida and dassl solver.", + /* FLAG_JACOBIAN_THREADS */ "[int default: 1] value specifies the number of threads for jacobian evaluation in dassl or ida.", /* FLAG_L */ "value specifies a time where the linearization of the model should be performed", /* FLAG_L_DATA_RECOVERY */ "emit data recovery matrices with model linearization", /* FLAG_LOG_FORMAT */ "value specifies the log format of the executable. -logFormat=text (default), -logFormat=xml or -logFormat=xmltcp", @@ -385,6 +387,9 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = { " Value specifies lvl for a warm start in ipopt: 1,2,3,...", /* FLAG_JACOBIAN */ " Select the calculation method for Jacobian used by the integration method:\n", + /* FLAG_JACOBIAN_THREADS */ + " Value specifies the number of threads for jacobian evaluation in dassl or ida." + " The value is an Integer with default value 1.", /* FLAG_L */ " Value specifies a time where the linearization of the model should be performed.", /* FLAG_L_DATA_RECOVERY */ @@ -586,6 +591,7 @@ const int FLAG_TYPE[FLAG_MAX] = { /* FLAG_IPOPT_MAX_ITER */ FLAG_TYPE_OPTION, /* FLAG_IPOPT_WARM_START */ FLAG_TYPE_OPTION, /* FLAG_JACOBIAN */ FLAG_TYPE_OPTION, + /* FLAG_JACOBIAN_THREADS */ FLAG_TYPE_OPTION, /* FLAG_L */ FLAG_TYPE_OPTION, /* FLAG_L_DATA_RECOVERY */ FLAG_TYPE_FLAG, /* FLAG_LOG_FORMAT */ FLAG_TYPE_OPTION, diff --git a/OMCompiler/SimulationRuntime/c/util/simulation_options.h b/OMCompiler/SimulationRuntime/c/util/simulation_options.h index 99852526bc6..93928551f35 100644 --- a/OMCompiler/SimulationRuntime/c/util/simulation_options.h +++ b/OMCompiler/SimulationRuntime/c/util/simulation_options.h @@ -96,6 +96,7 @@ enum _FLAG FLAG_IPOPT_MAX_ITER, FLAG_IPOPT_WARM_START, FLAG_JACOBIAN, + FLAG_JACOBIAN_THREADS, FLAG_L, FLAG_L_DATA_RECOVERY, FLAG_LOG_FORMAT, diff --git a/OMCompiler/configure.ac b/OMCompiler/configure.ac index d3dcdd77da2..1f3a90e3776 100644 --- a/OMCompiler/configure.ac +++ b/OMCompiler/configure.ac @@ -485,10 +485,10 @@ if test x$want_parjac = xyes; then USE_PARJAC="yes" RUNTIMECFLAGS="$RUNTIMECFLAGS $OMPCFLAGS -DUSE_PARJAC" else - AC_MSG_WARN([========= OpenMP is not available, not define par. jac]) + AC_MSG_WARN([========= OpenMP is not available, will evaluate jacobians sequentially]) fi + AC_MSG_RESULT([$USE_PARJAC]) fi -AC_MSG_RESULT([$USE_PARJAC]) AC_SUBST(USE_PARJAC) diff --git a/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug2764.mos b/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug2764.mos index b78c91c8d2a..3fa171b740a 100644 --- a/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug2764.mos +++ b/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug2764.mos @@ -91,11 +91,13 @@ readFile("modelDescription.tmp.xml"); // // // +// // // // // // +// // // // diff --git a/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug3049.mos b/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug3049.mos index d411c1972e8..b3e70fc1e65 100644 --- a/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug3049.mos +++ b/testsuite/openmodelica/fmi/ModelExchange/2.0/testBug3049.mos @@ -82,11 +82,13 @@ readFile("modelDescription.tmp.xml"); // // // +// // // // // // +// // // // diff --git a/testsuite/openmodelica/fmi/ModelExchange/2.0/testDisableDep.mos b/testsuite/openmodelica/fmi/ModelExchange/2.0/testDisableDep.mos index 29cf467b0d2..67408266740 100644 --- a/testsuite/openmodelica/fmi/ModelExchange/2.0/testDisableDep.mos +++ b/testsuite/openmodelica/fmi/ModelExchange/2.0/testDisableDep.mos @@ -96,11 +96,13 @@ readFile("modelDescription.tmp.xml"); // // // +// // // // // // +// // // // diff --git a/testsuite/openmodelica/fmi/ModelExchange/2.0/testDiscreteStructe.mos b/testsuite/openmodelica/fmi/ModelExchange/2.0/testDiscreteStructe.mos index 38986a3e3f7..ae39e02b3f2 100644 --- a/testsuite/openmodelica/fmi/ModelExchange/2.0/testDiscreteStructe.mos +++ b/testsuite/openmodelica/fmi/ModelExchange/2.0/testDiscreteStructe.mos @@ -93,11 +93,13 @@ readFile("modelDescription.tmp.xml"); // // // +// // // // // // +// // // //