OMCompiler/Compiler/NBackEnd/Util/NBASSC.mo
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * This file is part of OpenModelica. | ||
| 3 | * | ||
| 4 | * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), | ||
| 5 | * c/o Linköpings universitet, Department of Computer and Information Science, | ||
| 6 | * SE-58183 Linköping, Sweden. | ||
| 7 | * | ||
| 8 | * All rights reserved. | ||
| 9 | * | ||
| 10 | * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF AGPL VERSION 3 LICENSE OR | ||
| 11 | * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. | ||
| 12 | * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES | ||
| 13 | * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GNU AGPL | ||
| 14 | * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. | ||
| 15 | * | ||
| 16 | * The OpenModelica software and the OSMC (Open Source Modelica Consortium) | ||
| 17 | * Public License (OSMC-PL) are obtained from OSMC, either from the above | ||
| 18 | * address, from the URLs: | ||
| 19 | * http://www.openmodelica.org or | ||
| 20 | * https://github.com/OpenModelica/ or | ||
| 21 | * http://www.ida.liu.se/projects/OpenModelica, | ||
| 22 | * and in the OpenModelica distribution. | ||
| 23 | * | ||
| 24 | * GNU AGPL version 3 is obtained from: | ||
| 25 | * https://www.gnu.org/licenses/licenses.html#GPL | ||
| 26 | * | ||
| 27 | * This program is distributed WITHOUT ANY WARRANTY; without | ||
| 28 | * even the implied warranty of MERCHANTABILITY or FITNESS | ||
| 29 | * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH | ||
| 30 | * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL. | ||
| 31 | * | ||
| 32 | * See the full OSMC Public License conditions for more details. | ||
| 33 | * | ||
| 34 | */ | ||
| 35 | |||
| 36 | encapsulated package NBASSC | ||
| 37 | "file: NBASSC.mo | ||
| 38 | package: NBASSC | ||
| 39 | description: This file contains the functions which will perform analytical to structural singularity conversion. | ||
| 40 | " | ||
| 41 | public import DAE; | ||
| 42 | public import ExpressionBasics; | ||
| 43 | public import ExpressionDump; | ||
| 44 | public import ExpressionSimplify; | ||
| 45 | |||
| 46 | protected | ||
| 47 | // NF imports | ||
| 48 | import ComponentRef = NFComponentRef; | ||
| 49 | import Expression = NFExpression; | ||
| 50 | import NFFlatten.FunctionTreeImpl; | ||
| 51 | import NFFunction.Function; | ||
| 52 | import Operator = NFOperator; | ||
| 53 | import SimplifyExp = NFSimplifyExp; | ||
| 54 | import Type = NFType; | ||
| 55 | |||
| 56 | // Backend imports | ||
| 57 | import Differentiate = NBDifferentiate; | ||
| 58 | import NBDifferentiate.{DifferentiationType, DifferentiationArguments}; | ||
| 59 | import NBEquation.{Equation, EquationAttributes, EquationKind, EquationPointer, EquationPointers, Iterator}; | ||
| 60 | import Replacements = NBReplacements; | ||
| 61 | import Solve = NBSolve; | ||
| 62 | import NBSolve.Status; | ||
| 63 | import BVariable = NBVariable; | ||
| 64 | import NBVariable.{VariablePointers, VariablePointer, VarData}; | ||
| 65 | |||
| 66 | // Util imports | ||
| 67 | import UnorderedMap; | ||
| 68 | import UnorderedSet; | ||
| 69 | |||
| 70 | public | ||
| 71 | function main | ||
| 72 | "Main function that resolves cyclic alias sets using Bareiss elimination." | ||
| 73 | input list<Pointer<Equation>> eqns; | ||
| 74 | input list<ComponentRef> vars; | ||
| 75 | input Pointer<Integer> index; | ||
| 76 | output list<Pointer<Equation>> resolved_eqns = {}; | ||
| 77 | protected | ||
| 78 | array<list<Integer>> indices, values; | ||
| 79 | Integer num_crefs, num_eqns, num_nonzero_val; | ||
| 80 | array<Integer> op_modes, op_val1, op_val2, op_val3, op_val4; | ||
| 81 | Integer num_op, count_zero_row; | ||
| 82 | UnorderedMap<EquationPointer, Expression> lhs_map; | ||
| 83 | array<Expression> lhs_array; | ||
| 84 | Boolean singular; | ||
| 85 | algorithm | ||
| 86 | 5 | (indices, values, num_crefs, num_eqns, num_nonzero_val, lhs_map) := buildSparseRepresentation(eqns, vars); | |
| 87 | 5 | setMatrix(num_crefs, num_eqns, num_nonzero_val, indices, values); | |
| 88 | 5 | (indices, values) := performBareissElimination(indices, values); | |
| 89 | 5 | (num_op, op_modes, op_val1, op_val2, op_val3, op_val4, lhs_array) := applyRecordedOperations(lhs_map); | |
| 90 | // check if the matrix is singular | ||
| 91 | 5 | (singular, count_zero_row) := checkSingularity(indices, num_eqns); | |
| 92 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3 times.
|
5 | if singular then |
| 93 | 2 | tracebackZeroRows(eqns, num_eqns, count_zero_row, num_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 94 | else | ||
| 95 | // create new equations from the transformed matrix | ||
| 96 | 3 | resolved_eqns := createEquations(vars, index, indices, values, num_eqns, lhs_array); | |
| 97 | end if; | ||
| 98 | 3 | freeMatrix(); | |
| 99 | end main; | ||
| 100 | |||
| 101 | function buildSparseRepresentation | ||
| 102 | "Builds the sparse matrix representation from the equation system." | ||
| 103 | input list<Pointer<Equation>> eqns; | ||
| 104 | input list<ComponentRef> vars; | ||
| 105 | output array<list<Integer>> indices, values; | ||
| 106 | output Integer num_crefs, num_eqns, num_nonzero_val; | ||
| 107 | output UnorderedMap<EquationPointer, Expression> lhs_map = UnorderedMap.new<Expression>(Equation.hash, Equation.isEqualPtr); | ||
| 108 | protected | ||
| 109 | Boolean b = true; | ||
| 110 | list<ComponentRef> cref_lst; | ||
| 111 | Expression res, diff_res, expr; | ||
| 112 | DifferentiationArguments args; | ||
| 113 | Integer diff_res_int, eqn_index, var_index, var_index1, var_index2; | ||
| 114 | Tuple_Id id; | ||
| 115 | UnorderedMap<Tuple_Id,Integer> diffs = UnorderedMap.new<Integer>(Tuple_Id.hash, Tuple_Id.isEqual); | ||
| 116 | UnorderedSet<EquationPointer> int_eqns = UnorderedSet.new(Equation.hash, Equation.isEqualPtr); | ||
| 117 | UnorderedSet<ComponentRef> int_crefs = UnorderedSet.new(ComponentRef.hash, ComponentRef.isEqual); | ||
| 118 | UnorderedMap<EquationPointer,CrefLst> rows = UnorderedMap.new<CrefLst>(Equation.hash, Equation.isEqualPtr); | ||
| 119 | UnorderedMap<ComponentRef, Expression> replacements = UnorderedMap.new<Expression>(ComponentRef.hash, ComponentRef.isEqual); | ||
| 120 | list<Integer> lst_enum; | ||
| 121 | list<EquationPointer> lst_eqns; | ||
| 122 | list<ComponentRef> crefs_rows; | ||
| 123 | UnorderedMap<EquationPointer, Integer> enum_eqns; | ||
| 124 | UnorderedMap<ComponentRef, Integer> enum_crefs; | ||
| 125 | algorithm | ||
| 126 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 5 times.
|
29 | for eq_ptr in eqns loop |
| 127 | // find all crefs of vars in eqn | ||
| 128 | 24 | cref_lst := Equation.collectCrefs(Pointer.access(eq_ptr), function Equation.collectFromMap(check_map = UnorderedMap.fromLists(vars, vars, ComponentRef.hash, ComponentRef.isEqual))); | |
| 129 | 24 | res := Equation.getResidualExp(Pointer.access(eq_ptr)); | |
| 130 | 24 | args := Differentiate.DifferentiationArguments.default(NBDifferentiate.DifferentiationType.SIMPLE); | |
| 131 |
2/2✓ Branch 1 taken 48 times.
✓ Branch 2 taken 24 times.
|
72 | for cr in cref_lst loop |
| 132 | 48 | args.diffCref := cr; | |
| 133 | // differentiate eqn for cref | ||
| 134 | 48 | diff_res := SimplifyExp.simplify(Differentiate.differentiateExpression(res, args)); | |
| 135 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
|
48 | if Type.isReal(Expression.typeOf(diff_res)) then // if diff is real; check if its an integer of type real |
| 136 | ✗ | diff_res_int := realInt(Expression.realValue(diff_res)); | |
| 137 | ✗ | id := TUPLE_ID(eq_ptr, cr); | |
| 138 | elseif Type.isInteger(Expression.typeOf(diff_res)) then // if diff is integer | ||
| 139 | 48 | diff_res_int := Expression.integerValue(diff_res); | |
| 140 | 48 | id := TUPLE_ID(eq_ptr, cr); | |
| 141 | 48 | UnorderedMap.add(id, diff_res_int, diffs); | |
| 142 | else | ||
| 143 | // eqn is not part of linear system | ||
| 144 | b := false; | ||
| 145 | break; | ||
| 146 | end if; | ||
| 147 | end for; | ||
| 148 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | if b then |
| 149 | 24 | UnorderedSet.add(eq_ptr, int_eqns); | |
| 150 |
2/2✓ Branch 0 taken 48 times.
✓ Branch 1 taken 24 times.
|
72 | for cr in cref_lst loop |
| 151 | 48 | UnorderedSet.add(cr, int_crefs); | |
| 152 | 48 | UnorderedMap.add(cr, Expression.makeZero(Equation.getType(Pointer.access(eq_ptr))), replacements); | |
| 153 | end for; | ||
| 154 | 24 | UnorderedMap.add(eq_ptr, cref_lst, rows); | |
| 155 | 24 | expr := SimplifyExp.simplify(Expression.map(res, function Replacements.applySimpleExp(replacements = replacements))); | |
| 156 | 24 | UnorderedMap.add(eq_ptr, Expression.negate(expr), lhs_map); | |
| 157 | end if; | ||
| 158 | end for; | ||
| 159 | // enumerate sets | ||
| 160 | 5 | lst_enum := List.intRange(UnorderedSet.size(int_eqns)); | |
| 161 | 5 | enum_eqns := UnorderedMap.fromLists(eqns, lst_enum, Equation.hash, Equation.isEqualPtr); | |
| 162 | 5 | lst_enum := List.intRange(UnorderedSet.size(int_crefs)); | |
| 163 | 5 | enum_crefs := UnorderedMap.fromLists(vars, lst_enum, ComponentRef.hash, ComponentRef.isEqual); | |
| 164 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 165 | ✗ | print("Variable-to-column mapping:\n"+UnorderedMap.toString(enum_crefs, ComponentRef.toString, intString)+"\n"); | |
| 166 | end if; | ||
| 167 | 5 | indices := arrayCreate(UnorderedSet.size(int_eqns), {}); | |
| 168 | 5 | values := arrayCreate(UnorderedSet.size(int_eqns), {}); | |
| 169 | // create matrix elements for sparse matrix | ||
| 170 | 5 | lst_eqns := UnorderedSet.toList(int_eqns); | |
| 171 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 5 times.
|
29 | for eq_ptr in lst_eqns loop |
| 172 | 24 | crefs_rows := UnorderedMap.getSafe(eq_ptr, rows, sourceInfo()); | |
| 173 | 24 | var_index1 := UnorderedMap.getSafe(listGet(crefs_rows,1), enum_crefs, sourceInfo()); | |
| 174 | 24 | var_index2 := UnorderedMap.getSafe(listGet(crefs_rows,2), enum_crefs, sourceInfo()); | |
| 175 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 12 times.
|
24 | if var_index1 > var_index2 then |
| 176 | 12 | crefs_rows := listReverse(crefs_rows); | |
| 177 | end if; | ||
| 178 |
2/2✓ Branch 1 taken 48 times.
✓ Branch 2 taken 24 times.
|
72 | for cr in listReverse(crefs_rows) loop |
| 179 | 48 | eqn_index := UnorderedMap.getSafe(eq_ptr, enum_eqns, sourceInfo()); | |
| 180 | 48 | var_index := UnorderedMap.getSafe(cr, enum_crefs, sourceInfo()); | |
| 181 | 96 | indices[eqn_index] := var_index :: indices[eqn_index]; | |
| 182 | 96 | values[eqn_index] := UnorderedMap.getSafe(TUPLE_ID(eq_ptr,cr), diffs, sourceInfo()) :: values[eqn_index]; | |
| 183 | end for; | ||
| 184 | end for; | ||
| 185 | // determine sparse matrix dimensions | ||
| 186 | 5 | num_crefs := UnorderedSet.size(int_crefs); | |
| 187 | 5 | num_eqns := UnorderedSet.size(int_eqns); | |
| 188 | 5 | num_nonzero_val := UnorderedMap.size(diffs); | |
| 189 | end buildSparseRepresentation; | ||
| 190 | |||
| 191 | function performBareissElimination | ||
| 192 | "Performs the Bareiss elimination procedure on the system matrix." | ||
| 193 | input output array<list<Integer>> indices, values; | ||
| 194 | algorithm | ||
| 195 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 196 | ✗ | print("Sparse matrix before applying the Bareiss algorithm:\n"); | |
| 197 | ✗ | printMatrix(); | |
| 198 | end if; | ||
| 199 | 5 | bareiss(); | |
| 200 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 201 | ✗ | print("Sparse matrix after applying the Bareiss algorithm:\n"); | |
| 202 | ✗ | printMatrix(); | |
| 203 | ✗ | print("\n"); | |
| 204 | end if; | ||
| 205 | 5 | indices := arrayCreate(arrayLength(indices),{}); | |
| 206 | 5 | values := arrayCreate(arrayLength(values),{}); | |
| 207 | 5 | getMatrix(indices,values); | |
| 208 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 209 | ✗ | print("List indices:\n"); | |
| 210 | ✗ | for i in 1:arrayLength(indices) loop | |
| 211 | ✗ | print(List.toString(indices[i], intString) + "\n"); | |
| 212 | end for; | ||
| 213 | ✗ | print("\nList values:"); | |
| 214 | ✗ | for i in 1:arrayLength(values) loop | |
| 215 | ✗ | print(List.toString(values[i], intString) + "\n"); | |
| 216 | end for; | ||
| 217 | ✗ | print("\n"); | |
| 218 | end if; | ||
| 219 | end performBareissElimination; | ||
| 220 | |||
| 221 | function applyRecordedOperations | ||
| 222 | "Applies the recorded Bareiss operations to the left-hand side expressions." | ||
| 223 | input UnorderedMap<EquationPointer, Expression> lhs_map; | ||
| 224 | output Integer num_op; | ||
| 225 | output array<Integer> op_modes, op_val1, op_val2, op_val3, op_val4; | ||
| 226 | output array<Expression> lhs_array; | ||
| 227 | protected | ||
| 228 | array<Integer> nop; | ||
| 229 | Integer mode; | ||
| 230 | algorithm | ||
| 231 | 5 | nop := arrayCreate(1,-1); | |
| 232 | 5 | num_op := getNumberOfOperations(nop); | |
| 233 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 234 | ✗ | print("Number of operations: "+intString(num_op)+"\n"); | |
| 235 | end if; | ||
| 236 | // allocate operation storage | ||
| 237 | 5 | op_modes := arrayCreate(num_op,-1); | |
| 238 | 5 | op_val1 := arrayCreate(num_op,-1); | |
| 239 | 5 | op_val2 := arrayCreate(num_op,-1); | |
| 240 | 5 | op_val3 := arrayCreate(num_op,-1); | |
| 241 | 5 | op_val4 := arrayCreate(num_op,-1); | |
| 242 | // retrieve all recorded Bareiss operations from the runtime | ||
| 243 | 5 | getOperations(op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 244 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 245 | ✗ | print("All operations:\n"); | |
| 246 | ✗ | print(Array.toString(op_modes, intString)+"\n"); | |
| 247 | ✗ | print(Array.toString(op_val1, intString)+"\n"); | |
| 248 | ✗ | print(Array.toString(op_val2, intString)+"\n"); | |
| 249 | ✗ | print(Array.toString(op_val3, intString)+"\n"); | |
| 250 | ✗ | print(Array.toString(op_val4, intString)+"\n"); | |
| 251 | end if; | ||
| 252 | // apply operations on left-hand side | ||
| 253 | 5 | lhs_array := listArray(UnorderedMap.valueList(lhs_map)); | |
| 254 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 255 | ✗ | print("\nlhs array: "+Array.toString(lhs_array, Expression.toString)+"\n"); | |
| 256 | end if; | ||
| 257 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
38 | for i in 1:num_op loop |
| 258 | 33 | mode := op_modes[i]; | |
| 259 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 33 times.
|
33 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 260 | ✗ | print("current op_mode: "+intString(mode)+"\n"); | |
| 261 | end if; | ||
| 262 | _:= match mode | ||
| 263 | local | ||
| 264 | Expression tmp_val, tmp_val1, tmp_val2; | ||
| 265 | Integer gcd; | ||
| 266 | case 0 algorithm // mode for pivot-update operation | ||
| 267 | 51 | tmp_val1 := Expression.MULTARY( | |
| 268 | arguments = {Expression.makeInteger(op_val2[i]), lhs_array[op_val3[i]+1]}, | ||
| 269 | inv_arguments = {}, | ||
| 270 | operator = Operator.makeMul(Type.REAL())); | ||
| 271 | 17 | tmp_val1:= SimplifyExp.simplify(tmp_val1); | |
| 272 | 51 | tmp_val2 := Expression.MULTARY( | |
| 273 | arguments = {Expression.makeInteger(op_val4[i]), lhs_array[op_val1[i]+1]}, | ||
| 274 | inv_arguments = {}, | ||
| 275 | operator = Operator.makeMul(Type.REAL())); | ||
| 276 | 17 | tmp_val2:= SimplifyExp.simplify(tmp_val2); | |
| 277 | 17 | lhs_array[op_val3[i]+1]:= Expression.MULTARY( | |
| 278 | arguments = {tmp_val1}, | ||
| 279 | inv_arguments = {tmp_val2}, | ||
| 280 | operator = Operator.makeAdd(Type.REAL())); | ||
| 281 | 17 | lhs_array[op_val3[i]+1] := SimplifyExp.simplify(lhs_array[op_val3[i]+1]); | |
| 282 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
|
17 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 283 | ✗ | print("case 0, updated lhs_array: "+Array.toString(lhs_array, Expression.toString)+"\n"); | |
| 284 | end if; | ||
| 285 | then lhs_array; | ||
| 286 | case 1 algorithm // mode for swap-rows operation | ||
| 287 | 10 | tmp_val := lhs_array[op_val1[i]+1]; | |
| 288 | 10 | lhs_array[op_val1[i]+1] := lhs_array[op_val2[i]+1]; | |
| 289 | 10 | lhs_array[op_val2[i]+1] := tmp_val; | |
| 290 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
|
10 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 291 | ✗ | print("case 1, updated lhs_array: "+Array.toString(lhs_array, Expression.toString)+"\n"); | |
| 292 | end if; | ||
| 293 | then lhs_array; | ||
| 294 | case 2 algorithm // mode for gcd operation | ||
| 295 | 6 | gcd := op_val2[i]; | |
| 296 | 18 | lhs_array[op_val1[i]+1] := Expression.MULTARY( | |
| 297 | arguments = {lhs_array[op_val1[i]+1]}, | ||
| 298 | inv_arguments = {Expression.makeInteger(gcd)}, | ||
| 299 | operator = Operator.makeMul(Type.REAL())); | ||
| 300 | 6 | lhs_array[op_val1[i]+1] := SimplifyExp.simplify(lhs_array[op_val1[i]+1]); | |
| 301 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
|
6 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 302 | ✗ | print("case 2, updated lhs_array: "+Array.toString(lhs_array, Expression.toString)+"\n"); | |
| 303 | end if; | ||
| 304 | then lhs_array; | ||
| 305 | else lhs_array; | ||
| 306 | end match; | ||
| 307 | end for; | ||
| 308 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
5 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 309 | ✗ | print("final lhs_array: "+Array.toString(lhs_array, Expression.toString)+"\n\n"); | |
| 310 | end if; | ||
| 311 | end applyRecordedOperations; | ||
| 312 | |||
| 313 | function checkSingularity | ||
| 314 | "Detects singular matrices by checking for zero rows." | ||
| 315 | input array<list<Integer>> indices; | ||
| 316 | input Integer num_eqns; | ||
| 317 | output Boolean singular = false; | ||
| 318 | output Integer count_zero_row = 0; | ||
| 319 | algorithm | ||
| 320 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
29 | for i in 1:num_eqns loop |
| 321 |
2/2✓ Branch 2 taken 2 times.
✓ Branch 3 taken 22 times.
|
24 | if listLength(indices[i]) == 0 then |
| 322 | 2 | count_zero_row := count_zero_row + 1; | |
| 323 | end if; | ||
| 324 | end for; | ||
| 325 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
|
5 | if count_zero_row > 0 then |
| 326 | singular := true; | ||
| 327 | end if; | ||
| 328 | end checkSingularity; | ||
| 329 | |||
| 330 | function tracebackZeroRows | ||
| 331 | "Reconstructs the symbolic expression of each zero row and generates a detailed error message for singular systems." | ||
| 332 | input list<Pointer<Equation>> eqns; | ||
| 333 | input Integer num_eqns, count_zero_row, num_op; | ||
| 334 | input array<Integer> op_modes, op_val1, op_val2, op_val3, op_val4; | ||
| 335 | protected | ||
| 336 | Integer current_zero_row; | ||
| 337 | DAE.Exp exp_dae; | ||
| 338 | String eq_str = "", str_all = ""; | ||
| 339 | algorithm | ||
| 340 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 341 | ✗ | print("Number of zero rows: "+intString(count_zero_row)+"\n"); | |
| 342 | end if; | ||
| 343 | // process each zero row generated during elimination | ||
| 344 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
4 | for zero_row in 1:count_zero_row loop |
| 345 | 2 | current_zero_row := num_eqns - zero_row; | |
| 346 | // reconstruct the linear combination that produced the zero row. | ||
| 347 | 2 | exp_dae := traceEquation(current_zero_row, num_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 348 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 349 | ✗ | print("Reconstructed zero row: "+ExpressionBasics.printExpStr(exp_dae)+"\n"); | |
| 350 | end if; | ||
| 351 | // collect equation information for the error message. | ||
| 352 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10 | for eq in 1:num_eqns loop |
| 353 | 8 | eq_str := eq_str + "("+intString(eq-1)+"): " + Expression.toString(Util.getOption(Equation.getLHS(Pointer.access(listGet(eqns, eq))))) + " = " + Expression.toString(Util.getOption(Equation.getRHS(Pointer.access(listGet(eqns, eq))))) + "\n"; | |
| 354 | end for; | ||
| 355 | 2 | str_all := str_all + "The zero row in ("+ intString(current_zero_row) +") was produced by the following calculation: " + ExpressionBasics.printExpStr(exp_dae) + " with \n" + eq_str + "\n"; | |
| 356 | end for; | ||
| 357 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 358 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because sparse matrix is singular.\n" + str_all}); | |
| 359 | ✗ | fail(); | |
| 360 | else | ||
| 361 | 2 | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because sparse matrix is singular, for more information please use -d=dumpASSC.\n"}); | |
| 362 | 2 | fail(); | |
| 363 | end if; | ||
| 364 | end tracebackZeroRows; | ||
| 365 | |||
| 366 | function traceEquation | ||
| 367 | "Recursively reconstructs the expression of a matrix row by tracing back elimination operations." | ||
| 368 | input Integer current_row, last_op; | ||
| 369 | input array<Integer> op_modes, op_val1, op_val2, op_val3, op_val4; | ||
| 370 | output DAE.Exp exp_dae; | ||
| 371 | protected | ||
| 372 | list<DAE.Exp> expressions; | ||
| 373 | Boolean found; | ||
| 374 | Integer pre_op, pre_mode, row1, row2, factor1, factor2, gcd; | ||
| 375 | DAE.Exp pivot_exp, update_exp; | ||
| 376 | algorithm | ||
| 377 | 18 | (found, pre_op, pre_mode, row1, row2, factor1, factor2, gcd) := findLastOperation(current_row, last_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 378 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 10 times.
|
18 | if not found then |
| 379 | 8 | exp_dae := DAE.CREF(DAE.CREF_IDENT("("+intString(current_row)+")", DAE.T_REAL_DEFAULT, {}), DAE.T_REAL_DEFAULT); | |
| 380 | 8 | return; | |
| 381 | end if; | ||
| 382 | expressions := {}; | ||
| 383 | _:= match pre_mode | ||
| 384 | case 0 algorithm // mode for pivot-update operation | ||
| 385 | 6 | pivot_exp := traceEquation(row1, pre_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 386 | 6 | update_exp := traceEquation(row2, pre_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 387 | 6 | exp_dae := buildExpression(factor1, factor2, pivot_exp, update_exp); | |
| 388 | 6 | return; | |
| 389 | then pre_mode; | ||
| 390 | case 1 algorithm // mode for swap-rows operation | ||
| 391 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | if row1 == current_row then |
| 392 | 1 | exp_dae := traceEquation(row2, pre_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 393 | elseif row2 == current_row then | ||
| 394 | 1 | exp_dae := traceEquation(row1, pre_op, op_modes, op_val1, op_val2, op_val3, op_val4); | |
| 395 | end if; | ||
| 396 | 2 | return; | |
| 397 | then pre_mode; | ||
| 398 | case 2 algorithm // mode for gcd operation | ||
| 399 | //gcd := op_val2[i]; | ||
| 400 | 2 | exp_dae := DAE.BINARY(traceEquation(row1, pre_op, op_modes, op_val1, op_val2, op_val3, op_val4), DAE.DIV(DAE.T_REAL_DEFAULT), DAE.RCONST(gcd)); | |
| 401 | 2 | return; | |
| 402 | then pre_mode; | ||
| 403 | else pre_mode; | ||
| 404 | end match; | ||
| 405 | ✗ | return; | |
| 406 | end traceEquation; | ||
| 407 | |||
| 408 | function findLastOperation | ||
| 409 | "Finds the last elimination operation affecting the given row and returns the required information to reconstruct the previous state." | ||
| 410 | input Integer current_row, last_op; | ||
| 411 | input array<Integer> op_modes, op_val1, op_val2, op_val3, op_val4; | ||
| 412 | output Boolean found = false; | ||
| 413 | output Integer pre_op, pre_mode, row1, row2, factor1, factor2, gcd; | ||
| 414 | algorithm | ||
| 415 | (pre_op, pre_mode, row1, row2, factor1, factor2, gcd) := (0,0,0,0,0,0,0); | ||
| 416 | // trace back how each zero row was created | ||
| 417 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 14 times.
|
29 | for op in last_op:-1:1 loop // each operation backwards |
| 418 | // mode for pivot-update operation | ||
| 419 |
3/4✓ Branch 1 taken 6 times.
✓ Branch 2 taken 15 times.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
|
21 | if op_val3[op] == current_row and op_modes[op] == 0 then |
| 420 | pre_mode := 0; | ||
| 421 | 6 | row1 := op_val1[op]; | |
| 422 | row2 := op_val3[op]; | ||
| 423 | 6 | factor1 := op_val2[op]; | |
| 424 | 6 | factor2 := op_val4[op]; | |
| 425 | found := true; | ||
| 426 | 6 | pre_op := op - 1; | |
| 427 | 6 | return; | |
| 428 | // mode for swap-rows operation | ||
| 429 | elseif (op_val1[op] == current_row or op_val2[op] == current_row) and op_modes[op] == 1 then | ||
| 430 | pre_mode := 1; | ||
| 431 | row1 := op_val1[op]; | ||
| 432 | 2 | row2 := op_val2[op]; | |
| 433 | found := true; | ||
| 434 | 2 | pre_op := op - 1; | |
| 435 | 2 | return; | |
| 436 | // mode for gcd operation | ||
| 437 | elseif op_val1[op] == current_row and op_modes[op] == 2 then | ||
| 438 | pre_mode := 2; | ||
| 439 | row1 := op_val1[op]; | ||
| 440 | 2 | gcd := op_val2[op]; | |
| 441 | found := true; | ||
| 442 | 2 | pre_op := op - 1; | |
| 443 | 2 | return; | |
| 444 | end if; | ||
| 445 | end for; | ||
| 446 | end findLastOperation; | ||
| 447 | |||
| 448 | function buildExpression | ||
| 449 | "Constructs the symbolic expression resulting from a pivot-update operation." | ||
| 450 | input Integer factor_pivot, factor_update; | ||
| 451 | input DAE.Exp pivot_exp, update_exp; | ||
| 452 | output DAE.Exp exp_dae; | ||
| 453 | protected | ||
| 454 | DAE.Exp exp_dae_elem; | ||
| 455 | algorithm | ||
| 456 | // ToDo: switch to new simplify | ||
| 457 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
|
6 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 458 | ✗ | print("Step-by-step construction of the zero row:\n"); | |
| 459 | end if; | ||
| 460 | 6 | exp_dae := DAE.BINARY(DAE.BINARY(DAE.RCONST(factor_pivot), DAE.MUL(DAE.T_REAL_DEFAULT), update_exp), | |
| 461 | DAE.SUB(DAE.T_REAL_DEFAULT), | ||
| 462 | DAE.BINARY(DAE.RCONST(factor_update), DAE.MUL(DAE.T_REAL_DEFAULT), pivot_exp)); | ||
| 463 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
|
6 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 464 | ✗ | print(ExpressionBasics.printExpStr(exp_dae)+"\n"); | |
| 465 | end if; | ||
| 466 | 6 | (exp_dae, _) := ExpressionSimplify.simplify(exp_dae); | |
| 467 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 468 | ✗ | print("Simplified expression: "+ExpressionBasics.printExpStr(exp_dae)+"\n"); | |
| 469 | end if; | ||
| 470 | end buildExpression; | ||
| 471 | |||
| 472 | function createEquations | ||
| 473 | "Creates new equations according to the replacements and resolves cyclic alias equations." | ||
| 474 | input list<ComponentRef> vars; | ||
| 475 | input Pointer<Integer> index; | ||
| 476 | input array<list<Integer>> indices, values; | ||
| 477 | input Integer num_eqns; | ||
| 478 | input array<Expression> lhs_array; | ||
| 479 | output list<Pointer<Equation>> resolved_eqns = {}; | ||
| 480 | protected | ||
| 481 | Pointer<Equation> new_eq; | ||
| 482 | Expression cref_exp, sub_exp, rhs, lhs; | ||
| 483 | list<Integer> indices_list, values_list; | ||
| 484 | Status status; | ||
| 485 | Equation solved_eq; | ||
| 486 | algorithm | ||
| 487 | // iterate backwards due to upper-triangular dependency structure of the system | ||
| 488 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
19 | for i in num_eqns:-1:1 loop |
| 489 | 16 | rhs := Expression.makeInteger(0); | |
| 490 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
45 | for j in 1:listLength(indices[i]) loop |
| 491 | 29 | indices_list := indices[i]; | |
| 492 | 29 | cref_exp := Expression.fromCref(listGet(vars,listGet(indices_list, j)+1)); // indices are 0-based, Modelica lists are 1-based | |
| 493 | 29 | values_list := values[i]; | |
| 494 | 58 | sub_exp := Expression.MULTARY( | |
| 495 | arguments = {cref_exp, Expression.makeInteger(listGet(values_list, j))}, | ||
| 496 | inv_arguments = {}, | ||
| 497 | operator = Operator.makeMul(Type.REAL())); | ||
| 498 | 29 | sub_exp := SimplifyExp.simplify(sub_exp); | |
| 499 | 29 | rhs := Expression.MULTARY( | |
| 500 | arguments = {rhs, sub_exp}, | ||
| 501 | inv_arguments = {}, | ||
| 502 | operator = Operator.makeAdd(Expression.typeOf(sub_exp))); | ||
| 503 | end for; | ||
| 504 | 16 | rhs := SimplifyExp.simplify(rhs); | |
| 505 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
|
16 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 506 | ✗ | print("rhs: "+Expression.toString(rhs)+"\n"); | |
| 507 | end if; | ||
| 508 | 16 | new_eq := Equation.makeAssignment(rhs, lhs_array[i], index, NBEquation.TMP_STR, Iterator.EMPTY(), EquationAttributes.default(EquationKind.UNKNOWN, false)); | |
| 509 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
|
16 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 510 | ✗ | print("new_eq: "+Equation.toString(Pointer.access(new_eq))+"\n"); | |
| 511 | end if; | ||
| 512 | // solve equation for variable to eliminate cyclic dependencies | ||
| 513 | 16 | (solved_eq,status, _) := Solve.solveBody(Pointer.access(new_eq), listGet(vars,i), UnorderedMap.new<Function>(AbsynUtil.pathHash, AbsynUtil.pathEqual)); | |
| 514 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
|
16 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 515 | ✗ | print("solved_eq: "+Equation.toString(solved_eq)+"\n"); | |
| 516 | end if; | ||
| 517 | 16 | resolved_eqns := Pointer.create(solved_eq) :: resolved_eqns; | |
| 518 | end for; | ||
| 519 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | if Flags.isSet(Flags.DUMP_ASSC) then |
| 520 | ✗ | print("Number of equations: "+intString(listLength(resolved_eqns))+"\n"); | |
| 521 | ✗ | for eq_ptr in resolved_eqns loop | |
| 522 | ✗ | print("eq_ptr: "+Equation.toString(Pointer.access(eq_ptr))+"\n"); | |
| 523 | end for; | ||
| 524 | end if; | ||
| 525 | end createEquations; | ||
| 526 | |||
| 527 | 3 | function setMatrix | |
| 528 | input Integer nv "number of variables"; | ||
| 529 | input Integer ne "number of equations"; | ||
| 530 | input Integer nz "number of nonzero values"; | ||
| 531 | input array<list<Integer>> adj "adjacency matrix"; | ||
| 532 | input array<list<Integer>> val "value matrix"; | ||
| 533 | external "C" ASSC_setMatrix(nv,ne,nz,adj,val) annotation(Library = "omcruntime"); | ||
| 534 | end setMatrix; | ||
| 535 | |||
| 536 | 5 | function getMatrix | |
| 537 | input array<list<Integer>> adj "adjacency matrix"; | ||
| 538 | input array<list<Integer>> val "value matrix"; | ||
| 539 | external "C" ASSC_getMatrix(adj,val) annotation(Library = "omcruntime"); | ||
| 540 | end getMatrix; | ||
| 541 | |||
| 542 | 5 | function freeMatrix | |
| 543 | external "C" ASSC_freeMatrix() annotation(Library = "omcruntime"); | ||
| 544 | end freeMatrix; | ||
| 545 | |||
| 546 | 3 | function printMatrix | |
| 547 | external "C" ASSC_printMatrix() annotation(Library = "omcruntime"); | ||
| 548 | end printMatrix; | ||
| 549 | |||
| 550 | ✗ | function bareiss | |
| 551 | external "C" ASSC_bareiss() annotation(Library = "omcruntime"); | ||
| 552 | end bareiss; | ||
| 553 | |||
| 554 | 5 | function getNumberOfOperations | |
| 555 | input array<Integer> nop /* always size 1 */; | ||
| 556 | output Integer num; | ||
| 557 | external "C" num=ASSC_getNumberOfOperations(nop) annotation(Library = "omcruntime"); | ||
| 558 | end getNumberOfOperations; | ||
| 559 | |||
| 560 | 10 | function getOperations | |
| 561 | input array<Integer> op_modes; | ||
| 562 | input array<Integer> op_val1; | ||
| 563 | input array<Integer> op_val2; | ||
| 564 | input array<Integer> op_val3; | ||
| 565 | input array<Integer> op_val4; | ||
| 566 | external "C" ASSC_getOperations(op_modes, op_val1, op_val2, op_val3, op_val4) annotation(Library = "omcruntime"); | ||
| 567 | end getOperations; | ||
| 568 | |||
| 569 | protected | ||
| 570 | type CrefLst = list<ComponentRef>; | ||
| 571 | |||
| 572 | uniontype Tuple_Id | ||
| 573 | "tuple as key for UnorderedMap" | ||
| 574 | record TUPLE_ID | ||
| 575 | Pointer<Equation> eq_ptr; | ||
| 576 | ComponentRef cref; | ||
| 577 | end TUPLE_ID; | ||
| 578 | |||
| 579 | function toString | ||
| 580 | input Tuple_Id id; | ||
| 581 | output String str; | ||
| 582 | algorithm | ||
| 583 | 96 | str := Equation.toString(Pointer.access(id.eq_ptr)); | |
| 584 | 96 | str := BVariable.toString(BVariable.getVar(id.cref, sourceInfo())) + str; | |
| 585 | end toString; | ||
| 586 | |||
| 587 | function hash | ||
| 588 | "just hashes the id based on its string representation" | ||
| 589 | input Tuple_Id id; | ||
| 590 | output Integer hash; | ||
| 591 | algorithm | ||
| 592 | 96 | hash := stringHashDjb2(toString(id)); | |
| 593 | end hash; | ||
| 594 | |||
| 595 | function isEqual | ||
| 596 | input Tuple_Id id1; | ||
| 597 | input Tuple_Id id2; | ||
| 598 | output Boolean b; | ||
| 599 | algorithm | ||
| 600 |
2/4✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 48 times.
|
48 | b := Equation.isEqualPtr(id1.eq_ptr, id2.eq_ptr) and ComponentRef.isEqual(id1.cref, id2.cref); |
| 601 | end isEqual; | ||
| 602 | end Tuple_Id; | ||
| 603 | |||
| 604 | annotation(__OpenModelica_Interface="nbackend"); | ||
| 605 | end NBASSC; | ||
| 606 |