improved symbolic linear equation solving (via Jacobi) a bit

resources/3rdparty/cudd-3.0.0/cudd/cuddSat.c

@@ -843,11 +843,12 @@ Cudd_EqualSupNormRel(
     /* Check terminal cases. */
     if (f == g) return(1);
     if (Cudd_IsConstant(f) && Cudd_IsConstant(g)) {
-        if (ddAbs((cuddV(f) - cuddV(g))/cuddV(f)) < tolerance) {
+        CUDD_VALUE_TYPE absDiff = ddAbs((cuddV(f) - cuddV(g)));
+        if (absDiff/cuddV(f) < tolerance || absDiff < Cudd_ReadEpsilon(dd)) {
             return(1);
         } else {
             if (pr>0) {
-                (void) fprintf(dd->out,"Offending nodes:\n");
+                (void) fprintf(dd->out,"Offending nodes (wrt. precision %0.30f) with diff %0.30f:\n", Cudd_ReadEpsilon(dd), absDiff);
                 (void) fprintf(dd->out,
                                "f: address = %p\t value = %40.30f\n",
                                (void *) f, cuddV(f));

src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp

@@ -175,7 +175,7 @@ namespace storm {
 
                         // Solve the equation system.
                         std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
-                        storm::dd::Add<DdType, ValueType> result = solver->solveEquations(model.getManager().getConstant(0.5) * maybeStatesAdd, subvector);
+                        storm::dd::Add<DdType, ValueType> result = solver->solveEquations(model.getManager().getConstant(0.0) * maybeStatesAdd, subvector);
                         
                         return infinityStates.ite(model.getManager().getConstant(storm::utility::infinity<ValueType>()), result);
                     } else {

src/storm/settings/modules/CuddSettings.cpp

@@ -20,7 +20,7 @@ namespace storm {
             const std::string CuddSettings::reorderOptionName = "reorder";
             
             CuddSettings::CuddSettings() : ModuleSettings(moduleName) {
-                this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, true, "Sets the precision used by Cudd.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision up to which to constants are considered to be different.").setDefaultValueDouble(1e-15).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+                this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, true, "Sets the precision used by Cudd.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision up to which to constants are considered to be different.").setDefaultValueDouble(1e-16).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, maximalMemoryOptionName, true, "Sets the upper bound of memory available to Cudd in MB.").addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("value", "The memory available to Cudd (0 means unlimited).").setDefaultValueUnsignedInteger(4096).build()).build());
                 

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -110,7 +110,7 @@ namespace storm {
         
         template<typename ValueType>
         bool NativeLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
-            if(!this->cachedRowVector) {
+            if (!this->cachedRowVector) {
                 this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
             }
             
@@ -180,7 +180,7 @@ namespace storm {
                     std::swap(x, *currentX);
                 }
                 
-                if(!this->isCachingEnabled()) {
+                if (!this->isCachingEnabled()) {
                     clearCache();
                 }
                 
@@ -208,7 +208,7 @@ namespace storm {
                     result.swap(*this->cachedRowVector);
                 }
                 
-                if(!this->isCachingEnabled()) {
+                if (!this->isCachingEnabled()) {
                     clearCache();
                 }
             }

src/storm/solver/SymbolicLinearEquationSolver.cpp

@@ -37,7 +37,10 @@ namespace storm {
             
             storm::dd::Add<DdType, ValueType> lu = diagonal.ite(this->A.getDdManager().template getAddZero<ValueType>(), this->A);
             storm::dd::Add<DdType> diagonalAdd = diagonal.template toAdd<ValueType>();
-            storm::dd::Add<DdType, ValueType> dinv = diagonalAdd / (diagonalAdd * this->A);
+            storm::dd::Add<DdType, ValueType> diag = diagonalAdd.multiplyMatrix(this->A, this->columnMetaVariables);
+            
+            storm::dd::Add<DdType, ValueType> scaledLu = lu / diag;
+            storm::dd::Add<DdType, ValueType> scaledB = b / diag;
             
             // Set up additional environment variables.
             storm::dd::Add<DdType, ValueType> xCopy = x;
@@ -46,21 +49,23 @@ namespace storm {
             
             while (!converged && iterationCount < maximalNumberOfIterations) {
                 storm::dd::Add<DdType, ValueType> xCopyAsColumn = xCopy.swapVariables(this->rowColumnMetaVariablePairs);
-
-                storm::dd::Add<DdType, ValueType> tmp = lu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
-                tmp = b - tmp;
-                tmp = tmp.swapVariables(this->rowColumnMetaVariablePairs);
-                tmp = dinv.multiplyMatrix(tmp, this->columnMetaVariables);
+                storm::dd::Add<DdType, ValueType> tmp = scaledB - scaledLu.multiplyMatrix(xCopyAsColumn, this->columnMetaVariables);
                 
                 // Now check if the process already converged within our precision.
-                converged = xCopy.equalModuloPrecision(tmp, precision, relative);
+                converged = tmp.equalModuloPrecision(xCopy, precision, relative);
 
                 xCopy = tmp;
                 
                 // Increase iteration count so we can abort if convergence is too slow.
                 ++iterationCount;
             }
-                        
+            
+            if (converged) {
+                STORM_LOG_TRACE("Iterative solver converged in " << iterationCount << " iterations.");
+            } else {
+                STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterstions.");
+            }
+            
             return xCopy;
         }
         

src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp

@@ -43,8 +43,6 @@ namespace storm {
                 tmp += b;
                 
                 if (minimize) {
-                    // This is a hack and only here because of the lack of a suitable minAbstract/maxAbstract function
-                    // that can properly deal with a restriction of the choices.
                     tmp += illegalMaskAdd;
                     tmp = tmp.minAbstract(this->choiceVariables);
                 } else {
@@ -54,11 +52,14 @@ namespace storm {
                 // Now check if the process already converged within our precision.
                 converged = xCopy.equalModuloPrecision(tmp, precision, relative);
                 
-                // If the method did not converge yet, we prepare the x vector for the next iteration.
-                if (!converged) {
-                    xCopy = tmp;
-                }
+                xCopy = tmp;
                 
+                if (converged) {
+                    STORM_LOG_TRACE("Iterative solver converged in " << iterations << " iterations.");
+                } else {
+                    STORM_LOG_WARN("Iterative solver did not converge in " << iterations << " iterstions.");
+                }
+
                 ++iterations;
             }