From ca68297be7d409fa47aca740b55b33f9cac53965 Mon Sep 17 00:00:00 2001 From: Michele Scuttari Date: Thu, 31 Oct 2024 16:19:20 +0100 Subject: [PATCH] Remove matching information from IDA and KINSOL --- include/marco/Runtime/Solvers/IDA/Instance.h | 34 +-- .../marco/Runtime/Solvers/KINSOL/Instance.h | 34 +-- lib/Solvers/IDA/Instance.cpp | 222 ++++-------------- lib/Solvers/KINSOL/Instance.cpp | 221 ++++------------- 4 files changed, 92 insertions(+), 419 deletions(-) diff --git a/include/marco/Runtime/Solvers/IDA/Instance.h b/include/marco/Runtime/Solvers/IDA/Instance.h index e756726c..641268a9 100644 --- a/include/marco/Runtime/Solvers/IDA/Instance.h +++ b/include/marco/Runtime/Solvers/IDA/Instance.h @@ -79,8 +79,6 @@ class IDAInstance { /// Add the information about an equation that is handled by IDA. Equation addEquation(const int64_t *ranges, uint64_t rank, - Variable writtenVariable, - AccessFunction writeAccessFunction, const char *stringRepresentation); void addVariableAccess(Equation equation, Variable variableIndex, @@ -157,10 +155,6 @@ class IDAInstance { [[nodiscard]] uint64_t getEquationFlatSize(Equation equation) const; - [[nodiscard]] Variable getWrittenVariable(Equation equation) const; - - [[nodiscard]] AccessFunction getWriteAccessFunction(Equation equation) const; - [[nodiscard]] uint64_t getVariableRank(Variable variable) const; void @@ -224,17 +218,6 @@ class IDAInstance { bool idaSetMaxNumItersIC(); bool idaSetLineSearchOffIC(); - /// } - /// @name Utility functions - /// { - - /// Get the scalar equation writing to a certain scalar variable. - /// Warning: extremely slow, to be used only for debug purposes. - void getWritingEquation(Variable variable, - const std::vector &variableIndices, - Equation &equation, - std::vector &equationIndices) const; - /// } /// @name Debug functions /// { @@ -265,17 +248,6 @@ class IDAInstance { // The iteration ranges of the vectorized equations. std::vector equationRanges; - // The array variables written by the equations. - // The i-th position contains the information about the variable written - // by the i-th equation: the first element is the index of the IDA - // variable, while the second represents the ranges of the scalar - // variable. - std::vector> writeAccesses; - - // The order in which the equations must be processed when computing - // residuals and partial derivatives. - std::vector equationsProcessingOrder; - // The residual functions associated with the equations. // The i-th position contains the pointer to the residual function of the // i-th equation. @@ -300,6 +272,10 @@ class IDAInstance { // The dimensions list of each array variable. std::vector variablesDimensions; + // The offset of each array equation inside the flattened equations + // vector. + std::vector equationOffsets; + // Simulation times. realtype startTime; realtype endTime; @@ -389,7 +365,7 @@ RUNTIME_FUNC_DECL(idaAddVariableAccess, void, PTR(void), uint64_t, uint64_t, PTR(void)) RUNTIME_FUNC_DECL(idaAddEquation, uint64_t, PTR(void), PTR(int64_t), uint64_t, - uint64_t, PTR(void), PTR(void)) + PTR(void)) RUNTIME_FUNC_DECL(idaSetResidual, void, PTR(void), uint64_t, PTR(void)) diff --git a/include/marco/Runtime/Solvers/KINSOL/Instance.h b/include/marco/Runtime/Solvers/KINSOL/Instance.h index 313b3643..8559dc20 100644 --- a/include/marco/Runtime/Solvers/KINSOL/Instance.h +++ b/include/marco/Runtime/Solvers/KINSOL/Instance.h @@ -63,8 +63,6 @@ class KINSOLInstance { /// Add the information about an equation that is handled by KINSOL. Equation addEquation(const int64_t *ranges, uint64_t rank, - Variable writtenVariable, - AccessFunction writeAccessFunction, const char *stringRepresentation); void addVariableAccess(Equation equation, Variable variableIndex, @@ -117,10 +115,6 @@ class KINSOLInstance { [[nodiscard]] uint64_t getEquationFlatSize(Equation equation) const; - [[nodiscard]] Variable getWrittenVariable(Equation equation) const; - - [[nodiscard]] AccessFunction getWriteAccessFunction(Equation equation) const; - [[nodiscard]] uint64_t getVariableRank(Variable variable) const; void @@ -167,17 +161,6 @@ class KINSOLInstance { bool kinsolSetUserData(); bool kinsolSetJacobianFunction(); - /// } - /// @name Utility functions - /// { - - /// Get the scalar equation writing to a certain scalar variable. - /// Warning: extremely slow, to be used only for debug purposes. - void getWritingEquation(Variable variable, - const std::vector &variableIndices, - Equation &equation, - std::vector &equationIndices) const; - /// } /// @name Debug functions /// { @@ -206,17 +189,6 @@ class KINSOLInstance { // The iteration ranges of the vectorized equations. std::vector equationRanges; - // The array variables written by the equations. - // The i-th position contains the information about the variable written - // by the i-th equation: the first element is the index of the IDA - // variable, while the second represents the ranges of the scalar - // variable. - std::vector> writeAccesses; - - // The order in which the equations must be processed when computing - // residuals and partial derivatives. - std::vector equationsProcessingOrder; - // The residual functions associated with the equations. // The i-th position contains the pointer to the residual function of the // i-th equation. @@ -241,6 +213,10 @@ class KINSOLInstance { // The dimensions list of each array variable. std::vector variablesDimensions; + // The offset of each array equation inside the flattened equations + // vector. + std::vector equationOffsets; + // Variables vectors and values. N_Vector variablesVector; @@ -300,7 +276,7 @@ RUNTIME_FUNC_DECL(kinsolAddVariableAccess, void, PTR(void), uint64_t, uint64_t, PTR(void)) RUNTIME_FUNC_DECL(kinsolAddEquation, uint64_t, PTR(void), PTR(int64_t), - uint64_t, uint64_t, PTR(void), PTR(void)) + uint64_t, PTR(void)) RUNTIME_FUNC_DECL(kinsolSetResidual, void, PTR(void), uint64_t, PTR(void)) diff --git a/lib/Solvers/IDA/Instance.cpp b/lib/Solvers/IDA/Instance.cpp index f26c2681..fc1d4690 100644 --- a/lib/Solvers/IDA/Instance.cpp +++ b/lib/Solvers/IDA/Instance.cpp @@ -27,8 +27,9 @@ IDAInstance::IDAInstance() : startTime(simulation::getOptions().startTime), endTime(simulation::getOptions().endTime), timeStep(getOptions().timeStep) { - // Initially there is no variable in the instance. + // Initially there is are no variables or equations in the instance. variableOffsets.push_back(0); + equationOffsets.push_back(0); if (marco::runtime::simulation::getOptions().debug) { std::cerr << "[IDA] Instance created" << std::endl; @@ -220,8 +221,6 @@ Variable IDAInstance::addStateVariable(uint64_t rank, } Equation IDAInstance::addEquation(const int64_t *ranges, uint64_t equationRank, - Variable writtenVariable, - AccessFunction writeAccess, const char *stringRepresentation) { if (marco::runtime::simulation::getOptions().debug) { std::cerr << "[IDA] Adding equation"; @@ -235,17 +234,18 @@ Equation IDAInstance::addEquation(const int64_t *ranges, uint64_t equationRank, // Add the start and end dimensions of the current equation. MultidimensionalRange eqRanges = {}; + uint64_t flatSize = 1; for (size_t i = 0, e = equationRank * 2; i < e; i += 2) { int64_t begin = ranges[i]; int64_t end = ranges[i + 1]; eqRanges.push_back({begin, end}); + flatSize *= end - begin; } equationRanges.push_back(eqRanges); - - // Add the write access. - writeAccesses.emplace_back(writtenVariable, writeAccess); + size_t offset = equationOffsets.back(); + equationOffsets.push_back(offset + flatSize); // Return the index of the equation. Equation id = getNumOfVectorizedEquations() - 1; @@ -265,9 +265,6 @@ Equation IDAInstance::addEquation(const int64_t *ranges, uint64_t equationRank, } std::cerr << "]" << std::endl; - std::cerr << " - Written variable ID: " << writtenVariable << std::endl; - std::cerr << " - Write access function address: " - << reinterpret_cast(writeAccess) << std::endl; } return id; @@ -469,29 +466,6 @@ bool IDAInstance::initialize() { } } - // Determine the order in which the equations must be processed when - // computing residuals and jacobians. - assert(getNumOfVectorizedEquations() == writeAccesses.size()); - equationsProcessingOrder.resize(getNumOfVectorizedEquations()); - - for (size_t i = 0, e = getNumOfVectorizedEquations(); i < e; ++i) { - equationsProcessingOrder[i] = i; - } - - if (marco::runtime::simulation::getOptions().debug) { - std::cerr << "[IDA] Equations processing order: ["; - - for (size_t i = 0, e = equationsProcessingOrder.size(); i < e; ++i) { - if (i != 0) { - std::cerr << ", "; - } - - std::cerr << equationsProcessingOrder[i]; - } - - std::cerr << "]" << std::endl; - } - // Check that all the residual functions have been set. assert(residualFunctions.size() == getNumOfVectorizedEquations()); @@ -552,35 +526,19 @@ bool IDAInstance::initialize() { jacobianMatrixData.resize(scalarEquationsNumber); - for (Equation eq : equationsProcessingOrder) { + uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); + + for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { std::vector equationIndices; getEquationBeginIndices(eq, equationIndices); - Variable writtenVariable = getWrittenVariable(eq); - - uint64_t writtenVariableRank = getVariableRank(writtenVariable); - uint64_t writtenVariableArrayOffset = variableOffsets[writtenVariable]; - do { - std::vector writtenVariableIndices; - writtenVariableIndices.resize(writtenVariableRank, 0); - - AccessFunction writeAccessFunction = getWriteAccessFunction(eq); - - writeAccessFunction(equationIndices.data(), - writtenVariableIndices.data()); - - if (marco::runtime::simulation::getOptions().debug) { - std::cerr << " Variable indices: "; - printIndices(writtenVariableIndices); - std::cerr << std::endl; - } + uint64_t equationArrayOffset = equationOffsets[eq]; - uint64_t equationScalarVariableOffset = getVariableFlatIndex( - variablesDimensions[writtenVariable], writtenVariableIndices); + uint64_t equationScalarOffset = + getEquationFlatIndex(equationIndices, equationRanges[eq]); - uint64_t scalarEquationIndex = - writtenVariableArrayOffset + equationScalarVariableOffset; + uint64_t scalarEquationIndex = equationArrayOffset + equationScalarOffset; // Compute the column indexes that may be non-zeros. std::vector jacobianColumns = @@ -594,10 +552,6 @@ bool IDAInstance::initialize() { printIndices(equationIndices); std::cerr << std::endl; - std::cerr << " Variable indices: "; - printIndices(writtenVariableIndices); - std::cerr << std::endl; - std::cerr << " Scalar equation index: " << scalarEquationIndex << std::endl; @@ -911,28 +865,12 @@ int IDAInstance::residualFunction(realtype time, N_Vector variables, uint64_t equationRank = instance->getEquationRank(eq); assert(equationIndices.size() == equationRank); - Variable writtenVariable = instance->getWrittenVariable(eq); - - uint64_t writtenVariableArrayOffset = - instance->variableOffsets[writtenVariable]; + uint64_t equationArrayOffset = instance->equationOffsets[eq]; - uint64_t writtenVariableRank = - instance->getVariableRank(writtenVariable); + uint64_t equationScalarOffset = + getEquationFlatIndex(equationIndices, instance->equationRanges[eq]); - std::vector writtenVariableIndices(writtenVariableRank, 0); - - AccessFunction writeAccessFunction = - instance->getWriteAccessFunction(eq); - - writeAccessFunction(equationIndices.data(), - writtenVariableIndices.data()); - - uint64_t writtenVariableScalarOffset = - getVariableFlatIndex(instance->variablesDimensions[writtenVariable], - writtenVariableIndices); - - uint64_t offset = - writtenVariableArrayOffset + writtenVariableScalarOffset; + uint64_t offset = equationArrayOffset + equationScalarOffset; auto residualFn = instance->residualFunctions[eq]; auto *eqIndicesPtr = equationIndices.data(); @@ -978,29 +916,13 @@ int IDAInstance::jacobianMatrix(realtype time, realtype alpha, instance->equationsParallelIteration( [&](Equation eq, const std::vector &equationIndices, const JacobianSeedsMap &jacobianSeedsMap) { - Variable writtenVariable = instance->getWrittenVariable(eq); - - uint64_t writtenVariableArrayOffset = - instance->variableOffsets[writtenVariable]; - - uint64_t writtenVariableRank = - instance->getVariableRank(writtenVariable); - - std::vector writtenVariableIndices; - writtenVariableIndices.resize(writtenVariableRank, 0); - - AccessFunction writeAccessFunction = - instance->getWriteAccessFunction(eq); - - writeAccessFunction(equationIndices.data(), - writtenVariableIndices.data()); + uint64_t equationArrayOffset = instance->equationOffsets[eq]; - uint64_t writtenVariableScalarOffset = - getVariableFlatIndex(instance->variablesDimensions[writtenVariable], - writtenVariableIndices); + uint64_t equationScalarOffset = + getEquationFlatIndex(equationIndices, instance->equationRanges[eq]); uint64_t scalarEquationIndex = - writtenVariableArrayOffset + writtenVariableScalarOffset; + equationArrayOffset + equationScalarOffset; assert(scalarEquationIndex < instance->getNumOfScalarEquations()); @@ -1130,14 +1052,6 @@ uint64_t IDAInstance::getEquationFlatSize(Equation equation) const { return result; } -Variable IDAInstance::getWrittenVariable(Equation equation) const { - return writeAccesses[equation].first; -} - -AccessFunction IDAInstance::getWriteAccessFunction(Equation equation) const { - return writeAccesses[equation].second; -} - uint64_t IDAInstance::getVariableRank(Variable variable) const { return variablesDimensions[variable].rank(); } @@ -1279,7 +1193,7 @@ void IDAInstance::computeThreadChunks() { uint64_t processedEquations = 0; while (processedEquations < numOfVectorizedEquations) { - Equation equation = equationsProcessingOrder[processedEquations]; + Equation equation = processedEquations; uint64_t equationFlatSize = getEquationFlatSize(equation); uint64_t equationFlatIndex = 0; @@ -1952,42 +1866,6 @@ bool IDAInstance::idaSetLineSearchOffIC() { return retVal == IDA_SUCCESS; } -void IDAInstance::getWritingEquation( - Variable variable, const std::vector &variableIndices, - Equation &equation, std::vector &equationIndices) const { - bool found = false; - uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); - - for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { - Variable writtenVariable = getWrittenVariable(eq); - - if (writtenVariable == variable) { - std::vector writingEquationIndices; - getEquationBeginIndices(eq, writingEquationIndices); - - std::vector writtenVariableIndices( - getVariableRank(writtenVariable)); - - AccessFunction writeAccessFunction = getWriteAccessFunction(eq); - - do { - writeAccessFunction(writingEquationIndices.data(), - writtenVariableIndices.data()); - - if (writtenVariableIndices == variableIndices) { - assert(!found && "Multiple equations writing to the same variable"); - found = true; - equation = eq; - equationIndices = writingEquationIndices; - } - } while ( - advanceEquationIndices(writingEquationIndices, equationRanges[eq])); - } - } - - assert(found && "Writing equation not found"); -} - void IDAInstance::printVariablesVector(N_Vector variables) const { realtype *data = N_VGetArrayPointer(variables); uint64_t numOfArrayVariables = getNumOfArrayVariables(); @@ -2030,25 +1908,18 @@ void IDAInstance::printDerivativesVector(N_Vector derivatives) const { void IDAInstance::printResidualsVector(N_Vector residuals) const { realtype *data = N_VGetArrayPointer(residuals); - uint64_t numOfArrayVariables = getNumOfArrayVariables(); + uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); - for (Variable var = 0; var < numOfArrayVariables; ++var) { - std::vector variableIndices; - getVariableBeginIndices(var, variableIndices); + for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { + std::vector equationIndices; + getEquationBeginIndices(eq, equationIndices); do { - Equation eq; - std::vector equationIndices; - getWritingEquation(var, variableIndices, eq, equationIndices); - std::cerr << "eq " << eq << " "; printIndices(equationIndices); - std::cerr << " (writing to var " << var; - printIndices(variableIndices); - std::cerr << ")" << "\t" << std::fixed << std::setprecision(9) << *data - << "\n"; + std::cerr << "\t" << std::fixed << std::setprecision(9) << *data << "\n"; ++data; - } while (advanceVariableIndices(variableIndices, variablesDimensions[var])); + } while (advanceEquationIndices(equationIndices, equationRanges[eq])); } } @@ -2088,29 +1959,23 @@ void IDAInstance::printJacobianMatrix(SUNMatrix jacobianMatrix) const { std::cerr << std::endl; - // Print the rows containing the values. + // Print the matrix. + uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); uint64_t rowFlatIndex = 0; - for (Variable eqVar = 0; eqVar < numOfArrayVariables; ++eqVar) { - std::vector eqVarIndices; - getVariableBeginIndices(eqVar, eqVarIndices); + for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { + std::vector equationIndices; + getEquationBeginIndices(eq, equationIndices); do { - Equation eq; - std::vector equationIndices; - getWritingEquation(eqVar, eqVarIndices, eq, equationIndices); - std::cerr << "eq " << eq << " "; printIndices(equationIndices); - std::cerr << " (writing to var " << eqVar << " "; - printIndices(eqVarIndices); - std::cerr << ")"; uint64_t columnFlatIndex = 0; - for (Variable indVar = 0; indVar < numOfArrayVariables; ++indVar) { - std::vector indVarIndices; - getVariableBeginIndices(indVar, indVarIndices); + for (Variable var = 0; var < numOfArrayVariables; ++var) { + std::vector varIndices; + getVariableBeginIndices(var, varIndices); do { auto value = getCellFromSparseMatrix(jacobianMatrix, rowFlatIndex, @@ -2118,13 +1983,12 @@ void IDAInstance::printJacobianMatrix(SUNMatrix jacobianMatrix) const { std::cerr << "\t" << std::fixed << std::setprecision(9) << value; columnFlatIndex++; - } while ( - advanceVariableIndices(indVarIndices, variablesDimensions[indVar])); + } while (advanceVariableIndices(varIndices, variablesDimensions[var])); } std::cerr << std::endl; rowFlatIndex++; - } while (advanceVariableIndices(eqVarIndices, variablesDimensions[eqVar])); + } while (advanceEquationIndices(equationIndices, equationRanges[eq])); } } } // namespace marco::runtime::sundials::ida @@ -2262,17 +2126,13 @@ RUNTIME_FUNC_DEF(idaAddVariableAccess, void, PTR(void), uint64_t, uint64_t, // idaAddEquation static uint64_t idaAddEquation_i64(void *instance, int64_t *ranges, - uint64_t rank, uint64_t writtenVariable, - void *writeAccessFunction, - void *stringRepresentation) { + uint64_t rank, void *stringRepresentation) { return static_cast(instance)->addEquation( - ranges, rank, writtenVariable, - reinterpret_cast(writeAccessFunction), - static_cast(stringRepresentation)); + ranges, rank, static_cast(stringRepresentation)); } RUNTIME_FUNC_DEF(idaAddEquation, uint64_t, PTR(void), PTR(int64_t), uint64_t, - uint64_t, PTR(void), PTR(void)) + PTR(void)) //===---------------------------------------------------------------------===// // idaSetResidual diff --git a/lib/Solvers/KINSOL/Instance.cpp b/lib/Solvers/KINSOL/Instance.cpp index 5f33d3ef..2f30ff2d 100644 --- a/lib/Solvers/KINSOL/Instance.cpp +++ b/lib/Solvers/KINSOL/Instance.cpp @@ -23,8 +23,9 @@ using namespace marco::runtime::sundials::kinsol; namespace marco::runtime::sundials::kinsol { KINSOLInstance::KINSOLInstance() { - // Initially there is no variable in the instance. + // Initially there is are no variables or equations in the instance. variableOffsets.push_back(0); + equationOffsets.push_back(0); if (marco::runtime::simulation::getOptions().debug) { std::cerr << "[KINSOL] Instance created" << std::endl; @@ -114,8 +115,6 @@ Variable KINSOLInstance::addVariable(uint64_t rank, const uint64_t *dimensions, Equation KINSOLInstance::addEquation(const int64_t *ranges, uint64_t equationRank, - Variable writtenVariable, - AccessFunction writeAccess, const char *stringRepresentation) { if (marco::runtime::simulation::getOptions().debug) { std::cerr << "[KINSOL] Adding equation"; @@ -129,17 +128,18 @@ Equation KINSOLInstance::addEquation(const int64_t *ranges, // Add the start and end dimensions of the current equation. MultidimensionalRange eqRanges = {}; + uint64_t flatSize = 1; for (size_t i = 0, e = equationRank * 2; i < e; i += 2) { int64_t begin = ranges[i]; int64_t end = ranges[i + 1]; eqRanges.push_back({begin, end}); + flatSize *= end - begin; } equationRanges.push_back(eqRanges); - - // Add the write access. - writeAccesses.emplace_back(writtenVariable, writeAccess); + size_t offset = equationOffsets.back(); + equationOffsets.push_back(offset + flatSize); // Return the index of the equation. Equation id = getNumOfVectorizedEquations() - 1; @@ -159,9 +159,6 @@ Equation KINSOLInstance::addEquation(const int64_t *ranges, } std::cerr << "]" << std::endl; - std::cerr << " - Written variable ID: " << writtenVariable << std::endl; - std::cerr << " - Write access function address: " - << reinterpret_cast(writeAccess) << std::endl; } return id; @@ -335,29 +332,6 @@ bool KINSOLInstance::initialize() { N_VGetArrayPointer(residualScaleVector)[i] = 1; } - // Determine the order in which the equations must be processed when - // computing residuals and jacobians. - assert(getNumOfVectorizedEquations() == writeAccesses.size()); - equationsProcessingOrder.resize(getNumOfVectorizedEquations()); - - for (size_t i = 0, e = getNumOfVectorizedEquations(); i < e; ++i) { - equationsProcessingOrder[i] = i; - } - - if (marco::runtime::simulation::getOptions().debug) { - std::cerr << "[KINSOL] Equations processing order: ["; - - for (size_t i = 0, e = equationsProcessingOrder.size(); i < e; ++i) { - if (i != 0) { - std::cerr << ", "; - } - - std::cerr << equationsProcessingOrder[i]; - } - - std::cerr << "]" << std::endl; - } - // Check that all the residual functions have been set. assert(residualFunctions.size() == getNumOfVectorizedEquations()); @@ -403,35 +377,19 @@ bool KINSOLInstance::initialize() { jacobianMatrixData.resize(scalarEquationsNumber); - for (Equation eq : equationsProcessingOrder) { + uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); + + for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { std::vector equationIndices; getEquationBeginIndices(eq, equationIndices); - Variable writtenVariable = getWrittenVariable(eq); - - uint64_t writtenVariableRank = getVariableRank(writtenVariable); - uint64_t writtenVariableArrayOffset = variableOffsets[writtenVariable]; - do { - std::vector writtenVariableIndices; - writtenVariableIndices.resize(writtenVariableRank, 0); - - AccessFunction writeAccessFunction = getWriteAccessFunction(eq); - - writeAccessFunction(equationIndices.data(), - writtenVariableIndices.data()); - - if (marco::runtime::simulation::getOptions().debug) { - std::cerr << " Variable indices: "; - printIndices(writtenVariableIndices); - std::cerr << std::endl; - } + uint64_t equationArrayOffset = equationOffsets[eq]; - uint64_t equationScalarVariableOffset = getVariableFlatIndex( - variablesDimensions[writtenVariable], writtenVariableIndices); + uint64_t equationScalarOffset = + getEquationFlatIndex(equationIndices, equationRanges[eq]); - uint64_t scalarEquationIndex = - writtenVariableArrayOffset + equationScalarVariableOffset; + uint64_t scalarEquationIndex = equationArrayOffset + equationScalarOffset; // Compute the column indexes that may be non-zeros. std::vector jacobianColumns = @@ -445,10 +403,6 @@ bool KINSOLInstance::initialize() { printIndices(equationIndices); std::cerr << std::endl; - std::cerr << " Variable indices: "; - printIndices(writtenVariableIndices); - std::cerr << std::endl; - std::cerr << " Scalar equation index: " << scalarEquationIndex << std::endl; @@ -578,28 +532,12 @@ int KINSOLInstance::residualFunction(N_Vector variables, N_Vector residuals, uint64_t equationRank = instance->getEquationRank(eq); assert(equationIndices.size() == equationRank); - Variable writtenVariable = instance->getWrittenVariable(eq); - - uint64_t writtenVariableArrayOffset = - instance->variableOffsets[writtenVariable]; + uint64_t equationArrayOffset = instance->equationOffsets[eq]; - uint64_t writtenVariableRank = - instance->getVariableRank(writtenVariable); + uint64_t equationScalarOffset = + getEquationFlatIndex(equationIndices, instance->equationRanges[eq]); - std::vector writtenVariableIndices(writtenVariableRank, 0); - - AccessFunction writeAccessFunction = - instance->getWriteAccessFunction(eq); - - writeAccessFunction(equationIndices.data(), - writtenVariableIndices.data()); - - uint64_t writtenVariableScalarOffset = - getVariableFlatIndex(instance->variablesDimensions[writtenVariable], - writtenVariableIndices); - - uint64_t offset = - writtenVariableArrayOffset + writtenVariableScalarOffset; + uint64_t offset = equationArrayOffset + equationScalarOffset; auto residualFn = instance->residualFunctions[eq]; auto *eqIndicesPtr = equationIndices.data(); @@ -639,29 +577,13 @@ int KINSOLInstance::jacobianMatrix(N_Vector variables, N_Vector residuals, instance->equationsParallelIteration( [&](Equation eq, const std::vector &equationIndices, const JacobianSeedsMap &jacobianSeedsMap) { - Variable writtenVariable = instance->getWrittenVariable(eq); - - uint64_t writtenVariableArrayOffset = - instance->variableOffsets[writtenVariable]; - - uint64_t writtenVariableRank = - instance->getVariableRank(writtenVariable); - - std::vector writtenVariableIndices; - writtenVariableIndices.resize(writtenVariableRank, 0); - - AccessFunction writeAccessFunction = - instance->getWriteAccessFunction(eq); - - writeAccessFunction(equationIndices.data(), - writtenVariableIndices.data()); + uint64_t equationArrayOffset = instance->equationOffsets[eq]; - uint64_t writtenVariableScalarOffset = - getVariableFlatIndex(instance->variablesDimensions[writtenVariable], - writtenVariableIndices); + uint64_t equationScalarOffset = + getEquationFlatIndex(equationIndices, instance->equationRanges[eq]); uint64_t scalarEquationIndex = - writtenVariableArrayOffset + writtenVariableScalarOffset; + equationArrayOffset + equationScalarOffset; assert(scalarEquationIndex < instance->getNumOfScalarEquations()); @@ -782,14 +704,6 @@ uint64_t KINSOLInstance::getEquationFlatSize(Equation equation) const { return result; } -Variable KINSOLInstance::getWrittenVariable(Equation equation) const { - return writeAccesses[equation].first; -} - -AccessFunction KINSOLInstance::getWriteAccessFunction(Equation equation) const { - return writeAccesses[equation].second; -} - uint64_t KINSOLInstance::getVariableRank(Variable variable) const { return variablesDimensions[variable].rank(); } @@ -930,7 +844,7 @@ void KINSOLInstance::computeThreadChunks() { uint64_t processedEquations = 0; while (processedEquations < numOfVectorizedEquations) { - Equation equation = equationsProcessingOrder[processedEquations]; + Equation equation = processedEquations; uint64_t equationFlatSize = getEquationFlatSize(equation); uint64_t equationFlatIndex = 0; @@ -1274,42 +1188,6 @@ bool KINSOLInstance::kinsolSetJacobianFunction() { return retVal == KIN_SUCCESS; } -void KINSOLInstance::getWritingEquation( - Variable variable, const std::vector &variableIndices, - Equation &equation, std::vector &equationIndices) const { - bool found = false; - uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); - - for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { - Variable writtenVariable = getWrittenVariable(eq); - - if (writtenVariable == variable) { - std::vector writingEquationIndices; - getEquationBeginIndices(eq, writingEquationIndices); - - std::vector writtenVariableIndices( - getVariableRank(writtenVariable)); - - AccessFunction writeAccessFunction = getWriteAccessFunction(eq); - - do { - writeAccessFunction(writingEquationIndices.data(), - writtenVariableIndices.data()); - - if (writtenVariableIndices == variableIndices) { - assert(!found && "Multiple equations writing to the same variable"); - found = true; - equation = eq; - equationIndices = writingEquationIndices; - } - } while ( - advanceEquationIndices(writingEquationIndices, equationRanges[eq])); - } - } - - assert(found && "Writing equation not found"); -} - void KINSOLInstance::printVariablesVector(N_Vector variables) const { realtype *data = N_VGetArrayPointer(variables); uint64_t numOfArrayVariables = getNumOfArrayVariables(); @@ -1330,25 +1208,18 @@ void KINSOLInstance::printVariablesVector(N_Vector variables) const { void KINSOLInstance::printResidualsVector(N_Vector residuals) const { realtype *data = N_VGetArrayPointer(residuals); - uint64_t numOfArrayVariables = getNumOfArrayVariables(); + uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); - for (Variable var = 0; var < numOfArrayVariables; ++var) { - std::vector variableIndices; - getVariableBeginIndices(var, variableIndices); + for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { + std::vector equationIndices; + getEquationBeginIndices(eq, equationIndices); do { - Equation eq; - std::vector equationIndices; - getWritingEquation(var, variableIndices, eq, equationIndices); - std::cerr << "eq " << eq << " "; printIndices(equationIndices); - std::cerr << " (writing to var " << var; - printIndices(variableIndices); - std::cerr << ")" << "\t" << std::fixed << std::setprecision(9) << *data - << "\n"; + std::cerr << "\t" << std::fixed << std::setprecision(9) << *data << "\n"; ++data; - } while (advanceVariableIndices(variableIndices, variablesDimensions[var])); + } while (advanceEquationIndices(equationIndices, equationRanges[eq])); } } @@ -1368,29 +1239,23 @@ void KINSOLInstance::printJacobianMatrix(SUNMatrix jacobianMatrix) const { std::cerr << std::endl; - // Print the rows containing the values. + // Print the matrix. + uint64_t numOfVectorizedEquations = getNumOfVectorizedEquations(); uint64_t rowFlatIndex = 0; - for (Variable eqVar = 0; eqVar < numOfArrayVariables; ++eqVar) { - std::vector eqVarIndices; - getVariableBeginIndices(eqVar, eqVarIndices); + for (Equation eq = 0; eq < numOfVectorizedEquations; ++eq) { + std::vector equationIndices; + getEquationBeginIndices(eq, equationIndices); do { - Equation eq; - std::vector equationIndices; - getWritingEquation(eqVar, eqVarIndices, eq, equationIndices); - std::cerr << "eq " << eq << " "; printIndices(equationIndices); - std::cerr << " (writing to var " << eqVar << " "; - printIndices(eqVarIndices); - std::cerr << ")"; uint64_t columnFlatIndex = 0; - for (Variable indVar = 0; indVar < numOfArrayVariables; ++indVar) { - std::vector indVarIndices; - getVariableBeginIndices(indVar, indVarIndices); + for (Variable var = 0; var < numOfArrayVariables; ++var) { + std::vector varIndices; + getVariableBeginIndices(var, varIndices); do { auto value = getCellFromSparseMatrix(jacobianMatrix, rowFlatIndex, @@ -1398,13 +1263,12 @@ void KINSOLInstance::printJacobianMatrix(SUNMatrix jacobianMatrix) const { std::cerr << "\t" << std::fixed << std::setprecision(9) << value; columnFlatIndex++; - } while ( - advanceVariableIndices(indVarIndices, variablesDimensions[indVar])); + } while (advanceVariableIndices(varIndices, variablesDimensions[var])); } std::cerr << std::endl; rowFlatIndex++; - } while (advanceVariableIndices(eqVarIndices, variablesDimensions[eqVar])); + } while (advanceEquationIndices(equationIndices, equationRanges[eq])); } } } // namespace marco::runtime::sundials::kinsol @@ -1476,17 +1340,14 @@ RUNTIME_FUNC_DEF(kinsolAddVariableAccess, void, PTR(void), uint64_t, uint64_t, // kinsolAddEquation static uint64_t kinsolAddEquation_i64(void *instance, int64_t *ranges, - uint64_t rank, uint64_t writtenVariable, - void *writeAccessFunction, + uint64_t rank, void *stringRepresentation) { return static_cast(instance)->addEquation( - ranges, rank, writtenVariable, - reinterpret_cast(writeAccessFunction), - static_cast(stringRepresentation)); + ranges, rank, static_cast(stringRepresentation)); } RUNTIME_FUNC_DEF(kinsolAddEquation, uint64_t, PTR(void), PTR(int64_t), uint64_t, - uint64_t, PTR(void), PTR(void)) + PTR(void)) //===---------------------------------------------------------------------===// // kinsolSetResidual