Linux GNU 11.4.0 Code Coverage Report


Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 88.9% 806 / 0 / 907
Functions: -% 0 / 1 / 1
Branches: 76.0% 427 / 0 / 562

OMCompiler/Compiler/BackEnd/ResolveLoops.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 ResolveLoops
37 " file: ResolveLoops.mo
38 package: ResolveLoops
39 description: This package contains functions for the optimization module
40 resolveLoops."
41
42
43 import BackendDAE;
44 import DAE;
45
46 protected
47 import Array;
48 import AvlSetInt;
49 import BackendDAEUtil;
50 import BackendEquation;
51 import BackendVariable;
52 import BackendVarTransform;
53 import BackendDump;
54 import ComponentReference;
55 protected import ComponentReferenceBasics;
56 import Expression;
57 import ExpressionBasics;
58 import ExpressionSimplify;
59 import ExpressionSolve;
60 import Flags;
61 import HpcOmTaskGraph;
62 import List;
63 import Tearing;
64 import Util;
65
66 protected type IntList = list<Integer>;
67
68 public function resolveLoops "author:Waurich TUD 2013-12
69 traverses the equations and finds simple equations(i.e. linear functions
70 withcoefficients of 1 or -1). if these equations form loops, they will be
71 contracted.
72 This happens especially in eletrical models. Here, kirchhoffs voltage and
73 current law can be applied."
74 input BackendDAE.BackendDAE inDAE;
75 output BackendDAE.BackendDAE outDAE;
76 protected
77 BackendDAE.EqSystems eqSysts;
78 BackendDAE.Shared shared;
79 algorithm
80 1071 (eqSysts, shared, _) := List.mapFold2(inDAE.eqs, resolveLoops_main, inDAE.shared, 1);
81 1071 outDAE := BackendDAE.DAE(eqSysts, shared);
82 end resolveLoops;
83
84 protected function resolveLoops_main "author: Waurich TUD 2014-01
85 Collects the linear equations of the whole DAE. bipartite graphs of the
86 eqSystem can be output. All variables and equations which do not belong to a
87 loop will be removed. the loops will be analysed and resolved"
88 input BackendDAE.EqSystem inEqSys;
89 input BackendDAE.Shared inShared "unused, just for dumping graphml";
90 input Integer inSysIdx;
91 output BackendDAE.EqSystem outEqSys;
92 output BackendDAE.Shared outShared = inShared "unused";
93 output Integer outSysIdx;
94 algorithm
95 (outEqSys, outSysIdx) := matchcontinue inEqSys
96 local
97 Integer numSimpEqs, numVars;
98 array<Integer> eqMapArr,varMapArr,nonLoopEqMark,markLinEqVars;
99 list<Integer> eqMapping;
100 list<list<Integer>> partitions;
101 list<tuple<Boolean,String>> varAtts,eqAtts;
102 BackendDAE.Variables vars,simpVars;
103 BackendDAE.EquationArray eqs,simpEqs;
104 BackendDAE.EqSystem syst;
105 BackendDAE.AdjacencyMatrix m,mT,m_cut, mT_cut, m_after;
106 list<BackendDAE.Equation> simpEqLst;
107 list<BackendDAE.Var> simpVarLst;
108
109 case syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqs)
110 algorithm
111 1572 (m,_) := BackendDAEUtil.adjacencyMatrix(syst, BackendDAE.ABSOLUTE(), NONE(), BackendDAEUtil.isInitializationDAE(inShared));
112
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1572 times.
1572 if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
113 ✗ BackendDump.dumpBipartiteGraphEqSystem(syst,inShared, "whole System_"+intString(inSysIdx));
114 end if;
115 //BackendDump.dumpEqSystem(syst,"the complete DAE");
116
117 // get the linear equations and their vars
118 1572 markLinEqVars := arrayCreate(BackendVariable.varsSize(vars),-1);
119 1572 (simpEqLst,eqMapping,_,_,markLinEqVars,_) := BackendEquation.traverseEquationArray(eqs, getSimpleEquations, ({},{},1,vars,markLinEqVars,m));
120 1572 eqMapArr := listArray(eqMapping);
121 1572 (simpVarLst,varMapArr) := getSimpleEquationVariables(markLinEqVars,vars);
122
123 1572 simpEqs := BackendEquation.listEquation(simpEqLst);
124 1572 simpVars := BackendVariable.listVar1(simpVarLst);
125
126 // build the adjacency matrix for the linear equations
127 1572 numSimpEqs := listLength(simpEqLst);
128 1572 numVars := listLength(simpVarLst);
129 1572 (m,mT) := BackendDAEUtil.adjacencyMatrixDispatch(simpVars,simpEqs, BackendDAE.ABSOLUTE(), NONE(), BackendDAEUtil.isInitializationDAE(inShared));
130
131
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1572 times.
1572 if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
132 ✗ varAtts := List.threadMap(List.fill(false,numVars),List.fill("",numVars),Util.makeTuple);
133 ✗ eqAtts := List.threadMap(List.fill(false,numSimpEqs),List.fill("",numSimpEqs),Util.makeTuple);
134 ✗ BackendDump.dumpBipartiteGraphStrongComponent2(simpVars,simpEqs,m,varAtts,eqAtts,"rL_simpEqs_"+intString(inSysIdx));
135 end if;
136
137 //partition graph
138 1572 partitions := partitionBipartiteGraph(m,mT);
139 1572 partitions := List.filterOnTrue(partitions,List.hasSeveralElements);
140 //print("the partitions for system "+intString(inSysIdx)+" : \n"+stringDelimitList(List.map(partitions,HpcOmTaskGraph.intLstString),"\n")+"\n");
141
142 // cut the deadends (vars and eqs outside of the loops)
143 1572 m_cut := arrayCopy(m);
144 1572 mT_cut := arrayCopy(mT);
145 1572 (_,nonLoopEqMark) := resolveLoops_cutNodes(m_cut,mT_cut); // this is pretty memory intensive
146
147
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1572 times.
1572 if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
148 ✗ varAtts := List.threadMap(List.fill(false,numVars),List.fill("",numVars),Util.makeTuple);
149 ✗ eqAtts := List.threadMap(List.fill(false,numSimpEqs),List.fill("",numSimpEqs),Util.makeTuple);
150 ✗ BackendDump.dumpBipartiteGraphStrongComponent2(simpVars,simpEqs,m_cut,varAtts,eqAtts,"rL_loops_"+intString(inSysIdx));
151 end if;
152
153 // handle the partitions separately, resolve the loops in the partitions, insert the resolved equation
154 1572 eqs := resolveLoops_resolvePartitions(partitions,m_cut,mT_cut,m,mT,eqMapArr,varMapArr,eqs,vars,nonLoopEqMark);
155 1572 syst.orderedEqs := eqs;
156 //BackendDump.dumpEquationList(eqLst,"the complete DAE after resolving");
157
158
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1572 times.
1572 if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
159 // get the graphML for the resolved System
160 ✗ simpEqLst := BackendEquation.getList(eqMapping,eqs);
161 ✗ simpEqs := BackendEquation.listEquation(simpEqLst);
162 ✗ numSimpEqs := listLength(simpEqLst);
163 ✗ numVars := listLength(simpVarLst);
164 ✗ m_after := BackendDAEUtil.adjacencyMatrixDispatch(simpVars,simpEqs, BackendDAE.ABSOLUTE(),NONE(),BackendDAEUtil.isInitializationDAE(inShared));
165 ✗ varAtts := List.threadMap(List.fill(false,numVars),List.fill("",numVars),Util.makeTuple);
166 ✗ eqAtts := List.threadMap(List.fill(false,numSimpEqs),List.fill("",numSimpEqs),Util.makeTuple);
167 ✗ BackendDump.dumpBipartiteGraphStrongComponent2(simpVars,simpEqs,m_after,varAtts,eqAtts,"rL_after_"+intString(inSysIdx));
168 end if;
169
170 1572 syst := BackendDAEUtil.clearEqSyst(syst);
171 //BackendDump.dumpEqSystem(eqSys,"the complete DAE after");
172 1572 then (syst, inSysIdx+1);
173
174 ✗ else (inEqSys, inSysIdx+1);
175 end matchcontinue;
176 end resolveLoops_main;
177
178 protected function resolveLoops_resolvePartitions "author:Waurich TUD 2014-02
179 checks every partition for loops and resolves them if its worth to."
180 input list<list<Integer>> partitionsIn;
181 input BackendDAE.AdjacencyMatrix mIn;
182 input BackendDAE.AdjacencyMatrixT mTIn;
183 input BackendDAE.AdjacencyMatrix m_uncut;
184 input BackendDAE.AdjacencyMatrixT mT_uncut;
185 input array<Integer> eqMap;
186 input array<Integer> varMap;
187 input BackendDAE.EquationArray daeEqs;
188 input BackendDAE.Variables daeVars;
189 input array<Integer> nonLoopEqMark;
190 output BackendDAE.EquationArray daeEqsOut;
191 algorithm
192 daeEqsOut := match partitionsIn
193 local
194 Option<tuple<list<Integer>,BackendDAE.AdjacencyMatrix,list<list<Integer>>>> optStructureMapping;
195 list<Integer> partition, eqCrossLst, varCrossLst, mapIndices;
196 list<list<Integer>> rest, loops;
197 BackendDAE.EquationArray eqs;
198 BackendDAE.AdjacencyMatrix map;
199 case partition::rest
200 algorithm
201 // search the partitions for loops
202 969 partition := List.filter1OnTrue(partition,arrayIsZeroAt,nonLoopEqMark); //the eqs that are loops
203
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 969 times.
969 if listEmpty(partition) then
204 ✗ eqs := resolveLoops_resolvePartitions(rest,mIn,mTIn,m_uncut,mT_uncut,eqMap,varMap,daeEqs,daeVars,nonLoopEqMark);
205 else
206 //print("\nanalyse the partition "+stringDelimitList(List.map(partition,intString),",")+"\n");
207 969 (loops,eqCrossLst,varCrossLst,optStructureMapping) := resolveLoops_findLoops({partition},mIn,mTIn);
208 //print("the loops in this partition: \n"+stringDelimitList(List.map(loops,HpcOmTaskGraph.intLstString),"\n")+"\n");
209
210 // check if its worth to resolve the loops
211
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 969 times.
✓ Branch 2 taken 942 times.
✓ Branch 3 taken 27 times.
969 if isSome(optStructureMapping) then
212
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
27 SOME((mapIndices,map,loops)) := optStructureMapping;
213 27 loops := List.filter1OnTrueAndUpdate(loops,evaluateTripleLoop,updateTripleLoop,tripleLoopInfo(m_uncut,arrayLength(mT_uncut),mapIndices,map));
214 else
215 942 loops := List.filterOnFalse(loops,listEmpty);
216 942 loops := List.filter1OnTrue(loops,evaluateLoop,(m_uncut,mT_uncut,eqCrossLst));
217 end if;
218 //print("the loops that will be resolved: \n"+stringDelimitList(List.map(loops,HpcOmTaskGraph.intLstString),"\n")+"\n");
219 // resolve the loops
220 969 (eqs,_) := resolveLoops_resolveAndReplace(loops,eqCrossLst,varCrossLst,mIn,mTIn,eqMap,varMap,daeEqs,daeVars,{}); //KAB4
221 969 eqs := resolveLoops_resolvePartitions(rest,mIn,mTIn,m_uncut,mT_uncut,eqMap,varMap,eqs,daeVars,nonLoopEqMark);
222 end if;
223 then
224 eqs;
225 case {}
226 algorithm
227 then
228 daeEqs;
229 end match;
230 end resolveLoops_resolvePartitions;
231
232 protected function resolveLoops_cutNodes "author: Waurich TUD 2014-01
233 cut the deadend nodes from the partitions"
234 input BackendDAE.AdjacencyMatrix mIn;
235 input BackendDAE.AdjacencyMatrix mTIn;
236 output array<Integer> deadEndVarsMark;
237 output array<Integer> deadEndEqsMark;
238 algorithm
239 (deadEndVarsMark,deadEndEqsMark) := matchcontinue mTIn
240 local
241 Integer numVars, numEqs, idx;
242 list<Integer> loopVars, loopEqs, nonLoopVars;
243 case _
244 algorithm
245 // get the outer deadEnd variables
246 numVars := arrayLength(mTIn);
247 numEqs := arrayLength(mIn);
248 1572 nonLoopVars := List.filter2OnTrue(List.intRange(numVars),arrayEntryLengthIs,mTIn,1);
249 //print("nonLoopVars: \n"+stringDelimitList(List.map(nonLoopVars,intString)," / ")+"\n");
250
251 1572 deadEndVarsMark := arrayCreate(numVars,0);
252 1572 deadEndEqsMark := arrayCreate(numVars,0);
253
254
2/2
✓ Branch 1 taken 11617 times.
✓ Branch 2 taken 1572 times.
13189 for idx in nonLoopVars loop
255 11617 arrayUpdate(deadEndVarsMark,idx,1);
256 end for;
257
258
2/2
✓ Branch 0 taken 11617 times.
✓ Branch 1 taken 1572 times.
13189 for idx in nonLoopVars loop
259 // DFS
260 11617 markDeadEndsInBipartiteGraph(idx,mIn,mTIn,deadEndEqsMark,deadEndVarsMark);
261 end for;
262
263 idx := 1;
264
2/2
✓ Branch 0 taken 19239 times.
✓ Branch 1 taken 1572 times.
20811 while idx <= numVars loop
265 // non-loop var
266
2/2
✓ Branch 1 taken 13222 times.
✓ Branch 2 taken 6017 times.
19239 if arrayGet(deadEndVarsMark,idx) == 1 then
267 13222 arrayUpdate(mTIn,idx,{});
268 //loop var
269 else
270 6017 loopEqs := arrayGet(mTIn,idx);
271 6017 loopEqs := List.filter1OnTrue(loopEqs,arrayIsZeroAt,deadEndEqsMark); //the loop equations
272 6017 arrayUpdate(mTIn,idx,loopEqs);
273 end if;
274 19239 idx := idx+1;
275 end while;
276
277 idx := 1;
278
2/2
✓ Branch 0 taken 9287 times.
✓ Branch 1 taken 1572 times.
10859 while idx <= numEqs loop
279 // non-loop eq
280
2/2
✓ Branch 1 taken 2131 times.
✓ Branch 2 taken 7156 times.
9287 if arrayGet(deadEndEqsMark,idx) == 1 then
281 2131 arrayUpdate(mIn,idx,{});
282 //loop eq
283 else
284 7156 loopVars := arrayGet(mIn,idx);
285 7156 loopVars := List.filter1OnTrue(loopVars,arrayIsZeroAt,deadEndVarsMark); //the loop vars
286 7156 arrayUpdate(mIn,idx,loopVars);
287 end if;
288 9287 idx := idx+1;
289 end while;
290
291 then (deadEndVarsMark,deadEndEqsMark);
292 else
293 algorithm
294 ✗ Error.addInternalError("function resolveLoops_cutNodes failed", sourceInfo());
295 ✗ then
296 fail();
297 end matchcontinue;
298 end resolveLoops_cutNodes;
299
300 protected function arrayEntryLengthIs "author:Waurich TUD 2014-01
301 gets the indexed entry of the array and compares the length with the given value."
302 input Integer idx;
303 input array<list<Integer>> arr;
304 input Integer len;
305 output Boolean eqLen;
306 protected
307 list<Integer> entry;
308 Integer len1;
309 algorithm
310 19239 entry := arrayGet(arr,idx);
311 19239 len1 := listLength(entry);
312 19239 eqLen := intEq(len,len1);
313 end arrayEntryLengthIs;
314
315 protected function getSimpleEquations
316 "if the linear equation contains only variables with factor 1 or -1 except for the states."
317 input BackendDAE.Equation inEq;
318 input tuple<list<BackendDAE.Equation>, list<Integer>, Integer, BackendDAE.Variables, array<Integer>, BackendDAE.AdjacencyMatrix> inTpl;
319 output BackendDAE.Equation outEq = inEq;
320 output tuple<list<BackendDAE.Equation>, list<Integer>, Integer, BackendDAE.Variables, array<Integer>, BackendDAE.AdjacencyMatrix> outTpl;
321 protected
322 Boolean isSimple;
323 Integer idx;
324 BackendDAE.Equation eq;
325 BackendDAE.AdjacencyMatrix m;
326 BackendDAE.Variables vars;
327 array<Integer> markLinEqVars;
328 list<BackendDAE.Equation> eqLst;
329 list<Integer> idxMap;
330 algorithm
331 49246 (eqLst, idxMap, idx, vars, markLinEqVars, m) := inTpl;
332
4/4
✓ Branch 1 taken 47690 times.
✓ Branch 2 taken 1556 times.
✓ Branch 4 taken 47626 times.
✓ Branch 5 taken 64 times.
49246 if BackendEquation.isEquation(inEq) and not eqIsConst(inEq)/*simple assignments should not occur here, they cannot be improved any further*/ then
333 47626 (eq,(isSimple,_)) := BackendEquation.traverseExpsOfEquation(inEq,isAddOrSubExp,(true,vars));
334
2/2
✓ Branch 0 taken 9287 times.
✓ Branch 1 taken 38339 times.
47626 if isSimple then
335 eqLst := eq::eqLst;
336 idxMap := idx::idxMap;
337
2/2
✓ Branch 1 taken 28500 times.
✓ Branch 2 taken 9287 times.
37787 for varIdx in m[idx] loop
338 28500 arrayUpdate(markLinEqVars,intAbs(varIdx),1);
339 end for;
340 end if;
341 end if;
342 49246 outTpl := (eqLst,idxMap,idx+1,vars,markLinEqVars,m);
343 end getSimpleEquations;
344
345 protected function getSimpleEquationVariables"
346 Get the variables which occur in the linear equations and have been marked in the markLinEqVars array.
347 author:Waurich 2017-07"
348 input array<Integer> markLinEqVars;
349 input BackendDAE.Variables vars;
350 output list<BackendDAE.Var> simpVars={};
351 output array<Integer> varMapArr;
352 protected
353 Integer varIdx;
354 list<Integer> varMap;
355 algorithm
356 varMap := {};
357
2/2
✓ Branch 0 taken 1485 times.
✓ Branch 1 taken 87 times.
57384 for varIdx in 1:arrayLength(markLinEqVars) loop
358
2/2
✓ Branch 1 taken 19239 times.
✓ Branch 2 taken 36573 times.
55812 if markLinEqVars[varIdx] > 0 then
359 varMap := varIdx::varMap;
360 19239 simpVars := BackendVariable.getVarAt(vars, varIdx)::simpVars;
361 end if;
362 end for;
363 1572 varMapArr := listArray(varMap);
364 end getSimpleEquationVariables;
365
366 public function resolveLoops_findLoops "author:Waurich TUD 2014-02
367 gets the crossNodes for the partitions and searches for loops"
368 input list<list<Integer>> partitionsIn;
369 input BackendDAE.AdjacencyMatrix mIn; // the whole system of simpleEquations
370 input BackendDAE.AdjacencyMatrixT mTIn;
371 input Boolean findExactlyOneLoop=false;
372 output list<list<Integer>> loopsOut = {};
373 output list<Integer> crossEqsOut = {};
374 output list<Integer> crossVarsOut = {};
375 output Option<tuple<list<Integer>,BackendDAE.AdjacencyMatrix,list<list<Integer>>>> optStructureMapping = NONE();
376 protected
377 list<list<Integer>> loops, eqVars;
378 list<Integer> eqCrossLst, varCrossLst, partitionVars;
379 AvlSetInt.Tree set;
380 algorithm
381
2/2
✓ Branch 0 taken 5167 times.
✓ Branch 1 taken 5167 times.
10334 for partition in partitionsIn loop
382 try
383 // get the eqCrossNodes and varCrossNodes i.e. nodes with more than 2 edges
384 5167 eqVars := List.map1(partition,Array.getIndexFirst,mIn);
385 set := AvlSetInt.EMPTY();
386
2/2
✓ Branch 0 taken 21064 times.
✓ Branch 1 taken 5167 times.
26231 for vars in eqVars loop
387 21064 set := AvlSetInt.addList(set, vars);
388 end for;
389 5167 partitionVars := AvlSetInt.listKeys(set);
390 5167 eqCrossLst := List.fold2(partition,gatherCrossNodes,mIn,mTIn,{});
391 5167 varCrossLst := List.fold2(partitionVars,gatherCrossNodes,mTIn,mIn,{});
392
393 // search the partitions for loops
394 5167 (loops,optStructureMapping) := resolveLoops_findLoops2(partition,eqCrossLst,varCrossLst,mIn,mTIn,findExactlyOneLoop);
395
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5167 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
5167 if if findExactlyOneLoop then (not listEmpty(loops) and not listEmpty(loopsOut)) else false then
396 ✗ fail();
397 end if;
398 5167 loopsOut := listAppend(loops,loopsOut);
399
2/2
✓ Branch 0 taken 969 times.
✓ Branch 1 taken 4198 times.
5167 if isPresent(crossEqsOut) then
400 969 crossEqsOut := listAppend(eqCrossLst,crossEqsOut);
401 end if;
402
2/2
✓ Branch 0 taken 969 times.
✓ Branch 1 taken 4198 times.
5167 if isPresent(crossVarsOut) then
403 969 crossVarsOut := listAppend(varCrossLst,crossVarsOut);
404 end if;
405 else
406 ✗ return;
407 end try;
408 end for;
409 end resolveLoops_findLoops;
410
411 protected function resolveLoops_findLoops2 "author: Waurich TUD 2014-01
412 handles the given partition of eqs and vars depending whether there are only varCrossNodes, only EqCrossNodes, both of them or none of them."
413 input list<Integer> eqsIn;
414 input list<Integer> eqCrossLstIn;
415 input list<Integer> varCrossLstIn;
416 input BackendDAE.AdjacencyMatrix mIn; // the whole system of simpleEquations
417 input BackendDAE.AdjacencyMatrixT mTIn;
418 input Boolean findExactlyOneLoop;
419 output list<list<Integer>> loopsOut;
420 output Option<tuple<list<Integer>,BackendDAE.AdjacencyMatrix,list<list<Integer>>>> structureMapping;
421 algorithm
422 (loopsOut,structureMapping) := match(eqCrossLstIn, varCrossLstIn)
423 local
424 Boolean isNoSingleLoop;
425 list<Integer> eqCrossLst, subLoop, mapIndices;
426 list<list<Integer>> paths, allPaths, simpleLoops, tripleLoops, paths0, paths1, closedPaths, loopConnectors, connectedPaths;
427 AvlSetInt.Tree eqCrossSet;
428 BackendDAE.AdjacencyMatrix minAdj,map;
429 tuple<list<Integer>,BackendDAE.AdjacencyMatrix> mapping;
430 Option<tuple<list<Integer>,BackendDAE.AdjacencyMatrix,list<list<Integer>>>> optTripleMapping;
431 case(_::_, {})
432 algorithm //KAB1
433 //print("partition has only eqCrossNodes\n");
434 // get the paths between the crossEqNodes and order them according to their length
435 2210 allPaths := getPathTillNextCrossEq(eqCrossLstIn,mIn,mTIn,eqCrossLstIn,{},{});
436 2210 allPaths := List.sort(allPaths,List.listIsLonger);
437 //print("all paths: \n"+stringDelimitList(List.map(allPaths,HpcOmTaskGraph.intLstString)," / ")+"\n");
438 2210 paths1 := List.fold1(allPaths,getReverseDoubles,allPaths,{}); // all paths with just one direction
439 //UNUSED paths0 := List.unique(paths1); // only the paths between the eqs without concerning the vars in between
440
441 2210 simpleLoops := getDoubles(paths1,{}); // get 2 adjacent equations which form a simple loop i.e. they share 2 variables
442 //print("all simpleLoop-paths: \n"+stringDelimitList(List.map(simpleLoops,HpcOmTaskGraph.intLstString)," / ")+"\n");
443 2210 (_,paths,_) := List.intersection1OnTrue(paths1,simpleLoops,intLstIsEqual);
444
445 // special case to find more complex structures (arrays and triple loops)
446
2/2
✓ Branch 0 taken 1678 times.
✓ Branch 1 taken 532 times.
2210 if listEmpty(simpleLoops) then
447 1678 (eqCrossLst,paths1,mapping,minAdj) := findEqualPathStructure(eqCrossLstIn,paths1);
448 //print("crossNodes: " + HpcOmTaskGraph.intLstString(eqCrossLst) + "\n");
449 1678 (mapIndices,map) := mapping;
450 1678 (tripleLoops,paths0) := getTriples(eqCrossLst,minAdj);
451 1678 optTripleMapping := SOME((mapIndices,map,tripleLoops));
452 //print("all tripleLoop-paths after equal structure: \n"+stringDelimitList(List.map(tripleLoops,HpcOmTaskGraph.intLstString)," / ")+"\n");
453 else
454 optTripleMapping := NONE();
455 532 paths0 := List.sort(paths,List.listIsLonger); // solve the small loops first
456 532 (connectedPaths,loopConnectors) := connect2PathsToLoops(paths0,{},{});
457 532 loopConnectors := List.filter1OnTrue(loopConnectors,connectsLoops,simpleLoops);
458 532 simpleLoops := listAppend(simpleLoops,loopConnectors) annotation(__OpenModelica_DisableListAppendWarning=true);
459
460 //print("all simpleLoop-paths: \n"+stringDelimitList(List.map(simpleLoops,HpcOmTaskGraph.intLstString)," / ")+"\n");
461 532 subLoop := connectPathsToOneLoop(simpleLoops,{}); // try to build a a closed loop from these paths
462 532 isNoSingleLoop := listEmpty(subLoop);
463
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 521 times.
532 simpleLoops := if isNoSingleLoop then simpleLoops else {subLoop};
464 532 paths0 := listAppend(simpleLoops,connectedPaths);
465 532 paths0 := sortPathsAsChain(paths0);
466
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 532 times.
532 if findExactlyOneLoop then
467 ✗ if not listEmpty(paths0) then
468 ✗ {_} := paths0;
469 end if;
470 end if;
471 end if;
472
473 //print("all paths to be resolved: \n"+stringDelimitList(List.map(paths0,HpcOmTaskGraph.intLstString)," / ")+"\n");
474 2210 then
475 (paths0, optTripleMapping);
476 case({}, _::_)
477 algorithm
478 //print("partition has only varCrossNodes\n");
479 // get the paths between the crossVarNodes and order them according to their length
480 30 paths := getPathTillNextCrossEq(varCrossLstIn,mTIn,mIn,varCrossLstIn,{},{});
481 30 paths := List.sort(paths,List.listIsLonger);
482 30 paths := listReverse(paths);
483 //print("from all the paths: \n"+stringDelimitList(List.map(paths,HpcOmTaskGraph.intLstString)," / ")+"\n");
484
485 30 (paths0,paths1) := List.extract1OnTrue(paths,listLengthIs,listLength(List.last(paths)));
486 //print("the shortest paths: \n"+stringDelimitList(List.map(paths0,HpcOmTaskGraph.intLstString)," / ")+"\n");
487
488
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 4 times.
30 paths1 := if listEmpty(paths1) then paths0 else paths1;
489 30 closedPaths := List.map1(paths1,closePathDirectly,paths0);
490 30 closedPaths := List.fold1(closedPaths,getReverseDoubles,closedPaths,{}); // all paths with just one direction
491 30 closedPaths := List.map(closedPaths,List.unique);
492 30 closedPaths := List.map1(closedPaths,getEqNodesForVarLoop,mTIn);// get the eqs for these varLoops
493 //print("solve the smallest loops: \n"+stringDelimitList(List.map(closedPaths,HpcOmTaskGraph.intLstString)," / ")+"\n");
494
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 30 times.
30 if findExactlyOneLoop then
495 ✗ if not listEmpty(closedPaths) then
496 ✗ {_} := closedPaths;
497 end if;
498 end if;
499 then (closedPaths,NONE());
500 case({}, {})
501 algorithm
502 // no crossNodes
503 //print("no crossNodes\n");
504 subLoop := eqsIn;
505
2/2
✓ Branch 0 taken 969 times.
✓ Branch 1 taken 1647 times.
2616 for e in eqsIn loop
506
2/2
✓ Branch 1 taken 239 times.
✓ Branch 2 taken 730 times.
969 if listEmpty(mIn[e]) then
507 subLoop := {};
508 break;
509 end if;
510 end for;
511 then ({subLoop},NONE());
512 case(_::_, _::_)
513 algorithm
514 //print("there are both varCrossNodes and eqNodes\n");
515 //at least get paths of length 2 between eqCrossNodes
516
1/2
✓ Branch 0 taken 550 times.
✗ Branch 1 not taken.
62456 for i in 1:arrayLength(mIn) loop
517 61906 arrayUpdate(mIn, i, List.heapSortIntList(mIn[i]));
518 end for;
519
1/2
✓ Branch 0 taken 550 times.
✗ Branch 1 not taken.
104492 for i in 1:arrayLength(mTIn) loop
520 103942 arrayUpdate(mTIn, i, List.heapSortIntList(mTIn[i]));
521 end for;
522 550 eqCrossSet := AvlSetInt.addList(AvlSetInt.EMPTY(), eqCrossLstIn);
523 550 paths := getShortPathsBetweenEqCrossNodes(AvlSetInt.listKeysReverse(eqCrossSet), eqCrossSet, mIn, mTIn, {}, findExactlyOneLoop);
524 //
525 //print("GOT SOME NEW LOOPS: \n"+stringDelimitList(List.map(paths,HpcOmTaskGraph.intLstString)," / ")+"\n");
526 then (paths,NONE());
527 else
528 algorithm
529 ✗ Error.addInternalError("function resolveLoops_findLoops2 failed", sourceInfo());
530 ✗ then
531 fail();
532 end match;
533 end resolveLoops_findLoops2;
534
535 protected function findEqualPathStructure //KAB2
536 "author: kabdelhak FHB 2019-06
537 Resolves seemingly hard structures mostly formed through array connectors.
538 By finding equally connected crossNodes and merging them to one single
539 node, the loop becomes way simpler. All merged nodes will get treated the
540 same as the super node. The mapping contains the information which nodes
541 got merged and which index contains the super node."
542 input output list<Integer> crossNodes;
543 input output list<list<Integer>> uniquePaths;
544 output tuple<list<Integer>,BackendDAE.AdjacencyMatrix> mapping;
545 output BackendDAE.AdjacencyMatrix minAdj;
546 protected
547 list<Integer> mapIndices;
548 BackendDAE.AdjacencyMatrix map;
549 algorithm
550 1678 minAdj := getMinimalAdjacencyMatrix(crossNodes,uniquePaths);
551
4/4
✓ Branch 0 taken 3041 times.
✓ Branch 1 taken 1678 times.
✓ Branch 2 taken 3041 times.
✓ Branch 3 taken 1678 times.
4719 (minAdj,uniquePaths,mapIndices,map,crossNodes) := removeEqualPaths(crossNodes,minAdj,uniquePaths,{},arrayCreate(max(cn for cn in crossNodes),{}));
552 1678 mapping := (mapIndices,map);
553 end findEqualPathStructure;
554
555 protected function getMinimalAdjacencyMatrix
556 "author: kabdelhak FHB 2019-06
557 Returns the minimal adjacency matrix of connected equations in the partition.
558 NOTE: Rows AND Columns represent equations. No variables involved."
559 input list<Integer> crossNodes;
560 input list<list<Integer>> uniquePaths;
561 output BackendDAE.AdjacencyMatrix minAdj;
562 algorithm
563
4/4
✓ Branch 0 taken 3041 times.
✓ Branch 1 taken 1678 times.
✓ Branch 2 taken 3041 times.
✓ Branch 3 taken 1678 times.
4719 minAdj := arrayCreate(max(cn for cn in crossNodes),{});
564
2/2
✓ Branch 0 taken 1510 times.
✓ Branch 1 taken 1678 times.
3188 for path in uniquePaths loop
565 _ := match path
566 local
567 Integer a, b;
568 case a::{b} algorithm // a::rest ?
569 1451 minAdj := Array.consToElement(a, b, minAdj);
570 1451 minAdj := Array.consToElement(b, a, minAdj);
571 then 0;
572 else 1;
573 end match;
574 end for;
575
576 //sort
577
2/2
✓ Branch 0 taken 3041 times.
✓ Branch 1 taken 1678 times.
4719 for cn in crossNodes loop
578 3041 arrayUpdate(minAdj,cn,List.sort(arrayGet(minAdj,cn),intGt));
579 end for;
580 end getMinimalAdjacencyMatrix;
581
582 protected function removeEqualPaths
583 "author: kabdelhak FHB 2019-06
584 Helper function for findEqualPathStructure. It removes the superfluous nodes
585 and declares the index first as super node."
586 input list<Integer> crossNodes;
587 input output BackendDAE.AdjacencyMatrix minAdj;
588 input output list<list<Integer>> uniquePaths;
589 input output list<Integer> mapIndices;
590 input output BackendDAE.AdjacencyMatrix map;
591 output list<Integer> accCrossNodes = {};
592 protected
593 UnorderedMap<IntList, Integer> groups;
594 array<Integer> groupOf = arrayCreate(arrayLength(minAdj), 0);
595 array<Boolean> merged = arrayCreate(arrayLength(minAdj), false), collected = arrayCreate(arrayLength(minAdj), false);
596 Integer cn1, numGroups = 0, numMerged = 0;
597 list<Integer> row, nodes = crossNodes, rest, assigned, unassigned;
598 algorithm
599 1678 groups := UnorderedMap.new<Integer>(hashIntList, HpcOmTaskGraph.equalLists);
600
2/2
✓ Branch 0 taken 3041 times.
✓ Branch 1 taken 1678 times.
4719 for node in crossNodes loop
601 3041 row := arrayGet(minAdj, node);
602
2/2
✓ Branch 1 taken 2885 times.
✓ Branch 2 taken 156 times.
3041 if not UnorderedMap.contains(row, groups) then
603 2885 numGroups := numGroups + 1;
604 2885 UnorderedMap.add(row, numGroups, groups);
605 end if;
606 3041 arrayUpdate(groupOf, node, UnorderedMap.getOrFail(row, groups));
607 end for;
608
2/2
✓ Branch 0 taken 2885 times.
✓ Branch 1 taken 1678 times.
4563 while not listEmpty(nodes) loop
609 2885 cn1 :: rest := nodes;
610
2/2
✓ Branch 1 taken 1678 times.
✓ Branch 2 taken 1207 times.
2885 if not collected[cn1] then
611 1678 arrayUpdate(collected, cn1, true);
612 accCrossNodes := cn1 :: accCrossNodes;
613 end if;
614 assigned := {};
615 unassigned := {};
616
2/2
✓ Branch 0 taken 4337 times.
✓ Branch 1 taken 2885 times.
7222 for cn2 in rest loop
617
2/2
✓ Branch 2 taken 156 times.
✓ Branch 3 taken 4181 times.
4337 if groupOf[cn2] == groupOf[cn1] then
618 assigned := cn2 :: assigned;
619 156 arrayUpdate(minAdj, cn2, {});
620 156 arrayUpdate(merged, cn2, true);
621 156 numMerged := numMerged + 1;
622 else
623 unassigned := cn2 :: unassigned;
624
2/2
✓ Branch 1 taken 1274 times.
✓ Branch 2 taken 2907 times.
4181 if not collected[cn2] then
625 1274 arrayUpdate(collected, cn2, true);
626 accCrossNodes := cn2 :: accCrossNodes;
627 end if;
628 end if;
629 end for;
630
2/2
✓ Branch 0 taken 2745 times.
✓ Branch 1 taken 140 times.
2885 if not listEmpty(assigned) then
631 mapIndices := cn1 :: mapIndices;
632 140 map := Array.appendToElement(cn1, assigned, map);
633 end if;
634 nodes := unassigned;
635 end while;
636
7/8
✓ Branch 1 taken 175 times.
✓ Branch 2 taken 1335 times.
✓ Branch 3 taken 1510 times.
✓ Branch 4 taken 1678 times.
✓ Branch 5 taken 1335 times.
✓ Branch 6 taken 1678 times.
✓ Branch 7 taken 1678 times.
✗ Branch 8 not taken.
3188 uniquePaths := list(path for path guard not pathContainsMerged(path, merged) in uniquePaths);
637 // removing them one at a time reversed the list once per node
638
2/2
✓ Branch 0 taken 1554 times.
✓ Branch 1 taken 124 times.
1678 if intMod(numMerged, 2) == 1 then
639 124 uniquePaths := listReverse(uniquePaths);
640 end if;
641 end removeEqualPaths;
642
643 protected function hashIntList
644 input list<Integer> lst;
645 output Integer hash = 17;
646 algorithm
647
2/2
✓ Branch 0 taken 8540 times.
✓ Branch 1 taken 8967 times.
17507 for i in lst loop
648
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8540 times.
8540 hash := intMod(hash * 31 + i, 65599);
649 end for;
650 end hashIntList;
651
652 protected function pathContainsMerged
653 input list<Integer> path;
654 input array<Boolean> merged;
655 output Boolean c = false;
656 algorithm
657
2/2
✓ Branch 0 taken 2927 times.
✓ Branch 1 taken 1335 times.
4262 for n in path loop
658
5/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2927 times.
✓ Branch 2 taken 2885 times.
✓ Branch 3 taken 42 times.
✓ Branch 5 taken 175 times.
✓ Branch 6 taken 2710 times.
5854 if n <= arrayLength(merged) and merged[n] then
659 c := true;
660 175 return;
661 end if;
662 end for;
663 end pathContainsMerged;
664
665 protected function listContains
666 input list<Integer> lst;
667 input Integer int;
668 output Boolean res = false;
669 algorithm
670
2/2
✓ Branch 0 taken 626 times.
✓ Branch 1 taken 312 times.
938 for i in lst loop
671
2/2
✓ Branch 0 taken 73 times.
✓ Branch 1 taken 553 times.
626 if intEq(i,int) then
672 res := true;
673 73 return;
674 end if;
675 end for;
676 end listContains;
677
678 protected function hasSameIntSortedExcept
679 input list<Integer> inList1;
680 input list<Integer> inList2;
681 input Integer excl;
682 output Boolean rv = false;
683 protected
684 Integer i1, i2;
685 list<Integer> l1 = inList1, l2 = inList2;
686 algorithm
687
2/4
✓ Branch 0 taken 48378 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48378 times.
48378 if listEmpty(inList1) or listEmpty(inList2) then
688 ✗ return;
689 end if;
690 48378 i1::l1 := l1;
691 48378 i2::l2 := l2;
692 while true loop
693
2/2
✓ Branch 0 taken 64754 times.
✓ Branch 1 taken 108021 times.
172775 if i1 > i2 then
694
2/2
✓ Branch 0 taken 2932 times.
✓ Branch 1 taken 61822 times.
64754 if listEmpty(l2) then
695 2932 return;
696 end if;
697 61822 i2::l2 := l2;
698 elseif i1 < i2 then
699
2/2
✓ Branch 0 taken 2450 times.
✓ Branch 1 taken 43855 times.
46305 if listEmpty(l1) then
700 2450 return;
701 end if;
702 43855 i1::l1 := l1;
703 else
704
2/2
✓ Branch 0 taken 37973 times.
✓ Branch 1 taken 23743 times.
61716 if i1 <> excl then
705 rv := true;
706 37973 return;
707 end if;
708
4/4
✓ Branch 0 taken 20711 times.
✓ Branch 1 taken 3032 times.
✓ Branch 2 taken 1991 times.
✓ Branch 3 taken 18720 times.
23743 if listEmpty(l1) or listEmpty(l2) then
709 5023 return;
710 end if;
711 18720 i1::l1 := l1;
712 18720 i2::l2 := l2;
713 end if;
714 end while;
715 end hasSameIntSortedExcept;
716
717 protected function getShortPathsBetweenEqCrossNodes"find closedLoops between 2 eqCrossNode, no matter whether there are var cross nodes between them.
718 author: vwaurich TUD 12-2016"
719 input list<Integer> eqCrossLstIn;
720 input AvlSetInt.Tree eqCrossSet;
721 input BackendDAE.AdjacencyMatrix mIn "rows sorted ascending";
722 input BackendDAE.AdjacencyMatrixT mTIn "rows sorted ascending";
723 input list<list<Integer>> pathsIn;
724 input Boolean findExactlyOneLoop;
725 output list<list<Integer>> pathsOut = pathsIn;
726 protected
727 Integer hub;
728 list<Integer> adjVars, adjEqs, newPath;
729 list<list<Integer>> paths;
730 algorithm
731
2/2
✓ Branch 0 taken 14540 times.
✓ Branch 1 taken 550 times.
15090 for crossEq in eqCrossLstIn loop
732 paths := {};
733 14540 adjVars := arrayGet(mIn, crossEq);
734 // the row of a variable shared by most equations is not walked
735 14540 hub := longestRow(adjVars, mTIn);
736
2/2
✓ Branch 0 taken 48905 times.
✓ Branch 1 taken 14540 times.
63445 for adjVar in adjVars loop
737
2/2
✓ Branch 0 taken 14540 times.
✓ Branch 1 taken 34365 times.
48905 adjEqs := if adjVar == hub then eqsSharingVia(adjVars, hub, mIn, mTIn) else arrayGet(mTIn, adjVar);
738 //all adjEqs which are crossnodes as well
739
2/2
✓ Branch 0 taken 146940 times.
✓ Branch 1 taken 48905 times.
195845 for adjEq in adjEqs loop
740
4/4
✓ Branch 0 taken 48630 times.
✓ Branch 1 taken 98310 times.
✓ Branch 3 taken 252 times.
✓ Branch 4 taken 48378 times.
146940 if if adjEq > crossEq then (not AvlSetInt.hasKey(eqCrossSet, adjEq)) else true then
741 98562 continue;
742 end if;
743
2/2
✓ Branch 2 taken 37973 times.
✓ Branch 3 taken 10405 times.
48378 if hasSameIntSortedExcept(adjVars, arrayGet(mIn, adjEq), adjVar) then
744 newPath := adjEq::{crossEq};
745 37973 paths := List.unionElt(newPath, paths);
746
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 37973 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
37973 if if findExactlyOneLoop then (not listEmpty(pathsOut)) else false then
747 ✗ fail();
748 end if;
749 end if;
750 end for;
751 end for;
752 14540 pathsOut := listAppend(paths, pathsOut);
753 end for;
754 end getShortPathsBetweenEqCrossNodes;
755
756 protected function longestRow "walks the rows in lockstep, so a long row costs no more than the others"
757 input list<Integer> vars;
758 input BackendDAE.AdjacencyMatrixT mT;
759 output Integer var = 0;
760 protected
761 list<tuple<Integer, list<Integer>>> rows = list((v, arrayGet(mT, v)) for v in vars), left;
762 Integer v;
763 list<Integer> row;
764 algorithm
765
2/2
✓ Branch 0 taken 74878 times.
✓ Branch 1 taken 6772 times.
81650 while not listEmpty(rows) loop
766 74878 (var, _) := listHead(rows);
767
2/2
✓ Branch 1 taken 7768 times.
✓ Branch 2 taken 67110 times.
74878 if listEmpty(listRest(rows)) then
768 7768 return;
769 end if;
770 left := {};
771
2/2
✓ Branch 0 taken 210104 times.
✓ Branch 1 taken 67110 times.
277214 for r in rows loop
772 210104 (v, row) := r;
773
2/2
✓ Branch 0 taken 168967 times.
✓ Branch 1 taken 41137 times.
210104 if not listEmpty(row) then
774 168967 left := (v, listRest(row)) :: left;
775 end if;
776 end for;
777 67110 rows := listReverse(left);
778 end while;
779 end longestRow;
780
781 protected function eqsSharingVia "the equations containing var and another of vars, ascending like the rows of mT"
782 input list<Integer> vars;
783 input Integer var;
784 input BackendDAE.AdjacencyMatrix m;
785 input BackendDAE.AdjacencyMatrixT mT;
786 output list<Integer> eqs = {};
787 algorithm
788
2/2
✓ Branch 0 taken 48905 times.
✓ Branch 1 taken 14540 times.
63445 for v in vars loop
789
2/2
✓ Branch 0 taken 34365 times.
✓ Branch 1 taken 14540 times.
48905 if v <> var then
790
2/2
✓ Branch 1 taken 108629 times.
✓ Branch 2 taken 34365 times.
142994 for eq in arrayGet(mT, v) loop
791
2/2
✓ Branch 2 taken 80180 times.
✓ Branch 3 taken 28449 times.
108629 if sortedListContains(arrayGet(m, eq), var) then
792 eqs := eq :: eqs;
793 end if;
794 end for;
795 end if;
796 end for;
797 14540 eqs := List.sortedUnique(List.sort(eqs, intGt), intEq);
798 end eqsSharingVia;
799
800 protected function sortedListContains
801 input list<Integer> lst "ascending";
802 input Integer x;
803 output Boolean found = false;
804 algorithm
805
2/2
✓ Branch 0 taken 249343 times.
✓ Branch 1 taken 4775 times.
254118 for i in lst loop
806
2/2
✓ Branch 0 taken 103854 times.
✓ Branch 1 taken 145489 times.
249343 if i >= x then
807 103854 found := i == x;
808 103854 return;
809 end if;
810 end for;
811 end sortedListContains;
812
813 protected function connectsLoops "author:Waurich TUD 2014-02
814 checks if the given path connects 2 closed simple Loops"
815 input list<Integer> path;
816 input list<list<Integer>> allLoops;
817 output Boolean connected;
818 protected
819 Boolean b1, b2;
820 Integer startNode, endNode;
821 list<list<Integer>> loops1, loops2;
822 algorithm
823 151 startNode := listHead(path);
824 151 endNode := List.last(path);
825 // the startNode is connected to a loop
826 151 loops1 := List.filter1OnTrue(allLoops,firstInListIsEqual,startNode);
827 151 loops2 := List.filter1OnTrue(allLoops,lastInListIsEqual,startNode);
828
4/4
✓ Branch 0 taken 99 times.
✓ Branch 1 taken 52 times.
✓ Branch 2 taken 57 times.
✓ Branch 3 taken 42 times.
151 b1 := (not listEmpty(loops1)) or (not listEmpty(loops2));
829 // the endNode is connected to a loop
830 151 loops1 := List.filter1OnTrue(allLoops,firstInListIsEqual,endNode);
831 151 loops2 := List.filter1OnTrue(allLoops,lastInListIsEqual,endNode);
832
4/4
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 31 times.
✓ Branch 2 taken 65 times.
✓ Branch 3 taken 55 times.
151 b2 := (not listEmpty(loops1)) or (not listEmpty(loops2));
833 151 connected := b1 and b2;
834 end connectsLoops;
835
836 protected function connectPathsToOneLoop "author:Waurich TUD 2014-02
837 tries to connect various paths to one closed, simple loop"
838 input list<list<Integer>> allPathsIn;
839 input list<Integer> loopIn;
840 output list<Integer> loopOut;
841 algorithm
842 loopOut := matchcontinue(allPathsIn,loopIn)
843 local
844 Integer startNode, endNode;
845 list<Integer> path, nextPath, restPath;
846 list<list<Integer>> rest, nextPaths1, nextPaths2;
847 case(_,startNode::path)
848 algorithm
849 123 endNode := List.last(path);
850
2/2
✓ Branch 0 taken 112 times.
✓ Branch 1 taken 11 times.
123 true := intEq(startNode,endNode);
851 then
852 path;
853 case(_,startNode::_)
854 algorithm
855 // TODO: This makes a list of all matching paths when it seems to really
856 // only need the first matching. Same in the case below.
857 112 nextPaths1 := List.filter1OnTrue(allPathsIn, firstInListIsEqual, startNode);
858 112 nextPaths2 := List.filter1OnTrue(allPathsIn, lastInListIsEqual, startNode);
859 112 nextPaths2 := listAppend(nextPaths1,nextPaths2);
860 112 nextPath := listHead(nextPaths2);
861 53 rest := List.deleteMemberOnTrue(nextPath,allPathsIn, function List.isEqualOnTrue(inCompFunc = intEq));
862 53 nextPath := List.deleteMemberOnTrue(startNode,nextPath,intEq);
863 53 path := listAppend(nextPath,loopIn);
864 53 path := connectPathsToOneLoop(rest,path);
865 then
866 path;
867 case(path::rest,{})
868 algorithm
869
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 532 times.
532 startNode::restPath := path;
870 532 nextPaths1 := List.filter1OnTrue(rest, firstInListIsEqual, startNode);
871 532 nextPaths2 := List.filter1OnTrue(rest, lastInListIsEqual, startNode);
872 532 nextPaths2 := listAppend(nextPaths1,nextPaths2);
873 532 nextPath := listHead(nextPaths2);
874 70 rest := List.deleteMemberOnTrue(nextPath,rest, function List.isEqualOnTrue(inCompFunc = intEq));
875 70 path := listAppend(nextPath,restPath);
876 70 path := connectPathsToOneLoop(rest,path);
877 then
878 path;
879 else
880 algorithm
881 then
882 {};
883 end matchcontinue;
884 end connectPathsToOneLoop;
885
886 protected function resolveLoops_resolveAndReplace "author:Waurich TUD 2014-01
887 resolves a singleLoop. depending on whether there are only eqCrossNodes, varCrossNodes, both or none."
888 input list<list<Integer>> loopsIn;
889 input list<Integer> eqCrossLstIn;
890 input list<Integer> varCrossLstIn;
891 input BackendDAE.AdjacencyMatrix mIn;
892 input BackendDAE.AdjacencyMatrixT mTIn;
893 input array<Integer> eqMap;
894 input array<Integer> varMap;
895 input BackendDAE.EquationArray daeEqsIn;
896 input BackendDAE.Variables daeVarsIn;
897 input list<Integer> replEqsIn;
898 output BackendDAE.EquationArray daeEqsOut;
899 output list<Integer> replEqsOut;
900 algorithm
901 (daeEqsOut,replEqsOut) := match(loopsIn, eqCrossLstIn, varCrossLstIn)
902 local
903 Integer pos,eq1,eq2;
904 list<Integer> loop1, eqs, vars, crossEqs, crossVars, replEqs, loopVars, adjVars, m_row;
905 list<list<Integer>> rest, eqVars;
906 BackendDAE.Equation resolvedEq;
907 BackendDAE.EquationArray daeEqs;
908 case({}, _, _)
909 algorithm
910 then
911 (daeEqsIn,replEqsIn);
912 case(loop1::rest, _::crossEqs, {})
913 algorithm
914 // only eqCrossNodes
915 //print("only eqCrossNodes\n");
916 169 loop1 := List.unique(loop1);
917 169 (resolvedEq, m_row) := resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
918
919 // get the equation that will be replaced and the rest
920 169 (crossEqs,eqs,_) := List.intersection1OnTrue(loop1,eqCrossLstIn,intEq); // replace a crossEq in the loop
921 169 replEqs := List.intersectionOnTrue(replEqsIn,loop1,intEq); // just consider the already replaced equations in this loop
922
923 // first try to replace a non cross node, otherwise an already replaced eq, or if none of them is available take a crossnode (THIS IS NOT YET CLEAR)
924
2/2
✓ Branch 0 taken 29 times.
✓ Branch 1 taken 140 times.
169 if not listEmpty(eqs) then
925 29 pos := listHead(eqs);
926 elseif not listEmpty(replEqs) then
927 59 pos := listHead(replEqs);
928 elseif not listEmpty(crossEqs) then
929 81 pos := listHead(crossEqs);
930 else
931 pos := -1;
932 end if;
933
934 169 eqs := List.deleteMemberOnTrue(pos,loop1,intEq);
935 //print("contract eqs: "+stringDelimitList(List.map(eqs,intString),",")+" to eq "+intString(pos)+"\n");
936
937 // get the corresponding vars
938 169 eqVars := List.map1(loop1,Array.getIndexFirst,mIn);
939 169 vars := List.flatten(eqVars);
940 169 loopVars := doubleEntriesInLst(vars); // the vars in the loop
941 169 (_,adjVars,_) := List.intersection1OnTrue(vars,loopVars,intEq); // the vars adjacent to the loop
942
943 // update adjacencyMatrix
944 169 List.map2_0(loopVars,Array.updateIndexFirst,{},mTIn); //delete the vars in the loop
945 169 List.map2_0(adjVars,arrayGetDeleteInLst,loop1,mTIn); // remove the loop eqs from the adjacent vars
946 169 List.map2_0(adjVars,arrayGetAppendLst,{pos},mTIn); // redirect the adjacent vars to the replaced eq
947 169 List.map2_0(loop1,Array.updateIndexFirst,{},mIn); //delete the eqs in the loop
948 169 arrayUpdate(mIn,pos,adjVars); // redirect the replaced equation to the vars outside of the loops
949
950 // update remaining paths
951 169 rest := List.map2(rest,replaceContractedNodes,pos,eqs);
952 169 rest := List.unique(rest);
953 //print("the remaining paths: "+stringDelimitList(List.map(rest,HpcOmTaskGraph.intLstString),"\n")+"\n\n");
954
955 // replace Equation
956 //print("replace equation "+intString(pos)+"\n");
957 replEqs := pos::replEqsIn;
958 169 arrayUpdate(mIn,pos,m_row);
959 169 pos := arrayGet(eqMap,pos);
960 169 daeEqs := BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);
961
962 495 (daeEqs,replEqs) := resolveLoops_resolveAndReplace(rest,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqs,daeVarsIn,replEqs);
963 then
964 (daeEqs,replEqs);
965 case(loop1::rest, {}, _::crossVars)
966 algorithm
967 // only varCrossNodes
968 //print("only varCrossNodes\n");
969 70 loop1 := List.unique(loop1);
970 70 (resolvedEq, m_row) := resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
971
972 // get the equation that will be replaced and the rest
973 70 (replEqs,_,eqs) := List.intersection1OnTrue(replEqsIn,loop1,intEq); // just consider the already replaced equations in this loop
974
975 //priorize the not yet replaced equations
976 70 eqs := priorizeEqsWithVarCrosses(eqs,mIn,varCrossLstIn);
977 //print("priorized eqs: "+stringDelimitList(List.map(eqs,intString),",")+"\n");
978
979 // first try to replace a non cross node, otherwise an already replaced eq
980
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 47 times.
70 pos := if not listEmpty(replEqs) then listHead(replEqs) else -1;
981
1/2
✓ Branch 0 taken 70 times.
✗ Branch 1 not taken.
70 pos := if not listEmpty(eqs) then listHead(eqs) else pos;
982
983 70 eqs := List.deleteMemberOnTrue(pos,loop1,intEq);
984 //print("contract eqs: "+stringDelimitList(List.map(eqs,intString),",")+" to eq "+intString(pos)+"\n");
985
986 // get the corresponding vars
987 70 eqVars := List.map1(loop1,Array.getIndexFirst,mIn);
988 70 vars := List.flatten(eqVars);
989 70 loopVars := doubleEntriesInLst(vars); // the vars in the loop
990 70 (crossVars,loopVars,_) := List.intersection1OnTrue(loopVars,varCrossLstIn,intEq); // some crossVars have to remain
991 //print("loopVars: "+stringDelimitList(List.map(loopVars,intString),",")+"\n");
992
993 70 (_,adjVars,_) := List.intersection1OnTrue(vars,loopVars,intEq); // the vars adjacent to the loop
994 70 adjVars := listAppend(crossVars,adjVars);
995 70 adjVars := List.unique(adjVars);
996
997 // update adjacencyMatrix
998 70 List.map2_0(loopVars,Array.updateIndexFirst,{},mTIn); //delete the vars in the loop
999 70 List.map2_0(adjVars,arrayGetDeleteInLst,loop1,mTIn); // remove the loop eqs from the adjacent vars
1000 70 List.map2_0(adjVars,arrayGetAppendLst,{pos},mTIn); // redirect the adjacent vars to the replaced eq
1001 70 List.map2_0(loop1,Array.updateIndexFirst,{},mIn); //delete the eqs in the loop
1002 70 arrayUpdate(mIn,pos,adjVars); // redirect the replaced equation to the vars outside of the loops
1003
1004 // update remaining paths
1005 70 rest := List.map2(rest,replaceContractedNodes,pos,eqs);
1006 70 rest := List.unique(rest);
1007 //print("the remaining paths: "+stringDelimitList(List.map(rest,HpcOmTaskGraph.intLstString),"\n")+"\n\n");
1008
1009 // replace Equation
1010 //print("replace equation "+intString(pos)+"\n");
1011 replEqs := pos::replEqsIn;
1012 70 arrayUpdate(mIn,pos,m_row);
1013 70 pos := arrayGet(eqMap,pos);
1014 70 daeEqs := BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);
1015
1016 70 (daeEqs,replEqs) := resolveLoops_resolveAndReplace(rest,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqs,daeVarsIn,replEqs);
1017 then
1018 (daeEqs,replEqs);
1019 case(loop1::rest, {}, {})
1020 algorithm
1021 // single Loop
1022 62 loop1 := List.unique(loop1);
1023 //print("single loop\n");
1024 62 (resolvedEq, m_row) := resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
1025
1026 // update AdjacencyMatrix
1027 62 (_,crossEqs,_) := List.intersection1OnTrue(loop1,replEqsIn,intEq); // do not replace an already replaced Eq
1028
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 62 times.
62 pos::_ := crossEqs; // the equation that will be replaced = pos
1029 62 eqVars := List.map1(loop1,Array.getIndexFirst,mIn);
1030 62 vars := List.flatten(eqVars);
1031 //print("delete vars: "+stringDelimitList(List.map(vars,intString),",")+" in the eqs: "+stringDelimitList(List.map(crossEqs,intString),",")+"\n");
1032 62 List.map2_0(loop1,Array.updateIndexFirst,{},mIn); //delete the equations in the loop
1033 62 List.map2_0(vars,Array.updateIndexFirst,{},mTIn); //delete the vars from the loop
1034
1035 // replace Equation
1036 //print("replace equation "+intString(pos)+"\n");
1037 replEqs := pos::replEqsIn;
1038 62 arrayUpdate(mIn,pos,m_row);
1039 62 pos := arrayGet(eqMap,pos);
1040 62 daeEqs := BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);
1041
1042 62 (daeEqs,replEqs) := resolveLoops_resolveAndReplace(rest,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqs,daeVarsIn,replEqs);
1043 then
1044 (daeEqs,replEqs);
1045 case(loop1::rest, _::_, _::_)
1046 algorithm
1047 // both eqCrossNodes and varCrossNodes, at least try the small loops
1048 //print("both eqCrossNodes and varCrossNodes, loopLength"+intString(listLength(loop1))+"\n");
1049 194 loop1 := List.unique(loop1);
1050
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 194 times.
194 true := listLength(loop1) == 2;
1051 194 (resolvedEq, m_row) := resolveClosedLoop(loop1,mIn,mTIn,eqMap,varMap,daeEqsIn,daeVarsIn);
1052 //print("resolved eq to "+BackendDump.equationString(resolvedEq)+"\n");
1053
1054
2/2
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 180 times.
194 if eqIsConst(resolvedEq)then
1055 //replace the equations that had the most vars
1056
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 14 times.
14 {eq1, eq2} := loop1;
1057 //print("eq1 "+intString(eq1)+" has vars: "+stringDelimitList(List.map(arrayGet(mIn,eq1),intString),",")+"\n");
1058 //print("eq2 "+intString(eq2)+" has vars: "+stringDelimitList(List.map(arrayGet(mIn,eq2),intString),",")+"\n");
1059
1060 //TODO: check if the assigned cref occurs in this equation!!!!
1061
1062
2/2
✓ Branch 8 taken 12 times.
✓ Branch 9 taken 2 times.
14 if listLength(BackendEquation.equationVars(BackendEquation.get(daeEqsIn,arrayGet(eqMap,eq1)),daeVarsIn)) >= listLength(BackendEquation.equationVars(BackendEquation.get(daeEqsIn,arrayGet(eqMap,eq2)),daeVarsIn)) then
1063 pos := eq1;
1064 else
1065 pos := eq2;
1066 end if;
1067 replEqs := pos::replEqsIn;
1068 //print("contract eqs: "+stringDelimitList(List.map(loop1,intString),",")+" to eq "+intString(pos)+"\n");
1069 14 arrayUpdate(mIn,pos,m_row);
1070 14 pos := arrayGet(eqMap,pos);
1071 14 daeEqs := BackendEquation.setAtIndex(daeEqsIn,pos,resolvedEq);
1072 else
1073 replEqs := replEqsIn;
1074 daeEqs := daeEqsIn;
1075 end if;
1076 194 (daeEqs,replEqs) := resolveLoops_resolveAndReplace(rest,eqCrossLstIn,varCrossLstIn,mIn,mTIn,eqMap,varMap,daeEqs,daeVarsIn,replEqs);
1077 then
1078 (daeEqs,replEqs);
1079 end match;
1080 end resolveLoops_resolveAndReplace;
1081
1082 protected function eqIsConst"outputs true if the equation is a constant assignment.
1083 author: vwaurich TUD 2017-01"
1084 input BackendDAE.Equation eq;
1085 output Boolean b;
1086 algorithm
1087 b :=match eq
1088 case BackendDAE.EQUATION(exp=DAE.RCONST(),scalar=DAE.CREF())
1089 then true;
1090 case BackendDAE.EQUATION(exp=DAE.CREF(),scalar=DAE.RCONST())
1091 then true;
1092 else then false;
1093 end match;
1094 end eqIsConst;
1095
1096 protected function arrayIsZeroAt
1097 input Integer pos;
1098 input array<Integer> arr;
1099 output Boolean isZero;
1100 algorithm
1101 96493 isZero := intEq(0,arr[pos]);
1102 end arrayIsZeroAt;
1103
1104 protected function markDeadEndsInBipartiteGraph "author:Waurich TUD 2017-01
1105 creates marking arrays for the equations and variables.
1106 if the node is a non-loop not its marked with 1, otherwise with0
1107 update the adjacencymatrix.
1108 "
1109 input Integer varIdx; //deadEnd
1110 input BackendDAE.AdjacencyMatrix mIn; // the rows correspond to the primary nodes
1111 input BackendDAE.AdjacencyMatrixT mTIn;
1112 input array<Integer> deadEndEqs; //marks all equations which are not in a loop with 1, otherwise 0
1113 input array<Integer> deadEndVars; //marks all variables which are not in a loop with 1, otherwise 0
1114 protected
1115 Integer eqIdx, var = varIdx;
1116 list<Integer> adjEqs, adjVars;
1117 Boolean walk = true;
1118 algorithm
1119 while walk loop
1120 walk := false;
1121 13748 adjEqs := arrayGet(mTIn,var);
1122 13748 adjEqs := List.filter1OnTrue(adjEqs,arrayIsZeroAt,deadEndEqs); //the eqs that are not yet marked
1123
2/2
✓ Branch 1 taken 12207 times.
✓ Branch 2 taken 1541 times.
13748 if listLength(adjEqs) == 1 then // the var is a dead end var
1124 12207 eqIdx := listHead(adjEqs);
1125 12207 arrayUpdate(deadEndVars,var,1);
1126 12207 adjVars := arrayGet(mIn,eqIdx);
1127 12207 adjVars := List.filter1OnTrue(adjVars,arrayIsZeroAt,deadEndVars); //the vars that are not yet marked
1128
2/2
✓ Branch 1 taken 2131 times.
✓ Branch 2 taken 10076 times.
12207 if listLength(adjVars) == 1 then //the adjacent equation is a dead end eq
1129 2131 arrayUpdate(deadEndEqs,eqIdx,1);
1130 2131 var := listHead(adjVars);
1131 walk := true;
1132 end if;
1133 end if;
1134 end while;
1135 end markDeadEndsInBipartiteGraph;
1136
1137 protected function arrayGetDeleteInLst "deletes all entries given in delEntries from the indexed list<Integer> of the array"
1138 input Integer idx;
1139 input list<Integer> delEntries;
1140 input array<list<Integer>> arrIn;
1141 protected
1142 list<Integer> entry;
1143 algorithm
1144 330 entry := arrayGet(arrIn,idx);
1145 330 (_,entry,_) := List.intersection1OnTrue(entry,delEntries,intEq);
1146 330 arrayUpdate(arrIn,idx,entry);
1147 end arrayGetDeleteInLst;
1148
1149 protected function arrayGetAppendLst "appends appLst to the indexed list<Integer> of the array"
1150 input Integer idx;
1151 input list<Integer> appLst;
1152 input array<list<Integer>> arrIn;
1153 protected
1154 list<Integer> entry;
1155 algorithm
1156 570 entry := arrayGet(arrIn,idx);
1157 570 arrayUpdate(arrIn,idx,listAppend(entry,appLst));
1158 end arrayGetAppendLst;
1159
1160 protected function getReverseDoubles "author: Waurich TUD 2014-01
1161 fold function to get the reversed doubles in a list."
1162 input list<Integer> elem;
1163 input list<list<Integer>> elemLst;
1164 input list<list<Integer>> foldLstIn;
1165 output list<list<Integer>> foldLstOut;
1166 replaceable type ElementType subtypeof Any;
1167 algorithm
1168 foldLstOut := matchcontinue foldLstIn
1169 local
1170 list<Integer> elemR;
1171 list<list<Integer>> foldLst;
1172 case _
1173 algorithm
1174 5904 elemR := listReverse(elem);
1175 5904 elemR := List.getMember(elemR,elemLst);
1176 5898 foldLst := List.deleteMemberOnTrue(elem,foldLstIn,function List.isEqualOnTrue(inCompFunc = intEq));
1177 then
1178 elemR::foldLst;
1179 else
1180 algorithm
1181 then
1182 foldLstIn;
1183 end matchcontinue;
1184 end getReverseDoubles;
1185
1186 protected function getDoubles "author: Waurich TUD 2014-01
1187 function to get the reversed doubles in a list."
1188 input list<list<Integer>> elemLstIn;
1189 input list<list<Integer>> lstIn;
1190 output list<list<Integer>> lstOut = lstIn;
1191 replaceable type ElementType subtypeof Any;
1192 protected
1193 list<Integer> elem;
1194 list<list<Integer>> elemLst = elemLstIn;
1195 algorithm
1196
2/2
✓ Branch 0 taken 2862 times.
✓ Branch 1 taken 2210 times.
5072 while not listEmpty(elemLst) loop
1197 2862 elem::elemLst := elemLst;
1198
2/2
✓ Branch 1 taken 2227 times.
✓ Branch 2 taken 635 times.
2862 if listMember(elem,elemLst) then
1199 lstOut := elem::lstOut;
1200 end if;
1201 end while;
1202 end getDoubles;
1203
1204 protected function getTriples
1205 "author: kabdelhak FHB 2019-07
1206 function to get loops containing three eqCrossNodes."
1207 input list<Integer> crossNodes;
1208 input BackendDAE.AdjacencyMatrix minAdj;
1209 output list<list<Integer>> tripleLoops = {};
1210 output list<list<Integer>> allPaths = {};
1211 protected
1212 list<Integer> path1, path2, path3;
1213 algorithm
1214
2/2
✓ Branch 0 taken 2952 times.
✓ Branch 1 taken 1678 times.
4630 for c0 in crossNodes loop
1215 2952 path1 := arrayGet(minAdj,c0);
1216
2/2
✓ Branch 0 taken 2736 times.
✓ Branch 1 taken 2952 times.
5688 for c1 in path1 loop
1217
2/2
✓ Branch 0 taken 1411 times.
✓ Branch 1 taken 1325 times.
2736 if intGt(c1,c0) then
1218 1411 path2 := arrayGet(minAdj,c1);
1219
2/2
✓ Branch 0 taken 2286 times.
✓ Branch 1 taken 1411 times.
3697 for c2 in path2 loop
1220
2/2
✓ Branch 0 taken 385 times.
✓ Branch 1 taken 1901 times.
2286 if intGt(c2,c1) then
1221 385 path3 := arrayGet(minAdj,c2);
1222
2/2
✓ Branch 1 taken 73 times.
✓ Branch 2 taken 312 times.
385 if listContains(path3,c0) then
1223 tripleLoops := {c0,c1,c2}::tripleLoops;
1224 allPaths := {c1,c2}::allPaths;
1225 allPaths := {c0,c2}::allPaths;
1226 allPaths := {c0,c1}::allPaths;
1227 end if;
1228 end if;
1229 end for;
1230 end if;
1231 end for;
1232 end for;
1233 end getTriples;
1234
1235 protected function getEqNodesForVarLoop "author: Waurich TUD 2013-01
1236 fold function to get the eqs in a loop that is given by the varNodes."
1237 input list<Integer> varIdcs;
1238 input BackendDAE.AdjacencyMatrixT mTIn;
1239 output list<Integer> eqIdcs;
1240 protected
1241 list<list<Integer>> varEqLst;
1242 list<Integer> eqLst;
1243 algorithm
1244 87 varEqLst := List.map1(varIdcs,Array.getIndexFirst,mTIn); // get the eqNodes from these loops
1245 87 eqLst := List.flatten(varEqLst);
1246 87 eqIdcs := doubleEntriesInLst(eqLst);
1247 end getEqNodesForVarLoop;
1248
1249 protected function resolveClosedLoop "author:Waurich TUD 2014-02
1250 sums up all equations in a loop so that the variables shared by the equations disappear."
1251 input list<Integer> loopIn;
1252 input BackendDAE.AdjacencyMatrix m;
1253 input BackendDAE.AdjacencyMatrixT mT;
1254 input array<Integer> eqMap;
1255 input array<Integer> varMap;
1256 input BackendDAE.EquationArray daeEqsIn;
1257 input BackendDAE.Variables daeVarsIn;
1258 output BackendDAE.Equation eqOut;
1259 output list<Integer> m_row;
1260 protected
1261 Integer startEqIdx,startEqDaeIdx;
1262 list<Integer> loop1, restLoop;
1263 BackendDAE.Equation eq;
1264 algorithm
1265
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 495 times.
495 startEqIdx::restLoop := loopIn;
1266 495 startEqDaeIdx := arrayGet(eqMap,startEqIdx);
1267 495 loop1 := sortLoop(restLoop,m,mT,{startEqIdx});
1268
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 495 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
495 if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) and listLength(loop1) > 1 then
1269 ✗ print("solve the loop: " + List.toString(loop1, intString) + "\n");
1270 end if;
1271 495 eq := BackendEquation.get(daeEqsIn,startEqDaeIdx);
1272 495 (eqOut, m_row) := resolveClosedLoop2(eq,loop1,m, arrayGet(m,startEqIdx), eqMap,varMap,daeEqsIn,daeVarsIn);
1273 end resolveClosedLoop;
1274
1275 protected function resolveClosedLoop2 "author:Waurich TUD 2013-12"
1276 input output BackendDAE.Equation eq;
1277 input list<Integer> loopIn;
1278 input BackendDAE.AdjacencyMatrix m;
1279 input output list<Integer> m_row "row of eq in m";
1280 input array<Integer> eqMap;
1281 input array<Integer> varMap;
1282 input BackendDAE.EquationArray daeEqsIn;
1283 input BackendDAE.Variables daeVarsIn;
1284 algorithm
1285 (eq, m_row) := match loopIn
1286 local
1287 Boolean algSign;
1288 Integer eqIdx2;
1289 list<Integer> adjVars, adjVars1, adjVars2, restLoop, posVars, negVars, nonUnitVars;
1290 list<DAE.ComponentRef> adjCrefs;
1291 BackendDAE.Equation eq2, eq3, resolvedEq;
1292 BackendVarTransform.VariableReplacements replacements;
1293 case {_} then (eq, m_row);
1294 case _::eqIdx2::restLoop algorithm
1295 // the equation to add
1296 639 eq2 := BackendEquation.get(daeEqsIn, arrayGet(eqMap,eqIdx2));
1297
1298 // get the vars that are shared of the 2 equations
1299 639 adjVars1 := m_row;
1300 639 adjVars2 := arrayGet(m,eqIdx2);
1301 639 (adjVars, adjVars1, adjVars2) := List.intersection1OnTrue(adjVars1, adjVars2, intEq);
1302
1303 // Only shared variables with a +/-1 coefficient in BOTH equations can be cancelled by
1304 // adding/subtracting the equations (the cancellation below replaces them with zero).
1305 // A variable scaled by a non-unit factor (e.g. a state, or a zero-sequence current
1306 // n*i0 = sum(i)) must NOT be cancelled this way, otherwise its coefficient would be
1307 // silently dropped and the resolved equation would be wrong (#13292). Keep such vars.
1308 639 (adjVars, nonUnitVars) := List.splitOnTrue(adjVars, function varIsUnitCoeff(varMap = varMap, daeVarsIn = daeVarsIn, eq1 = eq, eq2 = eq2));
1309
1310 // split shared vars by sign
1311 639 (posVars, negVars) := List.splitOnTrue(adjVars, function varSign(varMap = varMap, daeVarsIn = daeVarsIn, eq1 = eq, eq2 = eq2));
1312 639 algSign := listLength(posVars) > listLength(negVars); // choose set with more canceling vars
1313
6/6
✓ Branch 0 taken 229 times.
✓ Branch 1 taken 410 times.
✓ Branch 2 taken 1214 times.
✓ Branch 3 taken 639 times.
✓ Branch 4 taken 1214 times.
✓ Branch 5 taken 639 times.
1853 adjCrefs := list(crefFromIndex(idx, varMap, daeVarsIn) for idx in (if algSign then posVars else negVars));
1314
1315 // cancelled crefs are removed from adjacency; non-cancellable shared vars are kept
1316
2/2
✓ Branch 0 taken 410 times.
✓ Branch 1 taken 229 times.
1278 m_row := List.flatten({adjVars1, adjVars2, nonUnitVars, (if algSign then negVars else posVars)});
1317
1318 // replace `cref` with zero to make the job easier for `simplify`
1319 639 replacements := BackendVarTransform.emptyReplacementsSized(listLength(adjCrefs));
1320
4/4
✓ Branch 0 taken 1214 times.
✓ Branch 1 taken 639 times.
✓ Branch 2 taken 1214 times.
✓ Branch 3 taken 639 times.
1853 replacements := BackendVarTransform.addReplacements(replacements, adjCrefs, list(Expression.createZeroExpression(ComponentReference.crefTypeFull(c)) for c in adjCrefs), NONE());
1321
3/6
✗ Branch 1 not taken.
✓ Branch 2 taken 639 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 639 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 639 times.
639 ({resolvedEq, eq3}, _) := BackendVarTransform.replaceEquations({eq, eq2}, replacements, NONE());
1322
1323 639 resolvedEq := sumUp2Equations(algSign,resolvedEq,eq3);
1324
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 639 times.
639 if Flags.isSet(Flags.RESOLVE_LOOPS_DUMP) then
1325 ✗ print("From eqs \n"+BackendDump.equationString(eq)+"\n"+BackendDump.equationString(eq2)+"\n");
1326 ✗ print("resolved the eq \n"+BackendDump.equationString(resolvedEq)+"\n\n");
1327 end if;
1328 639 then resolveClosedLoop2(resolvedEq,eqIdx2::restLoop,m, m_row, eqMap,varMap,daeEqsIn,daeVarsIn);
1329 end match;
1330 end resolveClosedLoop2;
1331
1332 protected function crefFromIndex
1333 input Integer varIdx;
1334 input array<Integer> varMap;
1335 input BackendDAE.Variables daeVarsIn;
1336 output DAE.ComponentRef cref;
1337 protected
1338 Integer daeVarIdx;
1339 BackendDAE.Var var;
1340 algorithm
1341 3719 daeVarIdx := arrayGet(varMap, varIdx);
1342 3719 var := BackendVariable.getVarAt(daeVarsIn, daeVarIdx);
1343 3719 cref := BackendVariable.varCref(var);
1344 end crefFromIndex;
1345
1346 protected function varSign
1347 input Integer index;
1348 input array<Integer> varMap;
1349 input BackendDAE.Variables daeVarsIn;
1350 input BackendDAE.Equation eq1;
1351 input BackendDAE.Equation eq2;
1352 output Boolean algSign;
1353 protected
1354 DAE.ComponentRef cref = crefFromIndex(index, varMap, daeVarsIn);
1355 algorithm
1356 1244 algSign := CRefIsPosOnRHS(cref, eq1) <> CRefIsPosOnRHS(cref, eq2) "check the algebraic signs"; // XOR
1357 end varSign;
1358
1359 protected function varIsUnitCoeff "author: #13292
1360 true if the variable (given by its loop-local index) occurs with a +/-1 coefficient in
1361 both equations and can therefore be cancelled by adding/subtracting the two equations."
1362 input Integer index;
1363 input array<Integer> varMap;
1364 input BackendDAE.Variables daeVarsIn;
1365 input BackendDAE.Equation eq1;
1366 input BackendDAE.Equation eq2;
1367 output Boolean isUnit;
1368 protected
1369 DAE.ComponentRef cref = crefFromIndex(index, varMap, daeVarsIn);
1370 algorithm
1371
4/4
✓ Branch 1 taken 1257 times.
✓ Branch 2 taken 4 times.
✓ Branch 4 taken 13 times.
✓ Branch 5 taken 1244 times.
1261 isUnit := crefHasUnitCoeff(cref, eq1) and crefHasUnitCoeff(cref, eq2);
1372 end varIsUnitCoeff;
1373
1374 protected function crefHasUnitCoeff "author: #13292
1375 true unless the cref appears scaled by a constant other than +/-1 in the equation."
1376 input DAE.ComponentRef cref;
1377 input BackendDAE.Equation eq;
1378 output Boolean isUnit;
1379 algorithm
1380 isUnit := match eq
1381 local
1382 DAE.Exp e1, e2;
1383 case BackendDAE.EQUATION(exp = e1, scalar = e2)
1384
3/4
✓ Branch 1 taken 2518 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17 times.
✓ Branch 5 taken 2501 times.
2518 then crefUnitCoeffInExp(e1, cref) and crefUnitCoeffInExp(e2, cref);
1385 else true;
1386 end match;
1387 end crefHasUnitCoeff;
1388
1389 protected function crefUnitCoeffInExp "author: #13292
1390 false if cref appears multiplied by a non-(+/-1) constant in exp; true otherwise
1391 (cref absent, or present with a +/-1 coefficient)."
1392 input DAE.Exp exp;
1393 input DAE.ComponentRef cref;
1394 output Boolean isUnit;
1395 algorithm
1396 isUnit := match exp
1397 local
1398 DAE.Exp e1, e2;
1399 DAE.ComponentRef c;
1400 case DAE.BINARY(exp1 = e1, operator = DAE.ADD(), exp2 = e2)
1401
4/4
✓ Branch 1 taken 2940 times.
✓ Branch 2 taken 10 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 2924 times.
2950 then crefUnitCoeffInExp(e1, cref) and crefUnitCoeffInExp(e2, cref);
1402 case DAE.BINARY(exp1 = e1, operator = DAE.SUB(), exp2 = e2)
1403
3/4
✓ Branch 1 taken 2797 times.
✓ Branch 2 taken 1 time.
✗ Branch 4 not taken.
✓ Branch 5 taken 2797 times.
2798 then crefUnitCoeffInExp(e1, cref) and crefUnitCoeffInExp(e2, cref);
1404 case DAE.UNARY(exp = e1)
1405 944 then crefUnitCoeffInExp(e1, cref);
1406 case DAE.BINARY(exp1 = DAE.CREF(componentRef = c), operator = DAE.MUL(), exp2 = e2)
1407 ✗ then not ComponentReferenceBasics.crefEqualNoStringCompare(cref, c) or Expression.isOne(e2) or Expression.isConstMinusOne(e2);
1408 case DAE.BINARY(exp1 = e1, operator = DAE.MUL(), exp2 = DAE.CREF(componentRef = c))
1409
4/6
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 105 times.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 17 times.
122 then not ComponentReferenceBasics.crefEqualNoStringCompare(cref, c) or Expression.isOne(e1) or Expression.isConstMinusOne(e1);
1410 else true;
1411 end match;
1412 end crefUnitCoeffInExp;
1413
1414 public function sortLoop "author:Waurich TUD 2014-01
1415 sorts the equations in a loop so that they are solved in a row."
1416 input list<Integer> loopIn;
1417 input BackendDAE.AdjacencyMatrix m;
1418 input BackendDAE.AdjacencyMatrixT mT;
1419 input list<Integer> sortLoopIn;
1420 output list<Integer> sortLoopOut;
1421 algorithm
1422 sortLoopOut := match(loopIn, sortLoopIn)
1423 local
1424 Integer start, next;
1425 list<Integer> rest, vars, eqs;
1426 list<list<Integer>> varEqs;
1427 case({}, _)
1428 algorithm
1429 495 then
1430 listReverse(sortLoopIn);
1431 case(_, start::_)
1432 algorithm
1433 639 vars := arrayGet(m,start);
1434 639 varEqs := List.map1(vars,Array.getIndexFirst,mT);
1435 639 eqs := List.flatten(varEqs);
1436 639 eqs := List.unique(eqs);
1437 639 eqs := List.intersectionOnTrue(eqs,loopIn,intEq);
1438
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 629 times.
639 if listEmpty(eqs) then
1439 10 next := listHead(loopIn);
1440 else
1441 629 next := listHead(eqs);
1442 end if;
1443 639 rest := List.deleteMemberOnTrue(next,loopIn,intEq);
1444 639 then sortLoop(rest,m,mT,next::sortLoopIn);
1445 end match;
1446 end sortLoop;
1447
1448 protected function closePathDirectly "author:Waurich TUD 2014-01
1449 tries to close the given path with the one of the paths from the list. It outputs the whole loop"
1450 input list<Integer> pathIn;
1451 input list<list<Integer>> pathLstIn;
1452 output list<Integer> pathOut;
1453 algorithm
1454 pathOut := matchcontinue pathLstIn
1455 local
1456 Boolean closed;
1457 Integer startNode,endNode;
1458 list<Integer> path;
1459 case _
1460 algorithm
1461 // the path is already closed
1462 180 startNode := listHead(pathIn);
1463 180 endNode := List.last(pathIn);
1464
2/2
✓ Branch 0 taken 164 times.
✓ Branch 1 taken 16 times.
180 true := intEq(startNode,endNode);
1465 then
1466 pathIn;
1467 case _
1468 algorithm
1469 // it is an open path
1470
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 164 times.
164 startNode::_ := pathIn;
1471 164 endNode := List.last(pathIn);
1472 164 path := findPathByEnds(pathLstIn,startNode,endNode);
1473 164 closed := not listEmpty(path);
1474
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 136 times.
164 path := if closed then path else {};
1475 164 path := listAppend(pathIn,path);
1476 164 path := List.unique(path);
1477 then
1478 path;
1479 else
1480 algorithm
1481 ✗ Error.addInternalError("function ResolveLoops.closePathDirectly failed", sourceInfo());
1482 ✗ then
1483 fail();
1484 end matchcontinue;
1485 end closePathDirectly;
1486
1487 protected function findPathByEnds "author:Waurich TUD 2014-01
1488 searches the list<list<Integer>> for the first list<Integer> on which the start and end Node fit."
1489 input list<list<Integer>> pathLstIn;
1490 input Integer startNodeIn;
1491 input Integer endNodeIn;
1492 output list<Integer> pathOut;
1493 algorithm
1494 pathOut := matchcontinue pathLstIn
1495 local
1496 Boolean b1, b2;
1497 Integer startNode,endNode;
1498 list<Integer> path;
1499 list<list<Integer>> pathLst;
1500 case path::pathLst
1501 algorithm
1502 824 startNode := listHead(path);
1503 b1 := intEq(startNode,endNodeIn);
1504 824 endNode := List.last(path);
1505 b2 := intEq(endNode,startNodeIn);
1506
2/2
✓ Branch 0 taken 688 times.
✓ Branch 1 taken 136 times.
824 path := if not(b1 and b2) then findPathByEnds(pathLst,startNodeIn,endNodeIn) else path;
1507 then
1508 path;
1509 case {}
1510 algorithm
1511 then
1512 {};
1513 else
1514 algorithm
1515 ✗ Error.addInternalError("function ResolveLoops.findPathByEnds failed", sourceInfo());
1516 ✗ then
1517 fail();
1518 end matchcontinue;
1519 end findPathByEnds;
1520
1521 protected function countDoubleEntriesInLst "author:hkiel
1522 get the number of entries in the list which occur multiple times.
1523 Is there an entry from lstIn found in checkLst, it will be countew as doubled"
1524 input list<Integer> lstIn;
1525 output Integer num = 0;
1526 input output list<Integer> checkLst;
1527 input output list<Integer> dupLst;
1528 replaceable type ElementType subtypeof Any;
1529 algorithm
1530 ✗ for elem in lstIn loop
1531 ✗ if listMember(elem, checkLst) then
1532 ✗ num := num + 1;
1533 ✗ if not listMember(elem, dupLst) then
1534 dupLst := elem::dupLst;
1535 end if;
1536 else
1537 checkLst := elem::checkLst;
1538 end if;
1539 end for;
1540 end countDoubleEntriesInLst;
1541
1542 protected function countDoubleEntriesInLstLst "author:hkiel
1543 get the number of entries in the list of lists which occur multiple times.
1544 Is there an entry from lstIn found in checkLst, it will be countew as doubled"
1545 input list<list<Integer>> lstIn;
1546 output Integer num = 0;
1547 input output list<Integer> checkLst;
1548 input output list<Integer> dupLst;
1549 algorithm
1550
2/2
✓ Branch 0 taken 2191 times.
✓ Branch 1 taken 950 times.
3141 for lst in lstIn loop
1551
2/2
✓ Branch 0 taken 12006 times.
✓ Branch 1 taken 2191 times.
14197 for elem in lst loop
1552
2/2
✓ Branch 1 taken 3490 times.
✓ Branch 2 taken 8516 times.
12006 if listMember(elem, checkLst) then
1553 3490 num := num + 1;
1554
2/2
✓ Branch 1 taken 3458 times.
✓ Branch 2 taken 32 times.
3490 if not listMember(elem, dupLst) then
1555 dupLst := elem::dupLst;
1556 end if;
1557 else
1558 checkLst := elem::checkLst;
1559 end if;
1560 end for;
1561 end for;
1562 end countDoubleEntriesInLstLst;
1563
1564 protected function doubleEntriesInLst "author:Waurich TUD 2014-01
1565 get the entries in the list which occur multiple times.
1566 Is there an entry from lstIn found in checkLst, it will be output as doubled"
1567 input list<Integer> lstIn;
1568 output list<Integer> doubleLst = {};
1569 protected
1570 list<Integer> checkLst = {};
1571 algorithm
1572
2/2
✓ Branch 0 taken 2402 times.
✓ Branch 1 taken 326 times.
2728 for i in lstIn loop
1573
2/2
✓ Branch 1 taken 885 times.
✓ Branch 2 taken 1517 times.
2402 if listMember(i, checkLst) then
1574 doubleLst := i::doubleLst;
1575 else
1576 checkLst := i::checkLst;
1577 end if;
1578 end for;
1579 end doubleEntriesInLst;
1580
1581 protected function getPathTillNextCrossEq "author:Waurich TUD 2013-12
1582 collects the paths from the given crossEq to the next."
1583 input list<Integer> checkEqCrossNodes; //these will be traversed
1584 input BackendDAE.AdjacencyMatrix mIn;
1585 input BackendDAE.AdjacencyMatrixT mTIn;
1586 input list<Integer> allEqCrossNodes;
1587 input list<list<Integer>> unfinPathsIn;
1588 input list<list<Integer>> eqPathsIn;
1589 output list<list<Integer>> eqPathsOut = eqPathsIn;
1590 protected
1591 Integer crossEq, lastEq, prevEq;
1592 list<Integer> adjVars, nextEqs, endEqs, unfinEqs, crossNodes = checkEqCrossNodes, pathStart;
1593 list<list<Integer>> paths, adjEqs, unfinPaths = unfinPathsIn;
1594 algorithm
1595 while true loop
1596
2/2
✓ Branch 0 taken 787 times.
✓ Branch 1 taken 6553 times.
7340 if not listEmpty(unfinPaths) then
1597 787 pathStart::unfinPaths := unfinPaths;
1598 787 lastEq := listHead(pathStart);
1599 787 prevEq := List.second(pathStart);
1600 787 adjVars := arrayGet(mIn,lastEq);
1601 787 adjEqs := List.map1(adjVars,Array.getIndexFirst,mTIn);
1602
4/4
✓ Branch 0 taken 1566 times.
✓ Branch 1 taken 787 times.
✓ Branch 2 taken 1566 times.
✓ Branch 3 taken 787 times.
2353 adjEqs := list(List.deleteMemberOnTrue(lastEq, eq, intEq) for eq in adjEqs); // REMARK: this works only if there are no varCrossNodes
1603 787 adjEqs := List.filterOnFalse(adjEqs,listEmpty);
1604 787 nextEqs := List.map(adjEqs,listHead);
1605 787 (nextEqs,_) := List.deleteMemberOnTrue(prevEq,nextEqs,intEq); //do not take the path back to the previous node
1606 787 (endEqs,unfinEqs,_) := List.intersection1OnTrue(nextEqs,allEqCrossNodes,intEq);
1607 787 paths := List.map1(endEqs,cons1,pathStart); //TODO: replace this stupid cons1
1608 787 eqPathsOut := listAppend(paths,eqPathsOut) annotation(__OpenModelica_DisableListAppendWarning=true);
1609 787 paths := List.map1(unfinEqs,cons1,pathStart);
1610 787 unfinPaths := listAppend(paths,unfinPaths) annotation(__OpenModelica_DisableListAppendWarning=true);
1611 elseif not listEmpty(crossNodes) then
1612 4313 crossEq::crossNodes := crossNodes;
1613 // check the next eqNode of the crossEq whether the paths is finished here or the path goes on to another crossEq
1614 4313 adjVars := arrayGet(mIn,crossEq);
1615 4313 adjEqs := List.map1(adjVars,Array.getIndexFirst,mTIn);
1616
4/4
✓ Branch 0 taken 13195 times.
✓ Branch 1 taken 4313 times.
✓ Branch 2 taken 13195 times.
✓ Branch 3 taken 4313 times.
17508 adjEqs := list(List.deleteMemberOnTrue(crossEq, eq, intEq) for eq in adjEqs); // REMARK: this works only if there are no varCrossNodes
1617 4313 adjEqs := List.filterOnFalse(adjEqs,listEmpty);
1618 4313 nextEqs := List.flatten(adjEqs);
1619 4313 (endEqs,unfinEqs,_) := List.intersection1OnTrue(nextEqs,allEqCrossNodes,intEq);
1620 4313 paths := List.map1(endEqs,cons1,{crossEq}); //TODO: replace this stupid cons1
1621 4313 eqPathsOut := listAppend(paths,eqPathsOut) annotation(__OpenModelica_DisableListAppendWarning=true);
1622 4313 paths := List.map1(unfinEqs,cons1,{crossEq});
1623 4313 unfinPaths := listAppend(paths,unfinPaths) annotation(__OpenModelica_DisableListAppendWarning=true);
1624 else
1625 2240 return;
1626 end if;
1627 end while;
1628 end getPathTillNextCrossEq;
1629
1630 protected function cons1
1631 input Integer elem;
1632 input list<Integer> lst;
1633 output list<Integer> outLst;
1634 algorithm
1635 outLst := elem::lst;
1636 end cons1;
1637
1638 protected function replaceContractedNodes "replaces the replNodes in the pathIn with the nodeIn"
1639 input list<Integer> pathIn;
1640 input Integer nodeIn;
1641 input list<Integer> replNodes;
1642 output list<Integer> pathOut;
1643 algorithm
1644 233 pathOut := List.map2(pathIn,replaceContractedNodes2,nodeIn,replNodes);
1645 end replaceContractedNodes;
1646
1647 protected function replaceContractedNodes2 "replaces the replNodes in the pathIn with the nodeIn"
1648 input Integer entryIn;
1649 input Integer nodeIn;
1650 input list<Integer> replNodes;
1651 output Integer entryOut;
1652 protected
1653 Boolean repl;
1654 algorithm
1655 727 repl := List.isMemberOnTrue(entryIn,replNodes,intEq);
1656
2/2
✓ Branch 0 taken 565 times.
✓ Branch 1 taken 162 times.
727 entryOut := if repl then nodeIn else entryIn;
1657 end replaceContractedNodes2;
1658
1659 protected function priorizeEqsWithVarCrosses "author:Waurich TUD 2014-02
1660 the equations with the least number of varCrossNodes are the best."
1661 input list<Integer> eqsIn;
1662 input BackendDAE.AdjacencyMatrix mIn;
1663 input list<Integer> varCrossLst;
1664 output list<Integer> eqsOut;
1665 protected
1666 array<list<Integer>> priorities; //[0]eqs with no adjVarCross, [1] eqs with one adjVarCross, [2]rest
1667 algorithm
1668 70 priorities := arrayCreate(3,{});
1669
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 70 times.
310 for eq in eqsIn loop
1670 240 priorizeEqsWithVarCrosses2(eq, mIn, varCrossLst, priorities);
1671 end for;
1672 70 eqsOut := List.flatten(arrayList(priorities));
1673 end priorizeEqsWithVarCrosses;
1674
1675 protected function priorizeEqsWithVarCrosses2
1676 input Integer eq;
1677 input BackendDAE.AdjacencyMatrix mIn;
1678 input list<Integer> varCrossLst;
1679 input array<list<Integer>> priorities;
1680 protected
1681 list<Integer> eqVars,crossVars;
1682 algorithm
1683 240 eqVars := arrayGet(mIn,eq);
1684 240 crossVars := List.intersectionOnTrue(eqVars,varCrossLst,intEq);
1685
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 180 times.
240 if listEmpty(crossVars) then
1686 60 arrayGetAppendLst(1,{eq},priorities);
1687 elseif List.hasOneElement(crossVars) then
1688 140 arrayGetAppendLst(2,{eq},priorities);
1689 else
1690 40 arrayGetAppendLst(3,{eq},priorities);
1691 end if;
1692 end priorizeEqsWithVarCrosses2;
1693
1694 protected function evaluateLoop
1695 input list<Integer> loopIn;
1696 input tuple<BackendDAE.AdjacencyMatrix,BackendDAE.AdjacencyMatrixT,list<Integer>> tplIn;
1697 output Boolean resolve = true;
1698 protected
1699 Boolean r1,r2;
1700 Integer numInLoop,numOutLoop;
1701 list<Integer> eqCrossLst,chk={},dup={};
1702 list<list<Integer>> eqVars;
1703 BackendDAE.AdjacencyMatrix m;
1704 algorithm
1705
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 950 times.
950 if not intEq(Flags.getConfigInt(Flags.RESHUFFLE),3) then
1706 950 (m,_,eqCrossLst) := tplIn;
1707 950 eqVars := List.map1(loopIn,Array.getIndexFirst,m);
1708
1709 // check if its worth to resolve the loop. Therefore compare the amount of vars in and outside the loop
1710 950 (numInLoop, chk, dup) := countDoubleEntriesInLstLst(eqVars, chk, dup);
1711 950 numOutLoop := listLength(chk)-listLength(dup);
1712
4/4
✓ Branch 0 taken 488 times.
✓ Branch 1 taken 462 times.
✓ Branch 2 taken 15 times.
✓ Branch 3 taken 473 times.
950 r1 := intGe(numInLoop,numOutLoop-1) and intLe(numInLoop,6);
1713 950 r2 := intGe(numInLoop,numOutLoop-2);
1714
1/2
✓ Branch 1 taken 950 times.
✗ Branch 2 not taken.
950 r1 := if intEq(Flags.getConfigInt(Flags.RESHUFFLE),1) then r1 else false;
1715
1/2
✓ Branch 1 taken 950 times.
✗ Branch 2 not taken.
950 resolve := if intEq(Flags.getConfigInt(Flags.RESHUFFLE),2) then r2 else r1;
1716 end if;
1717 end evaluateLoop;
1718
1719 protected uniontype TripleLoopInfo
1720 record TRIPLE_LOOP_INFO
1721 BackendDAE.AdjacencyMatrix m;
1722 list<Integer> mapIndices;
1723 BackendDAE.AdjacencyMatrix map;
1724 array<Integer> count "occurrences of each variable in the merged nodes' rows";
1725 Integer entries, distinct, singles "of those rows";
1726 end TRIPLE_LOOP_INFO;
1727 end TripleLoopInfo;
1728
1729 protected function tripleLoopInfo
1730 "Every triple loop is evaluated together with the rows of all merged nodes,
1731 so those are counted once."
1732 input BackendDAE.AdjacencyMatrix m;
1733 input Integer numVars;
1734 input list<Integer> mapIndices;
1735 input BackendDAE.AdjacencyMatrix map;
1736 output TripleLoopInfo info;
1737 protected
1738 array<Integer> count = arrayCreate(numVars, 0);
1739 Integer entries = 0, distinct = 0, singles = 0;
1740 algorithm
1741
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 27 times.
35 for i in mapIndices loop
1742
2/2
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 8 times.
22 for j in arrayGet(map,i) loop
1743 14 (entries, distinct, singles) := countRow(arrayGet(m,j), count, entries, distinct, singles, 1);
1744 end for;
1745 end for;
1746 27 info := TRIPLE_LOOP_INFO(m, mapIndices, map, count, entries, distinct, singles);
1747 end tripleLoopInfo;
1748
1749 protected function countRow
1750 input list<Integer> row;
1751 input array<Integer> count;
1752 input output Integer entries, distinct, singles;
1753 input Integer delta "1 to add the row, -1 to remove it again";
1754 protected
1755 Integer c;
1756 algorithm
1757
2/2
✓ Branch 0 taken 1141 times.
✓ Branch 1 taken 338 times.
1479 for v in row loop
1758 1141 c := count[v] + delta;
1759 1141 arrayUpdate(count, v, c);
1760 1141 entries := entries + delta;
1761
2/2
✓ Branch 0 taken 601 times.
✓ Branch 1 taken 540 times.
1141 if delta > 0 then
1762
2/2
✓ Branch 0 taken 431 times.
✓ Branch 1 taken 170 times.
601 if c == 1 then
1763 431 distinct := distinct + 1;
1764 431 singles := singles + 1;
1765 elseif c == 2 then
1766 170 singles := singles - 1;
1767 end if;
1768 else
1769
2/2
✓ Branch 0 taken 378 times.
✓ Branch 1 taken 162 times.
540 if c == 0 then
1770 378 distinct := distinct - 1;
1771 378 singles := singles - 1;
1772 elseif c == 1 then
1773 162 singles := singles + 1;
1774 end if;
1775 end if;
1776 end for;
1777 end countRow;
1778
1779 protected function evaluateTripleLoop
1780 "author:kabdelhak FHB 2019-07
1781 Special case for loops containing three eqCrossNodes"
1782 input list<Integer> loopIn;
1783 input TripleLoopInfo info;
1784 output Boolean resolve = true;
1785 protected
1786 Boolean r1,r2;
1787 Integer entries, distinct, singles, numInLoop, numOutLoop;
1788 algorithm
1789
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 54 times.
54 if not intEq(Flags.getConfigInt(Flags.RESHUFFLE),3) then
1790 54 (entries, distinct, singles) := (info.entries, info.distinct, info.singles);
1791
2/2
✓ Branch 0 taken 162 times.
✓ Branch 1 taken 54 times.
216 for j in loopIn loop
1792 162 (entries, distinct, singles) := countRow(arrayGet(info.m,j), info.count, entries, distinct, singles, 1);
1793 end for;
1794
2/2
✓ Branch 0 taken 162 times.
✓ Branch 1 taken 54 times.
216 for j in loopIn loop
1795 162 countRow(arrayGet(info.m,j), info.count, 0, 0, 0, -1);
1796 end for;
1797 54 numInLoop := entries - distinct;
1798 54 numOutLoop := singles;
1799
1800 // check if its worth to resolve the loop. Therefore compare the amount of vars in and outside the loop
1801
3/4
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 46 times.
54 r1 := intGe(numInLoop,numOutLoop-1) and intLe(numInLoop,10);
1802 54 r2 := intGe(numInLoop,numOutLoop-2);
1803
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
54 r1 := if intEq(Flags.getConfigInt(Flags.RESHUFFLE),1) then r1 else false;
1804
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
54 resolve := if intEq(Flags.getConfigInt(Flags.RESHUFFLE),2) then r2 else r1;
1805 end if;
1806 end evaluateTripleLoop;
1807
1808 protected function updateTripleLoop
1809 "author:kabdelhak FHB 2019-07
1810 Update function for special case including for loops containing three eqCrossNodes"
1811 input output list<Integer> loopFull;
1812 input TripleLoopInfo info;
1813 algorithm
1814
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
46 for i in info.mapIndices loop
1815 ✗ loopFull := listAppend(arrayGet(info.map,i),loopFull);
1816 end for;
1817 end updateTripleLoop;
1818
1819 protected function sumUp2Equations "author:Waurich TUD 2013-12
1820 sums up or subtracts 2 equations, depending on the boolean (true=+, false =-)"
1821 input Boolean sumUp;
1822 input BackendDAE.Equation eq1;
1823 input BackendDAE.Equation eq2;
1824 output BackendDAE.Equation eqOut;
1825 protected
1826 DAE.Exp exp1, exp2, exp3, exp4;
1827 algorithm
1828
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 639 times.
639 BackendDAE.EQUATION(exp=exp1, scalar=exp2) := eq1;
1829
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 639 times.
639 BackendDAE.EQUATION(exp=exp3, scalar=exp4) := eq2;
1830 639 exp1 := sumUp2Expressions(sumUp, exp1, exp3);
1831 639 exp2 := sumUp2Expressions(sumUp, exp2, exp4);
1832 639 exp2 := sumUp2Expressions(false, exp2, exp1);
1833 639 (exp2, _) := ExpressionSimplify.simplify(exp2);
1834 639 exp1 := Expression.createZeroExpression(Expression.typeof(exp2));
1835 639 eqOut := BackendDAE.EQUATION(exp1, exp2, DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN);
1836 639 eqOut := simplifyZeroAssignment(eqOut);
1837 end sumUp2Equations;
1838
1839 protected function simplifyZeroAssignment
1840 input BackendDAE.Equation eIn;
1841 output BackendDAE.Equation eOut;
1842 algorithm
1843 eOut := match eIn
1844 local
1845 DAE.Exp e;
1846 DAE.ElementSource source;
1847 BackendDAE.EquationAttributes attr;
1848 case BackendDAE.EQUATION(exp = DAE.RCONST(0.0), scalar = DAE.BINARY(exp1 = DAE.RCONST(_), operator = DAE.MUL(), exp2 = e as DAE.CREF()), source=source, attr=attr)
1849 10 then BackendDAE.EQUATION(DAE.RCONST(0.0),e,source, attr);
1850 case BackendDAE.EQUATION(scalar = DAE.RCONST(0.0), exp = DAE.BINARY(exp1 = DAE.RCONST(_), operator = DAE.MUL(), exp2 = e as DAE.CREF()), source=source, attr=attr)
1851 ✗ then BackendDAE.EQUATION(DAE.RCONST(0.0),e,source, attr);
1852 else
1853 then eIn;
1854 end match;
1855 end simplifyZeroAssignment;
1856
1857 protected function CRefIsPosOnRHS "author:Waurich TUD 2013-12
1858 checks if the given cref occurs with a positiv algebraic sign on the right
1859 hand side of the equation. if its on the left hand side the algebraic sign
1860 has to be negated."
1861 input DAE.ComponentRef crefIn;
1862 input BackendDAE.Equation eqIn;
1863 output Boolean isPos;
1864 algorithm
1865 isPos := matchcontinue eqIn
1866 local
1867 Boolean exists1, sign1, sign2;
1868 DAE.Exp e1, e2;
1869
1870 case BackendDAE.EQUATION(exp=e1, scalar=e2) algorithm
1871 2488 (exists1, sign1) := expIsCref(e1, crefIn);
1872 2488 (_, sign2) := expIsCref(e2, crefIn);
1873
2/2
✓ Branch 0 taken 738 times.
✓ Branch 1 taken 1750 times.
2488 sign1 := if exists1 then not sign1 else sign2;
1874 then sign1;
1875
1876 else algorithm
1877 ✗ print("add a case to CRefIsPosOnRHS"+BackendDump.equationString(eqIn)+"\n");
1878 ✗ then fail();
1879 end matchcontinue;
1880 end CRefIsPosOnRHS;
1881
1882 protected function expIsCref "author: Waurich TUD 2013-12
1883 checks if the cref is in the exp.
1884 if it occurs with a plus sign then true if its with a minus sign then false."
1885 input DAE.Exp expIn;
1886 input DAE.ComponentRef crefIn;
1887 output Boolean isInExp;
1888 output Boolean algSign;
1889 algorithm
1890 (isInExp,algSign) := match expIn
1891 local
1892 Real r;
1893 Boolean sameCref,sign, sign1, sign2, exists, exists1, exists2;
1894 DAE.ComponentRef cref;
1895 DAE.Exp exp1, exp2;
1896 case DAE.CREF(componentRef=cref)
1897 algorithm
1898 // just a cref
1899 9097 sameCref := ComponentReferenceBasics.crefEqualNoStringCompare(crefIn,cref);
1900 then
1901 (sameCref,true);
1902 case DAE.BINARY(exp1=exp1, operator = DAE.SUB(), exp2=exp2)
1903 algorithm
1904 //exp1-exp2
1905 2790 (exists1,sign1) := expIsCref(exp1,crefIn);
1906 2790 (exists2,sign2) := expIsCref(exp2,crefIn);
1907 2790 sign2 := boolNot(sign2);
1908 2790 exists := boolOr(exists1,exists2);
1909
4/4
✓ Branch 0 taken 1148 times.
✓ Branch 1 taken 1642 times.
✓ Branch 2 taken 687 times.
✓ Branch 3 taken 461 times.
2790 sign := exists1 and sign1;
1910
2/2
✓ Branch 0 taken 2030 times.
✓ Branch 1 taken 760 times.
2790 sign := if exists2 then sign2 else sign;
1911 then
1912 (exists,sign);
1913 case DAE.BINARY(exp1=exp1, operator = DAE.ADD(), exp2=exp2)
1914 algorithm
1915 //exp1+exp2
1916 2900 (exists1,sign1) := expIsCref(exp1,crefIn);
1917 2900 (exists2,sign2) := expIsCref(exp2,crefIn);
1918 2900 exists := boolOr(exists1,exists2);
1919
3/4
✓ Branch 0 taken 900 times.
✓ Branch 1 taken 2000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 900 times.
2900 sign := exists1 and sign1;
1920
2/2
✓ Branch 0 taken 1486 times.
✓ Branch 1 taken 1414 times.
2900 sign := if exists2 then sign2 else sign;
1921 then
1922 (exists,sign);
1923 case DAE.BINARY(exp1=exp1 as DAE.CREF(), operator = DAE.MUL(), exp2=DAE.RCONST(r))
1924 algorithm
1925 //exp1*rconst
1926 ✗ (exists,_) := expIsCref(exp1,crefIn);
1927 ✗ sign := r > 0;
1928 then
1929 (exists,sign);
1930 case DAE.BINARY(exp1=exp1 as DAE.RCONST(r), operator = DAE.MUL(), exp2=DAE.CREF())
1931 algorithm
1932 //rconst*exp2
1933 90 (exists,_) := expIsCref(exp1,crefIn);
1934 90 sign := r > 0;
1935 then
1936 (exists,sign);
1937 case DAE.UNARY(operator=DAE.UMINUS(),exp=exp1)
1938 algorithm
1939 // -(exp)
1940 940 (exists,sign) := expIsCref(exp1,crefIn);
1941 940 sign := boolNot(sign);
1942 then
1943 (exists,sign);
1944 case DAE.RCONST()
1945 algorithm
1946 // constant
1947 then
1948 (false,false);
1949 case DAE.ICONST()
1950 algorithm
1951 // constant
1952 then
1953 (false,false);
1954 else
1955 algorithm
1956 ✗ print("add a case to expIsCref:"+ExpressionBasics.printExpStr(expIn)+"\n");
1957 then
1958 (false,false);
1959 end match;
1960 end expIsCref;
1961
1962 protected function listLengthIs
1963 input list<Integer> lst;
1964 input Integer value;
1965 output Boolean bOut;
1966 algorithm
1967 360 bOut := intEq(listLength(lst),value);
1968 end listLengthIs;
1969
1970 public function partitionBipartiteGraph "author: Waurich TUD 2013-12
1971 checks if there are independent subgraphs in the BIPARTITE graph. the given
1972 indeces refer to the equation indeces (rows in the adjacencyMatrix)."
1973 input BackendDAE.AdjacencyMatrix m;
1974 input BackendDAE.AdjacencyMatrixT mT;
1975 output list<list<Integer>> partitions;
1976 protected
1977 Integer numEqs, numVars;
1978 array<Integer> markEqs, markVars;
1979 algorithm
1980 numEqs := arrayLength(m);
1981 numVars := arrayLength(mT);
1982
1983
2/2
✓ Branch 0 taken 2568 times.
✓ Branch 1 taken 2140 times.
4708 if numEqs == 0 or numVars == 0 then
1984 2568 partitions := {{}};
1985 else
1986 2140 markEqs := arrayCreate(numEqs,-1);
1987 2140 markVars := arrayCreate(numVars,-1);
1988 2140 (_,partitions) := colorNodePartitions(m,mT,{1},markEqs,markVars,1,{});
1989 end if;
1990 end partitionBipartiteGraph;
1991
1992 protected function colorNodePartitions "author:Waurich TUD 2013-12
1993 helper for partitionsGraph1. Traverse the graph in a BFS manner.
1994 mark all visited nodes, gather partitions, color mark-arrays"
1995 input BackendDAE.AdjacencyMatrix m;
1996 input BackendDAE.AdjacencyMatrixT mT;
1997 input list<Integer> checkNextIn;
1998 input array<Integer> markEqs;
1999 input array<Integer> markVars;
2000 input Integer currNumberIn;
2001 input list<list<Integer>> partitionsIn;
2002 input Integer nextIndex = 1;
2003 output Integer currNumberOut;
2004 output list<list<Integer>> partitionsOut;
2005 protected
2006 Integer eq, next_index;
2007 list<Integer> rest, vars, eqs, part;
2008 list<list<Integer>> restPart, partitions;
2009 algorithm
2010 (currNumberOut,partitionsOut) := match checkNextIn
2011 //found no unassigned eqnode
2012 case {0}
2013
1/2
✓ Branch 0 taken 2140 times.
✗ Branch 1 not taken.
2140 then (currNumberIn - 1, partitionsIn);
2014
2015 case eq::rest
2016 algorithm
2017 //check unassigned node
2018
2/2
✓ Branch 1 taken 44184 times.
✓ Branch 2 taken 4761 times.
48945 if arrayGetIsNotPositive(eq,markEqs) then
2019 //mark this eq and add to partition
2020 44184 arrayUpdate(markEqs, eq, currNumberIn);
2021
2/2
✓ Branch 0 taken 2140 times.
✓ Branch 1 taken 42044 times.
44184 if listEmpty(partitionsIn) then
2022 partitions := {{eq}};
2023 else
2024 42044 part::restPart := partitionsIn;
2025 part := eq::part;
2026 partitions := part::restPart;
2027 end if;
2028
2029 // get adjacent equation nodes
2030 44184 vars := arrayGet(m,eq);
2031
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44184 times.
44184 true := not listEmpty(vars);
2032
2033 //all vars that havent been traversed
2034 44184 vars := List.filter1OnTrue(vars,arrayGetIsNotPositive,markVars);
2035 44184 List.map2_0(vars,Array.updateIndexFirst,currNumberIn,markVars);
2036
2037 //all eqs that havent been traversed
2038 44184 eqs := List.fold1(vars,getArrayEntryAndAppend,mT,{});
2039 44184 eqs := List.filter1OnTrue(eqs,arrayGetIsNegative,markEqs); // all new equations which havent been queued
2040 44184 List.map2_0(eqs,Array.updateIndexFirst,0,markEqs);
2041
2042 // check them later
2043 44184 rest := listAppend(rest,eqs) annotation(__OpenModelica_DisableListAppendWarning=true);
2044 else
2045 //the node has been investigated already
2046 partitions := partitionsIn;
2047 end if;
2048 48945 then
2049 colorNodePartitions(m,mT,rest,markEqs,markVars,currNumberIn,partitions,nextIndex);
2050
2051 case {}
2052 algorithm
2053 //nothing left in this partition
2054 eq := 0;
2055 next_index := nextIndex;
2056
2057 // Search for the next unmarked equation, starting from the next unsearched index.
2058
2/2
✓ Branch 0 taken 12254 times.
✓ Branch 1 taken 527 times.
46324 for i in nextIndex:arrayLength(markEqs) loop
2059
2/2
✓ Branch 1 taken 10641 times.
✓ Branch 2 taken 33543 times.
44184 if markEqs[i] == -1 then
2060 eq := i;
2061 10641 next_index := i + 1;
2062 10641 break;
2063 end if;
2064 end for;
2065 12781 then
2066 colorNodePartitions(m,mT,{eq},markEqs,markVars,currNumberIn+1,{}::partitionsIn, next_index);
2067
2068 end match;
2069 end colorNodePartitions;
2070
2071 protected function arrayGetIsNotPositive" outputs true if the indexed entry is not zero."
2072 input Integer idx;
2073 input array<Integer> arrayIn;
2074 output Boolean isNonZero;
2075 algorithm
2076 168397 isNonZero := arrayGet(arrayIn,idx) <= 0;
2077 end arrayGetIsNotPositive;
2078
2079 protected function arrayGetIsNegative" outputs true if the indexed entry is negative."
2080 input Integer idx;
2081 input array<Integer> arrayIn;
2082 output Boolean isNonZero;
2083 algorithm
2084 119452 isNonZero := arrayGet(arrayIn,idx) < 0;
2085 end arrayGetIsNegative;
2086
2087 protected function getArrayEntryAndAppend
2088 input Integer entry;
2089 input BackendDAE.AdjacencyMatrixT m;
2090 input list<Integer> lstIn;
2091 output list<Integer> lstOut;
2092 protected
2093 list<Integer> lst;
2094 algorithm
2095 73071 lst := arrayGet(m,entry);
2096 73071 lstOut := listAppend(lst,lstIn);
2097 end getArrayEntryAndAppend;
2098
2099 protected function gatherCrossNodes "author: Waurich TUD 2014-02
2100 checks if the indexed node has more than 2 neighbours (its a crossroad)."
2101 input Integer idx;
2102 input BackendDAE.AdjacencyMatrix m;
2103 input BackendDAE.AdjacencyMatrix mT;
2104 input list<Integer> lstIn;
2105 output list<Integer> lstOut;
2106 protected
2107 Boolean isCross;
2108 Integer num;
2109 list<Integer> row;
2110 algorithm
2111 // the node has more than 2 edges, it might be a crossnode
2112 55794 row := arrayGet(m,idx);
2113 55794 num := listLength(row);
2114 isCross := intGt(num,2);
2115
2/2
✓ Branch 0 taken 24976 times.
✓ Branch 1 taken 30818 times.
55794 lstOut := if isCross then idx::lstIn else lstIn;
2116 end gatherCrossNodes;
2117
2118 protected function isAddOrSubExp
2119 input DAE.Exp inExp;
2120 input tuple<Boolean,BackendDAE.Variables> inTuple;
2121 output DAE.Exp outExp;
2122 output tuple<Boolean,BackendDAE.Variables> outTuple;
2123 algorithm
2124 (outExp,outTuple) := match(inExp,inTuple)
2125 local
2126 Boolean b;
2127 BackendDAE.Variables vars;
2128 DAE.Exp exp,exp1,exp2;
2129 DAE.ComponentRef cref;
2130 case (DAE.CREF(componentRef=cref),(true,vars))
2131 algorithm
2132 //x, but reject array elements with non-constant indices since
2133 //resolveLoops cannot map them to a single scalar variable
2134 61188 b := Expression.subscriptConstants(ComponentReferenceBasics.crefSubs(cref));
2135
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 61177 times.
61199 then (inExp,(b,vars));
2136
2137 case (DAE.UNARY(exp=exp1),(true,vars))
2138 algorithm
2139 // (-x)
2140 2874 (_,(b,_)) := isAddOrSubExp(exp1,(true,vars));
2141
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 2835 times.
2913 then (inExp,(b,vars));
2142
2143 case (DAE.RCONST(),(true,vars)) // maybe we have to remove this, because this is just for kirchhoffs current law
2144 algorithm
2145 //const.
2146 5124 then (inExp,(true,vars));
2147
2148 case (DAE.BINARY(exp1 = exp1,operator = DAE.ADD(),exp2 = exp2),(true,vars))
2149 algorithm
2150 //x + y
2151 15412 (_,(b,_)) := isAddOrSubExp(exp1,(true,vars));
2152
2/2
✓ Branch 0 taken 4681 times.
✓ Branch 1 taken 10731 times.
20093 (_,(b,_)) := isAddOrSubExp(exp2,(b,vars));
2153
2/2
✓ Branch 0 taken 7689 times.
✓ Branch 1 taken 7723 times.
23101 then (inExp,(b,vars));
2154
2155 case (DAE.BINARY(exp1=exp1,operator = DAE.SUB(),exp2=exp2),(true,vars))
2156 algorithm
2157 //x - y
2158 8580 (_,(b,_)) := isAddOrSubExp(exp1,(true,vars));
2159
2/2
✓ Branch 0 taken 1294 times.
✓ Branch 1 taken 7286 times.
9874 (_,(b,_)) := isAddOrSubExp(exp2,(b,vars));
2160
2/2
✓ Branch 0 taken 1741 times.
✓ Branch 1 taken 6839 times.
10321 then (inExp,(b,vars));
2161
2162 case (DAE.BINARY(exp1=DAE.CREF(componentRef=cref),operator = DAE.MUL(),exp2=exp2),(true,vars))
2163 algorithm
2164 //state*const. is allowed
2165
4/4
✓ Branch 1 taken 442 times.
✓ Branch 2 taken 7978 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 439 times.
8420 b := BackendVariable.isState(cref, vars) and Expression.isConst(exp2);
2166 8420 then (inExp,(b,vars));
2167
2168 case (DAE.BINARY(exp1 = exp1,operator = DAE.MUL(),exp2=DAE.CREF(componentRef=cref)),(true,vars))
2169 algorithm
2170 //const*state. is allowed
2171
4/4
✓ Branch 1 taken 3367 times.
✓ Branch 2 taken 2284 times.
✓ Branch 4 taken 2508 times.
✓ Branch 5 taken 859 times.
5651 b := Expression.isConst(exp1) and BackendVariable.isState(cref, vars);
2172 5651 then (inExp,(b,vars));
2173 else
2174 algorithm
2175 38861 then
2176 (inExp,(false,Util.tuple22(inTuple)));
2177 end match;
2178 end isAddOrSubExp;
2179
2180 protected function sumUp2Expressions "author:Waurich TUD 2013-12
2181 sums up or subtracts 2 expressions, depending on the boolen (true=+, false =-)"
2182 input Boolean sumUp;
2183 input DAE.Exp exp1;
2184 input DAE.Exp exp2;
2185 output DAE.Exp expOut;
2186 protected
2187 DAE.Operator op;
2188 DAE.Type ty;
2189 algorithm
2190 ty := DAE.T_REAL_DEFAULT;
2191
2/2
✓ Branch 0 taken 820 times.
✓ Branch 1 taken 1097 times.
1917 op := if sumUp then DAE.ADD(ty) else DAE.SUB(ty);
2192 1917 expOut := DAE.BINARY(exp1,op,exp2);
2193 1917 (expOut,_) := ExpressionSimplify.simplify(expOut);
2194 end sumUp2Expressions;
2195
2196 protected function intLstIsEqual
2197 input list<Integer> lst1;
2198 input list<Integer> lst2;
2199 output Boolean bOut;
2200 algorithm
2201 1741 bOut := List.isEqualOnTrue(lst1,lst2,intEq);
2202 end intLstIsEqual;
2203
2204 protected function sortPathsAsChain "author: Waurich TUD 2014-01
2205 sorts the paths, so that the endNode of the next Path is an endNode of one of
2206 all already sorted path.
2207 the contractedNodes represent the endNodes of the already sorted path"
2208 input list<list<Integer>> pathsIn;
2209 output list<list<Integer>> pathsOut;
2210 algorithm
2211 pathsOut := matchcontinue pathsIn
2212 local
2213 list<list<Integer>> pathLst;
2214 case {}
2215 then
2216 {};
2217 case _
2218 algorithm
2219 532 pathLst := sortPathsAsChain1(pathsIn,0,0,{});
2220 then
2221 pathLst;
2222 else
2223 algorithm
2224 then
2225 pathsIn;
2226 end matchcontinue;
2227 end sortPathsAsChain;
2228
2229 protected function sortPathsAsChain1 "author: Waurich TUD 2014-01
2230 sorts the paths, so that the endNode of the next Path is an endNode of one of
2231 all already sorted path.
2232 the contractedNodes represent the endNodes of the already sorted path"
2233 input list<list<Integer>> pathsIn;
2234 input Integer firstNode;
2235 input Integer lastNode;
2236 input list<list<Integer>> sortedPathsIn;
2237 output list<list<Integer>> sortedPathsOut;
2238 algorithm
2239 sortedPathsOut := matchcontinue(pathsIn,firstNode,lastNode,sortedPathsIn)
2240 local
2241 Integer startNode,endNode;
2242 list<Integer> path;
2243 list<list<Integer>> rest, paths1, paths2, allPaths, sortedPaths;
2244 case({},_,_,_)
2245 algorithm
2246 then
2247 sortedPathsIn;
2248 case(_,-1,-1,_)
2249 algorithm
2250 then
2251 sortedPathsIn;
2252 case(path::rest,_,_,{})
2253 algorithm
2254 // the first node
2255 532 startNode := listHead(path);
2256 532 endNode := List.last(path);
2257 532 sortedPaths := sortPathsAsChain1(rest,startNode,endNode,{path});
2258 then
2259 sortedPaths;
2260 case(_,_,_,_)
2261 algorithm
2262 // check if theres a path that continues the endNode
2263 120 paths1 := List.filter1OnTrue(pathsIn, firstInListIsEqual, lastNode);
2264 120 paths2 := List.filter1OnTrue(pathsIn, lastInListIsEqual, lastNode);
2265 120 allPaths := listAppend(paths1,paths2);
2266
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 87 times.
120 false := listEmpty(allPaths);
2267 87 path := listHead(allPaths);
2268
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
87 endNode := if not listEmpty(allPaths) then List.last(path) else -1;
2269
2/2
✓ Branch 0 taken 67 times.
✓ Branch 1 taken 20 times.
87 endNode := if not listEmpty(paths2) then listHead(path) else -1;
2270 87 rest := List.deleteMemberOnTrue(path,pathsIn,function List.isEqualOnTrue(inCompFunc = intEq));
2271 87 sortedPaths := listAppend(sortedPathsIn,{path});
2272 87 sortedPaths := sortPathsAsChain1(rest,firstNode,endNode,sortedPaths);
2273
2274 then
2275 sortedPaths;
2276 case(_,_,_,_)
2277 algorithm
2278 // check if theres a path that continues the startNode
2279 33 paths1 := List.filter1OnTrue(pathsIn, firstInListIsEqual, firstNode);
2280 33 paths2 := List.filter1OnTrue(pathsIn, lastInListIsEqual, firstNode);
2281 33 allPaths := listAppend(paths1,paths2);
2282
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 6 times.
33 false := listEmpty(allPaths);
2283 6 path := listHead(allPaths);
2284
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
6 startNode := if not listEmpty(allPaths) then List.last(path) else -1;
2285
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 5 times.
6 startNode := if not listEmpty(paths2) then listHead(path) else -1;
2286 6 rest := List.deleteMemberOnTrue(path,pathsIn,function List.isEqualOnTrue(inCompFunc = intEq));
2287 sortedPaths := path::sortedPathsIn;
2288 6 sortedPaths := sortPathsAsChain1(rest,startNode,lastNode,sortedPaths);
2289 then
2290 sortedPaths;
2291 else
2292 algorithm// TODO: this case just put another, unconnectable path to the front of the list.
2293 //this path is currently only appendable through the startNode but it has to be also appendable throught the endNode.
2294 //this might have no effect because in those partitions is either one long path or none.
2295
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
27 path::rest := pathsIn;
2296 sortedPaths := path::sortedPathsIn;
2297 27 startNode := listHead(path);
2298 27 sortedPaths := sortPathsAsChain1(rest,startNode,lastNode,sortedPaths);
2299 then
2300 sortedPaths;
2301 end matchcontinue;
2302 end sortPathsAsChain1;
2303
2304 protected function firstInListIsEqual "author:Waurich TUD 2014-01
2305 checks if the first element in a list is equal to the given value"
2306 input list<Integer> lstIn;
2307 input Integer value;
2308 output Boolean isEq;
2309 protected
2310 Integer first;
2311 algorithm
2312 1623 first := listHead(lstIn);
2313 1623 isEq := intEq(first,value);
2314 end firstInListIsEqual;
2315
2316 protected function lastInListIsEqual "author:Waurich TUD 2014-01
2317 checks if the last element in a list is equal to the given value"
2318 input list<Integer> lstIn;
2319 input Integer value;
2320 output Boolean isEq;
2321 protected
2322 Integer last;
2323 algorithm
2324 1234 last := List.last(lstIn);
2325 1234 isEq := intEq(last,value);
2326 end lastInListIsEqual;
2327
2328 protected function connect2PathsToLoops "author:Waurich TUD 2014-01
2329 connects 2 paths to a closed loop"
2330 input list<list<Integer>> pathsIn;
2331 input list<list<Integer>> loopsIn; //empty input
2332 input list<list<Integer>> restPathsIn; // empt input
2333 output list<list<Integer>> pathsOut = loopsIn;
2334 output list<list<Integer>> restPathsOut = restPathsIn;
2335 protected
2336 Boolean closedALoop;
2337 Integer startNode, endNode;
2338 list<Integer> path;
2339 list<list<Integer>> rest = pathsIn, endPaths, startPaths, newLoops;
2340 algorithm
2341
2/2
✓ Branch 0 taken 81 times.
✓ Branch 1 taken 451 times.
532 if listEmpty(pathsIn) then
2342 pathsOut := {};
2343 restPathsOut := {};
2344 451 return;
2345 end if;
2346
2347
2/2
✓ Branch 0 taken 154 times.
✓ Branch 1 taken 81 times.
235 while not listEmpty(rest) loop
2348 154 path::rest := rest;
2349 154 startNode := listHead(path);
2350 154 endNode := List.last(path);
2351
2352
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 154 times.
154 if intEq(startNode,endNode) then
2353 // the loop closes itself
2354 pathsOut := path::pathsOut;
2355 elseif listEmpty(rest) then
2356 restPathsOut := path::restPathsOut;
2357 else
2358 // check if there is another path that closes the Loop. if not: put the path to the restPaths
2359 73 startPaths := List.filter1OnTrue(rest,firstInListIsEqual,startNode);
2360 73 startPaths := List.filter1OnTrue(startPaths,lastInListIsEqual,endNode);
2361 73 endPaths := List.filter1OnTrue(rest,firstInListIsEqual,endNode);
2362 73 endPaths := List.filter1OnTrue(endPaths,lastInListIsEqual,startNode);
2363 73 endPaths := listAppend(startPaths,endPaths);
2364 73 closedALoop := not listEmpty(endPaths);
2365
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 70 times.
73 newLoops := if closedALoop then connectPaths(path,endPaths) else {};
2366
2/2
✓ Branch 0 taken 70 times.
✓ Branch 1 taken 3 times.
73 if not closedALoop then
2367 restPathsOut := path::restPathsOut;
2368 end if;
2369 73 pathsOut := listAppend(newLoops,pathsOut);
2370 end if;
2371 end while;
2372 end connect2PathsToLoops;
2373
2374 protected function connectPaths "author:Waurich TUD 2014-02
2375 connects a given paths with the closing paths i.e. delete the first and last
2376 node of the path and append it to the given paths"
2377 input list<Integer> pathIn;
2378 input list<list<Integer>> closingPaths;
2379 output list<list<Integer>> loopsOut;
2380 protected
2381 list<Integer> path;
2382 algorithm
2383
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 _::path := pathIn;
2384 3 path := List.stripLast(path);
2385 3 loopsOut := List.map1(closingPaths,listAppend,path);
2386 end connectPaths;
2387
2388
2389 //____________________________________________________
2390 //reshuffle systems of equations, not yet finished
2391 //____________________________________________________
2392
2393 public function reshuffling_post
2394 input BackendDAE.BackendDAE inDAE;
2395 output BackendDAE.BackendDAE outDAE;
2396 protected
2397 BackendDAE.EqSystems eqSystems;
2398 algorithm
2399
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
1 if Flags.isSet(Flags.RESHUFFLE_POST) then
2400 //print("RESHUFFLING\n");
2401 //BackendDump.dumpBackendDAE(inDAE,"INDAE");
2402 1 eqSystems := List.map1(inDAE.eqs,reshuffling_post0, inDAE.shared);
2403 1 outDAE := BackendDAE.DAE(eqSystems, inDAE.shared);
2404 //BackendDump.dumpBackendDAE(outDAE,"OUTDAE");
2405 else
2406 outDAE := inDAE;
2407 end if;
2408 end reshuffling_post;
2409
2410 protected function reshuffling_post0 "author: waurich TUD 2014-09"
2411 input BackendDAE.EqSystem isyst;
2412 input BackendDAE.Shared shared;
2413 output BackendDAE.EqSystem osyst;
2414 protected
2415 BackendDAE.StrongComponents comps;
2416 algorithm
2417
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
4 BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(comps=comps)):=isyst;
2418 4 osyst := List.fold1(comps,reshuffling_post1,shared,isyst);
2419 end reshuffling_post0;
2420
2421 protected function reshuffling_post1
2422 input BackendDAE.StrongComponent compIn;
2423 input BackendDAE.Shared shared;
2424 input BackendDAE.EqSystem systIn;
2425 output BackendDAE.EqSystem systOut;
2426 algorithm
2427 systOut := matchcontinue compIn
2428 local
2429 list<Integer> vIdcs,eqIdcs;
2430 BackendDAE.EqSystem eqSys;
2431 BackendDAE.JacobianType jacType;
2432 Option<list<tuple<Integer, Integer, BackendDAE.Equation>>> ojac;
2433 case BackendDAE.EQUATIONSYSTEM(eqns=eqIdcs, vars=vIdcs, jac=BackendDAE.FULL_JACOBIAN(ojac), jacType=jacType as BackendDAE.JAC_LINEAR())
2434 algorithm
2435 1 (eqSys,_) := reshuffling_post2(eqIdcs, vIdcs, systIn, shared, ojac, jacType);
2436 then eqSys;
2437 else
2438 then systIn;
2439 end matchcontinue;
2440 end reshuffling_post1;
2441
2442 protected function reshuffling_post2 ""
2443 input list<Integer> eqIdcs;
2444 input list<Integer> varIdcs;
2445 input BackendDAE.EqSystem dae;
2446 input BackendDAE.Shared shared;
2447 input Option<list<tuple<Integer, Integer, BackendDAE.Equation>>> ojac;
2448 input BackendDAE.JacobianType jacType;
2449 output BackendDAE.EqSystem daeOut;
2450 output Boolean outRunMatching;
2451 protected
2452 Integer size;
2453 list<list<Integer>> resEqs;
2454 array<Integer> ass1, ass2, ass1Sys, ass2Sys;
2455 list<tuple<Boolean,String>> varAtts,eqAtts;
2456 BackendDAE.EquationArray eqs,daeEqs;
2457 BackendDAE.Variables vars, daeVars;
2458 BackendDAE.EqSystem subSys;
2459 BackendDAE.AdjacencyMatrixEnhanced me, meT;
2460 BackendDAE.AdjacencyMatrix m;
2461 AvlTreePathFunction.Tree funcs;
2462 list<BackendDAE.Equation> eqLst,eqsInLst;
2463 list<BackendDAE.Var> varLst;
2464 algorithm
2465 //prepare everything
2466 1 size := listLength(varIdcs);
2467
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
1 BackendDAE.EQSYSTEM(orderedVars=daeVars, orderedEqs=daeEqs, matching=BackendDAE.MATCHING(ass1=ass1Sys,ass2=ass2Sys)) := dae;
2468 1 funcs := BackendDAEUtil.getFunctions(shared);
2469 1 eqLst := BackendEquation.getList(eqIdcs,daeEqs);
2470 1 eqs := BackendEquation.listEquation(eqLst);
2471 1 varLst := List.map1r(varIdcs, BackendVariable.getVarAt, daeVars);
2472 1 vars := BackendVariable.listVar1(varLst);
2473 1 subSys := BackendDAEUtil.createEqSystem(vars, eqs);
2474 1 (me,meT,_,_) := BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(subSys,shared,false);
2475 2 (_,m,_,_,_) := BackendDAEUtil.getAdjacencyMatrixScalar(subSys,BackendDAE.SOLVABLE(),SOME(BackendDAEUtil.getFunctions(shared)), BackendDAEUtil.isInitializationDAE(shared));
2476 1 ass1 := arrayCreate(size,-1);
2477 1 ass2 := arrayCreate(size,-1);
2478
2479 // dump system as graphML
2480 1 varAtts := List.threadMap(List.fill(false,listLength(varLst)),List.map(eqIdcs,intString),Util.makeTuple);
2481 1 eqAtts := List.threadMap(List.fill(false,listLength(eqLst)),List.map(varIdcs,intString),Util.makeTuple);
2482 1 BackendDump.dumpBipartiteGraphStrongComponent2(vars,eqs,m,varAtts,eqAtts,"shuffle_pre");
2483
2484 //start reshuffling
2485 1 resEqs := reshuffling_post3_selectShuffleEqs(me,meT);
2486 //print("selected equation pairs: "+stringDelimitList(List.map(resEqs,HpcOmTaskGraph.intLstString)," | ")+"\n");
2487
2488 1 eqsInLst := reshuffling_post4_resolveAndReplace(resEqs,eqLst,varLst,me,meT);
2489
2490 // dump system as graphML
2491 //subSys := BackendDAEUtil.createEqSystem(vars, replEqs);
2492 //(me2,_,_,_) := BackendDAEUtil.getAdjacencyMatrixEnhancedScalar(subSys,shared,false);
2493 //BackendDump.dumpBipartiteGraphStrongComponentSolvable(vars,replEqs,me2,varAtts,eqAtts,"shuffle_post");
2494
2495 // the new eqSystem
2496 1 daeEqs := List.threadFold(eqIdcs,eqsInLst,BackendEquation.setAtIndexFirst,daeEqs);
2497 1 daeOut := BackendDAEUtil.setEqSystEqs(dae, daeEqs);
2498 1 daeOut := BackendDAEUtil.setEqSystMatching(daeOut, BackendDAE.MATCHING(ass1Sys, ass2Sys, {}));
2499
2500 2 (daeOut,_,_,_,_) := BackendDAEUtil.getAdjacencyMatrixScalar(daeOut, BackendDAE.NORMAL(), SOME(funcs), BackendDAEUtil.isInitializationDAE(shared));
2501
2502 outRunMatching := true;
2503 end reshuffling_post2;
2504
2505 protected function reshuffling_post3_selectShuffleEqs
2506 input BackendDAE.AdjacencyMatrixEnhanced me;
2507 input BackendDAE.AdjacencyMatrixEnhanced meT;
2508 output list<list<Integer>> resolveEqs;
2509 algorithm
2510 resolveEqs := matchcontinue meT
2511 local
2512 array<Boolean> bArr;
2513 list<Integer> suitableEqs;
2514 list<list<Integer>> eqPairs;
2515 case _
2516 algorithm
2517 1 bArr := Array.map1(me,chooseEquation,meT);
2518 1 (_,suitableEqs) := List.filter1OnTrueSync(arrayList(bArr),boolEq,true,List.intRange(arrayLength(me)));
2519 //print("suitableEqs: \n"+stringDelimitList(List.map(suitableEqs,intString)," / ")+"\n");
2520 1 eqPairs := List.map2(suitableEqs,getEqPairs,me,meT);
2521 1 eqPairs := List.filterOnTrue(eqPairs,List.hasSeveralElements);
2522 then eqPairs;
2523 else
2524 algorithm
2525 ✗ print("reshuffling_post3_selectShuffleEqs failed!\n");
2526 then {};
2527 end matchcontinue;
2528 end reshuffling_post3_selectShuffleEqs;
2529
2530 protected function reshuffling_post4_resolveAndReplace
2531 input list<list<Integer>> resolveEqLst;
2532 input list<BackendDAE.Equation> unassEqsIn;
2533 input list<BackendDAE.Var> unassVarsIn;
2534 input BackendDAE.AdjacencyMatrixEnhanced me;
2535 input BackendDAE.AdjacencyMatrixEnhanced meT;
2536 output list<BackendDAE.Equation> unassEqsOut;
2537 algorithm
2538 unassEqsOut := matchcontinue resolveEqLst
2539 local
2540 Integer maxNum, replEqIdx;
2541 list<Integer> numOfAdjVars, resolveEqs;
2542 list<list<Integer>> rest;
2543 list<BackendDAE.Equation> unassEqs;
2544 BackendDAE.Equation resolvedEq;
2545 case {}
2546 then unassEqsIn;
2547 case resolveEqs::rest
2548 algorithm
2549 2 resolvedEq := resolveEquations(NONE(),resolveEqs,me,meT,unassEqsIn,unassVarsIn);
2550 //BackendDump.dumpEquationList({resolvedEq},"resolvedEq");
2551
2552 //replace a former equation
2553 2 numOfAdjVars := List.map(List.map1(resolveEqs,Array.getIndexFirst,me),listLength);
2554 2 maxNum := List.fold(numOfAdjVars,intMax,listHead(numOfAdjVars));
2555 2 replEqIdx := listGet(resolveEqs,List.position(maxNum,numOfAdjVars));
2556 //BackendDump.dumpEquationList(unassEqsIn," not updated unassEqs");
2557 2 unassEqs := List.replaceAt(resolvedEq,replEqIdx,unassEqsIn);
2558 //print("replace equation "+intString(replEqIdx)+"\n");
2559 //BackendDump.dumpEquationList(unassEqs,"updated unassEqs");
2560 2 then reshuffling_post4_resolveAndReplace(rest,unassEqs,unassVarsIn,me,meT);
2561 else
2562 algorithm
2563 ✗ print("reshuffling_post4_resolveAndReplace failed!\n");
2564 ✗ then fail();
2565 end matchcontinue;
2566 end reshuffling_post4_resolveAndReplace;
2567
2568 protected function getEqPairs
2569 input Integer eq;
2570 input BackendDAE.AdjacencyMatrixEnhanced me;
2571 input BackendDAE.AdjacencyMatrixEnhanced meT;
2572 output list<Integer> lstOut;
2573 protected
2574 list<Integer> vars,eqs;
2575 algorithm
2576 2 vars := List.map(arrayGet(me,eq),Util.tuple31);
2577 //print("vars: \n"+stringDelimitList(List.map(vars,intString)," / ")+"\n");
2578 2 eqs := List.map(List.flatten(List.map1(vars,Array.getIndexFirst,meT)),Util.tuple31);
2579 //print("eqs: \n"+stringDelimitList(List.map(eqs,intString)," / ")+"\n");
2580 2 eqs := getDoublicates(eqs);
2581 //print("eqs: \n"+stringDelimitList(List.map(eqs,intString)," / ")+"\n");
2582 2 lstOut := List.consOnTrue(not listMember(eq, eqs), eq, eqs);
2583 end getEqPairs;
2584
2585 protected function chooseEquation
2586 input list<BackendDAE.AdjacencyMatrixElementEnhancedEntry> row;
2587 input BackendDAE.AdjacencyMatrixEnhanced meT;
2588 output Boolean chooseThis;
2589 protected
2590 Boolean b1,b2,b3;
2591 list<Integer> vars,eqs,numEqs;
2592 list<list<Integer>> eqLst;
2593 algorithm
2594 35 vars := List.map(row,Util.tuple31);
2595 35 b1 := intEq(listLength(row),2); // only two variables
2596 35 eqLst := List.mapList((List.map1(vars,Array.getIndexFirst,meT)),Util.tuple31);
2597 35 numEqs := List.map(eqLst,listLength);
2598 35 b3 := List.applyAndFold1(numEqs,boolOr,intEq,2,false); // at least one adjacent variable hast only 2 adj equations
2599 35 eqs := List.flatten(eqLst);
2600 35 b2 := intEq(listLength(eqs),listLength(List.unique(eqs))+2);
2601
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 33 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
35 b1 := b1 and b2 and b3;
2602
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
2 chooseThis := b1 and List.applyAndFold(row,boolAnd,isSolvable,true);
2603 end chooseEquation;
2604
2605 protected function getDoublicates
2606 input list<Integer> lstIn; // only positive Integer
2607 output list<Integer> lstOut;
2608 protected
2609 Integer max;
2610 array<Integer> arr;
2611 algorithm
2612 2 max := List.fold(lstIn,intMax,listHead(lstIn));
2613 2 arr := arrayCreate(max,-1);
2614 2 List.map1_0(lstIn,getDoublicates2,arr);
2615 2 (_,lstOut) := List.filter1OnTrueSync(arrayList(arr),intGe,1,List.intRange(arrayLength(arr)));
2616 end getDoublicates;
2617
2618 protected function getDoublicates2
2619 input Integer idx;
2620 input array<Integer> arr;
2621 protected
2622 Integer entry;
2623 algorithm
2624 22 entry := arrayGet(arr,idx);
2625 22 arrayUpdate(arr,idx,entry+1);
2626 end getDoublicates2;
2627
2628 protected function isSolvable
2629 input BackendDAE.AdjacencyMatrixElementEnhancedEntry entry;
2630 output Boolean solvable;
2631 algorithm
2632 4 solvable := not Tearing.unsolvable({entry});
2633 end isSolvable;
2634
2635 public function resolveEquations
2636 input Option<BackendDAE.Equation> eq;
2637 input list<Integer> loopIn;
2638 input BackendDAE.AdjacencyMatrixEnhanced me;
2639 input BackendDAE.AdjacencyMatrixEnhanced meT;
2640 input list<BackendDAE.Equation> eqsIn;
2641 input list<BackendDAE.Var> varsIn;
2642 output BackendDAE.Equation eqOut;
2643 protected
2644 algorithm
2645 eqOut := matchcontinue(eq, loopIn)
2646 local
2647 Integer startEq,nextEq,sharedVar;
2648 list<Integer> rest,vars1,vars2,numEqs;
2649 BackendDAE.Equation eq1,eq2;
2650 BackendDAE.Var var;
2651 BackendDAE.EquationAttributes attr;
2652 DAE.Exp lhs1, lhs2, rhs1, rhs2 ,varExp,eqExp;
2653 DAE.ElementSource source;
2654 case(SOME(eq1), {})
2655 algorithm
2656 // resolved the whole cycle
2657 then eq1;
2658 case(NONE(), startEq::rest)
2659 algorithm
2660 // start resolving the first 2 equations
2661
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 nextEq::rest := rest;
2662 2 vars1 := List.map(arrayGet(me,startEq),Util.tuple31);
2663 2 vars2 := List.map(arrayGet(me,nextEq),Util.tuple31);
2664 2 vars1 := List.intersectionOnTrue(vars1,vars2,intEq);
2665 2 numEqs := List.map(List.map1(vars1,Array.getIndexFirst,meT),listLength);
2666 2 (_,vars1) := List.filter1OnTrueSync(numEqs,intEq,2,vars1);
2667 2 sharedVar := listHead(vars1);
2668 2 eq1 := listGet(eqsIn,startEq);
2669 2 eq2 := listGet(eqsIn,nextEq);
2670 2 var := listGet(varsIn,sharedVar);
2671 2 varExp := Expression.crefExp(BackendVariable.varCref(var));
2672 //BackendDump.dumpEquationList({eq1},"eq1");
2673 //BackendDump.dumpEquationList({eq2},"eq2");
2674 //BackendDump.dumpVarList({var},"var");
2675
2676
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 BackendDAE.EQUATION(exp=lhs1,scalar=rhs1,source=source,attr=attr) := eq1;
2677
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 BackendDAE.EQUATION(exp=lhs2,scalar=rhs2) := eq2;
2678 2 (eqExp,_) := ExpressionSolve.solve(lhs1,rhs1,varExp);
2679 //BackendDump.dumpEquationList({eq1},"solved Eq");
2680
2681 2 (lhs2,_) := Expression.replaceExp(lhs2,varExp,eqExp);
2682 2 (rhs2,_) := Expression.replaceExp(rhs2,varExp,eqExp);
2683 2 (lhs2,_) := ExpressionSimplify.simplify(lhs2);
2684 2 (rhs2,_) := ExpressionSimplify.simplify(rhs2);
2685 2 eq2 := BackendDAE.EQUATION(lhs2,rhs2,source,attr);
2686 //BackendDump.dumpEquationList({eq2},"resolved Eq");
2687 2 then resolveEquations(SOME(eq2),rest,me,meT,eqsIn,varsIn);
2688 else
2689 algorithm
2690 ✗ print("resolveEquations failed!\n");
2691 ✗ then fail();
2692 end matchcontinue;
2693 end resolveEquations;
2694
2695 // =============================================================================
2696 // section for postOptModule >solveLinearSystem<<
2697 //
2698 // solve linear system of equations (A x = b)
2699 // =============================================================================
2700
2701 public function solveLinearSystem
2702 input BackendDAE.BackendDAE inDAE;
2703 output BackendDAE.BackendDAE outDAE;
2704 protected
2705 Integer maxSize = Flags.getConfigInt(Flags.MAX_SIZE_FOR_SOLVE_LINIEAR_SYSTEM);
2706 Boolean b = 1 < maxSize;
2707 algorithm
2708
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if b then
2709 2 (outDAE,_) := BackendDAEUtil.mapEqSystemAndFold(inDAE, solveLinearSystem0, (false,1,maxSize));
2710 else
2711 outDAE := inDAE;
2712 end if;
2713 end solveLinearSystem;
2714
2715 protected function solveLinearSystem0
2716 input BackendDAE.EqSystem isyst;
2717 input BackendDAE.Shared inShared;
2718 input tuple<Boolean,Integer,Integer> inTpl;
2719 output BackendDAE.EqSystem osyst;
2720 output BackendDAE.Shared outShared;
2721 output tuple<Boolean,Integer,Integer> outTpl;
2722 protected
2723 BackendDAE.StrongComponents comps;
2724 algorithm
2725
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(comps=comps)) := isyst;
2726 3 (osyst, outShared, outTpl) := solveLinearSystem1(isyst, inShared, comps, inTpl);
2727 end solveLinearSystem0;
2728
2729 protected function solveLinearSystem1
2730 input BackendDAE.EqSystem isyst;
2731 input BackendDAE.Shared ishared;
2732 input BackendDAE.StrongComponents inComps;
2733 input tuple<Boolean,Integer,Integer> inTpl;
2734 output BackendDAE.EqSystem osyst = isyst;
2735 output BackendDAE.Shared oshared = ishared;
2736 output tuple<Boolean,Integer,Integer> outTpl;
2737 protected
2738 Boolean b;
2739 Boolean runMatching;
2740 list<Integer> ii = {};
2741 Integer offset, maxSize;
2742 algorithm
2743 3 (runMatching, offset, maxSize) := inTpl;
2744
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 3 times.
42 for comp in inComps loop
2745 39 (osyst,oshared,b,ii,offset) := solveLinearSystem2(osyst,oshared,comp,ii, offset, maxSize);
2746
4/4
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 30 times.
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 2 times.
39 runMatching := runMatching or b;
2747 end for;
2748
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 outTpl := (runMatching, offset,maxSize);
2749
2750
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
3 if runMatching then
2751 osyst := match osyst
2752 local
2753 BackendDAE.Variables vars;
2754 BackendDAE.EquationArray eqns;
2755 BackendDAE.EqSystem syst;
2756 case syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns)
2757 algorithm
2758 // remove empty entries from vars/eqns
2759 3 eqns := List.fold(ii,BackendEquation.delete,eqns);
2760 3 syst.orderedVars := BackendVariable.listVar1(BackendVariable.varList(vars));
2761 3 syst.orderedEqs := BackendEquation.listEquation(BackendEquation.equationList(eqns));
2762 3 then
2763 BackendDAEUtil.clearEqSyst(syst);
2764 end match;
2765 end if;
2766 end solveLinearSystem1;
2767
2768 protected function solveLinearSystem2
2769 input BackendDAE.EqSystem isyst;
2770 input BackendDAE.Shared ishared;
2771 input BackendDAE.StrongComponent comp;
2772 input list<Integer> ii;
2773 input Integer offset;
2774 input Integer maxSize;
2775 output BackendDAE.EqSystem osyst;
2776 output BackendDAE.Shared oshared;
2777 output Boolean outRunMatching;
2778 output list<Integer> oi;
2779 output Integer offset_;
2780 algorithm
2781 (osyst,oshared,outRunMatching, oi, offset_):=
2782 matchcontinue (isyst,ishared,comp)
2783 local
2784 BackendDAE.Variables vars;
2785 BackendDAE.EquationArray eqns;
2786 list<BackendDAE.Equation> eqn_lst;
2787 list<BackendDAE.Var> var_lst;
2788 list<Integer> eindex,vindx;
2789 list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
2790 BackendDAE.EqSystem syst;
2791 BackendDAE.Shared shared;
2792 Integer toffset;
2793
2794 case ( syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns), shared,
2795 (BackendDAE.EQUATIONSYSTEM( eqns=eindex, vars=vindx, jac=BackendDAE.FULL_JACOBIAN(SOME(jac)), jacType=BackendDAE.JAC_LINEAR()))
2796 )
2797 algorithm
2798 2 eqn_lst := BackendEquation.getList(eindex,eqns);
2799 2 var_lst := List.map1r(vindx, BackendVariable.getVarAt, vars);
2800
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
2 true := listLength(var_lst) <= maxSize;
2801
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
2 ({},_) := List.splitOnTrue(var_lst, BackendVariable.isStateVar) "TODO: fix BackendDAEUtil.getEqnSysRhs for x and der(x)";
2802 2 (syst,shared, toffset) := solveLinearSystem3(syst,shared,eqn_lst,eindex,var_lst,vindx,jac,offset);
2803 2 then (syst,shared,true, listAppend(eindex, ii), toffset);
2804 else (isyst,ishared,false, ii, offset);
2805 end matchcontinue;
2806 end solveLinearSystem2;
2807
2808 protected function solveLinearSystem3
2809 input BackendDAE.EqSystem inSyst;
2810 input BackendDAE.Shared ishared;
2811 input list<BackendDAE.Equation> eqn_lst;
2812 input list<Integer> eqn_indxs;
2813 input list<BackendDAE.Var> var_lst;
2814 input list<Integer> var_indxs;
2815 input list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
2816 input Integer offset;
2817 output BackendDAE.EqSystem osyst;
2818 output BackendDAE.Shared oshared;
2819 output Integer offset_;
2820 algorithm
2821 (osyst,oshared, offset_):=
2822 match (inSyst, ishared)
2823 local
2824 BackendDAE.Variables vars;
2825 BackendDAE.EquationArray eqns;
2826 list<DAE.Exp> beqs;
2827 list<DAE.ComponentRef> names;
2828 AvlTreePathFunction.Tree funcs;
2829 BackendDAE.Shared shared;
2830 BackendDAE.EqSystem syst;
2831 Integer n;
2832
2833 case ( syst as BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns),
2834 shared as BackendDAE.SHARED(functionTree=funcs) )
2835 algorithm
2836 2 (beqs, _) := BackendDAEUtil.getEqnSysRhs( BackendEquation.listEquation(eqn_lst),
2837 BackendVariable.listVar1(var_lst), SOME(funcs) );
2838 2 beqs := listReverse(beqs);
2839 2 n := listLength(beqs);
2840 2 names := List.map(var_lst, BackendVariable.varCref);
2841 2 (eqns, vars, n, shared) := solveLinearSystem4(beqs, jac, names, var_lst, n, eqns, vars, offset, shared);
2842 4 syst.orderedVars := vars; syst.orderedEqs := eqns;
2843 2 syst := BackendDAEUtil.setEqSystMatrices(syst);
2844 //eqns = List.fold(eqn_indxs,BackendEquation.delete,eqns);
2845
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 then
2846 (syst, shared, n);
2847 end match;
2848 end solveLinearSystem3;
2849
2850 protected function solveLinearSystem4
2851 "
2852 author: Vitalij Ruge
2853 "
2854 input list<DAE.Exp> b_lst;
2855 input list<tuple<Integer, Integer, BackendDAE.Equation>> jac;
2856 input list<DAE.ComponentRef> cr_x;
2857 input list<BackendDAE.Var> var_lst;
2858 input Integer n;
2859 input BackendDAE.EquationArray ieqns;
2860 input BackendDAE.Variables ivars;
2861 input Integer offset;
2862 input BackendDAE.Shared ishared;
2863 output BackendDAE.EquationArray oeqns = ieqns;
2864 output BackendDAE.Variables ovars = ivars;
2865 output Integer offset_ = offset + 1;
2866 output BackendDAE.Shared oshared = ishared;
2867 protected
2868 array<DAE.Exp> R;
2869 array<DAE.Exp> Qb = arrayCreate(n,DAE.RCONST(0.0));
2870 array<DAE.Exp> b = arrayCreate(n,DAE.RCONST(0.0));
2871 array<DAE.Exp> A = arrayCreate(n*n,DAE.RCONST(0.0));
2872 array<DAE.Exp> ax = arrayCreate(n,DAE.RCONST(0.0));
2873 array<DAE.Exp> scaled_x = arrayCreate(n,DAE.RCONST(0.0));
2874 array<DAE.Exp> scaleA = arrayCreate(n,DAE.RCONST(0.0));
2875
2876 DAE.Exp a;
2877 Integer m, ii, jj, mm;
2878 list<DAE.Exp> x_lst = List.map(cr_x, Expression.crefExp);
2879 list<BackendDAE.Var> vars = var_lst;
2880 BackendDAE.Var var;
2881 BackendDAE.Equation eqn;
2882 list<tuple<Integer, Integer, BackendDAE.Equation>> jac_ = jac;
2883 algorithm
2884 2 mm := listLength(jac);
2885 //A
2886
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
406 for i in 1:mm loop
2887
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 404 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 404 times.
404 (jj,ii,BackendDAE.RESIDUAL_EQUATION(exp = a)) :: jac_ := jac_; // jac(1) = a11, jac(2)=a12,.., jac(n+1) = an1
2888 404 m := ii + (jj-1)*n;
2889 404 (a, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$A$" + intString(m), offset, oeqns, ovars, oshared);
2890 404 arrayUpdate(A,m,a);
2891 end for;
2892
2893 // calc scale A
2894
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
24 for i in 1:n loop
2895 22 m := (i-1)*n;
2896
2/2
✓ Branch 0 taken 404 times.
✓ Branch 1 taken 22 times.
426 a := Expression.makeSum1(list(Expression.makeAbs(arrayGet(A, m+j)) for j in 1:n));
2897 22 arrayUpdate(scaleA,i,a);
2898 end for;
2899
2900 // update A scaling
2901
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
24 for i in 1:n loop
2902 22 m := (i-1)*n;
2903 404 for j in 1:n loop
2904 404 a := arrayGet(A,j+m);
2905
1/2
✓ Branch 1 taken 404 times.
✗ Branch 2 not taken.
404 if not Expression.isZero(a) then
2906 404 a := Expression.expDiv(a, arrayGet(scaleA,j));
2907 404 (a, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$sA$" + intString(i + (j-1)*n), offset, oeqns, ovars, oshared);
2908 404 arrayUpdate(A, j+m, a);
2909 end if;
2910 end for;
2911 end for;
2912
2913 // b
2914 m := 1;
2915
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 2 times.
24 for b_ in b_lst loop
2916 22 (a, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(b_, "QR$b$" + intString(m), offset, oeqns, ovars, oshared);
2917 22 arrayUpdate(b, m, a);
2918 22 m := m + 1;
2919 end for;
2920 //qrDecomposition3(b, n, false, "b");
2921
2922 // x
2923 m := 1;
2924
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 2 times.
24 for xx in x_lst loop
2925
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
22 var :: vars := vars;
2926
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 22 times.
22 if BackendVariable.isStateVar(var) then
2927 ✗ arrayUpdate(ax,m,Expression.expDer(xx));
2928 else
2929 22 arrayUpdate(ax,m,xx);
2930 end if;
2931 22 m := m + 1;
2932 end for;
2933
2934 //rescale x
2935 // scale_x = x/factor
2936
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
24 for i in 1:n loop
2937 22 a := Expression.expMul(arrayGet(ax,i), arrayGet(scaleA,i));
2938 22 (a, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$sx$" + intString(i), offset, oeqns, ovars, oshared);
2939 22 arrayUpdate(scaled_x,i,a);
2940 end for;
2941
2942 //qrDecomposition3(ax, n, false, "x");
2943
2944 // A*x = b -> R*x = Q'b
2945 //(R, Qb, oeqns, ovars, oshared) := qrDecomposition(A, n, b, oeqns, ovars, offset, oshared);
2946 2 (R, Qb, oeqns, ovars, oshared) := qrDecompositionHouseholder(A, n, b, oeqns, ovars, offset, oshared);
2947
2948 // R*x = Q'*b where x is scaled
2949
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 for i in n:-1:1 loop
2950 22 m := (i-1)*n;
2951
2/2
✓ Branch 0 taken 213 times.
✓ Branch 1 taken 22 times.
235 a := Expression.makeSum1(list(Expression.expMul(arrayGet(R, m + j), arrayGet(scaled_x, j)) for j in i:n));
2952 22 eqn := BackendDAE.EQUATION(a, arrayGet(Qb,i), DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN);
2953 22 eqn := BackendEquation.solveEquation(eqn, arrayGet(scaled_x,i), NONE());
2954 22 oeqns := BackendEquation.add(eqn, oeqns);
2955 end for;
2956
2957 end solveLinearSystem4;
2958
2959 protected function qrDecompositionHouseholder
2960 "
2961 QR-Decomposition based on Householder
2962 author: Vitalij Ruge
2963 "
2964 input array<DAE.Exp> A;
2965 input Integer n;
2966 input array<DAE.Exp> ib;
2967 input BackendDAE.EquationArray ieqns;
2968 input BackendDAE.Variables ivars;
2969 input Integer offset;
2970 input BackendDAE.Shared ishared;
2971 output array<DAE.Exp> R = A;
2972 output array<DAE.Exp> b = ib;
2973 output BackendDAE.EquationArray oeqns = ieqns;
2974 output BackendDAE.Variables ovars = ivars;
2975 output BackendDAE.Shared oshared = ishared;
2976
2977 protected
2978 array<DAE.Exp> cA = arrayCreate(n,DAE.RCONST(0.0)) "column of A";
2979 array<DAE.Exp> v = arrayCreate(n,DAE.RCONST(0.0)) "vec in dyadic tensor";
2980 DAE.Exp alpha;
2981 DAE.Exp y1;
2982 DAE.Exp h,h2;
2983 DAE.Exp e1,e2,e;
2984 Integer m "cuurrent size";
2985 Integer nn = n-1;
2986 Integer idxVars = 1 "index for tmp vars";
2987 Integer shift;
2988 algorithm
2989 //qrDecomposition3(A,n,true,"A");
2990
2991
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
22 for iter in 1:nn loop
2992 20 m := n - iter + 1;
2993
2994 //first column
2995 20 qrGet_cA(A,iter,1,n,v);
2996 y1 := arrayGet(v,1);
2997
2998 20 alpha := qrCalc_alpha(v,y1,m);
2999 20 (alpha, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(alpha, "QR$a$" + intString(iter), offset, oeqns, ovars, oshared);
3000
3001 // calc v
3002 20 e := Expression.expAdd(y1,alpha);
3003 20 (e, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(e, "QR$y1$" + intString(iter), offset, oeqns, ovars, oshared);
3004 20 arrayUpdate(v,1, e);
3005
3006 // helper for const factor
3007 20 h := Expression.expAdd(y1,alpha);
3008 20 h := Expression.expMul(alpha,h);
3009 20 h := Expression.negate(h);
3010 20 (h, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(h, "QR$h$" + intString(iter), offset, oeqns, ovars, oshared);
3011
3012 20 shift := (iter-1)*n + iter;
3013 //update R
3014 20 arrayUpdate(R, shift, Expression.negate(alpha));
3015
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
211 for j in 2:m loop
3016 191 arrayUpdate(R, shift + (j-1)*n, DAE.RCONST(0.0));
3017 end for;
3018
3019
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
211 for col in 2:m loop
3020 191 qrGet_cA(A,iter,col,n,cA);
3021 191 h2 := Expression.makeScalarProduct(v, cA);
3022 191 h2 := Expression.expDiv(h2,h);
3023 191 (h2, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(h2, "QR$h2$" + intString(idxVars), offset, oeqns, ovars, oshared);
3024 191 idxVars := idxVars + 1;
3025 //vec add
3026 2662 for j in 1:m loop
3027 2662 e1 := arrayGet(cA,j);
3028 2662 e2 := arrayGet(v,j);
3029 2662 e := Expression.expAdd(e1,Expression.expMul(h2, e2));
3030 2662 (e, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(e, "QR$R$" + intString(idxVars), offset, oeqns, ovars, oshared);
3031 2662 idxVars := idxVars + 1;
3032 //update A
3033 2662 arrayUpdate(A, shift + (j-1)*n + col-1, e);
3034 end for;
3035 end for;
3036 //update b
3037
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
231 for j in 1:m loop
3038 211 arrayUpdate(cA,j, arrayGet(b,iter-1 + j));
3039 end for;
3040
3041 20 h2 := Expression.makeScalarProduct(v, cA);
3042 20 h2 := Expression.expDiv(h2, h);
3043 20 (h2, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(h2, "QR$b_$" + intString(idxVars), offset, oeqns, ovars, oshared);
3044 20 idxVars := idxVars + 1;
3045
3046 //vec add
3047
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
231 for j in 1:m loop
3048 211 e1 := arrayGet(cA, j);
3049 211 e2 := arrayGet(v,j);
3050 211 e := Expression.expAdd(e1,Expression.expMul(h2, e2));
3051 211 e := Expression.expand(e);
3052 211 e := ExpressionSimplify.simplify2(e);
3053 211 (e, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(e, "QR$b$" + intString(idxVars), offset, oeqns, ovars, oshared);
3054 211 idxVars := idxVars + 1;
3055 //update b
3056 211 arrayUpdate(b, iter-1+j, e);
3057 end for;
3058
3059 end for;
3060 //qrDecomposition3(A,n,true,"R");
3061
3062 end qrDecompositionHouseholder;
3063
3064 protected function qrGet_cA
3065 "
3066 helper for QR-Decomposition based on Householder
3067 return column j in A for iteration iter
3068 author: Vitalij Ruge
3069 "
3070 input array<DAE.Exp> A;
3071 input Integer iter "iteration";
3072 input Integer j "column";
3073 input Integer n "size";
3074 input array<DAE.Exp> cA "output";
3075 protected
3076 Integer shift = (iter-1)*n + iter + j - 1;
3077 Integer m = n-iter+1;
3078 algorithm
3079
3080
1/2
✓ Branch 0 taken 211 times.
✗ Branch 1 not taken.
3084 for i in 1:m loop
3081 2873 arrayUpdate(cA, i, arrayGet(A,shift+(i-1)*n));
3082 end for;
3083
3084
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 189 times.
1522 for i in m+1:n loop
3085 1311 arrayUpdate(cA, i, DAE.RCONST(0.0));
3086 end for;
3087
3088 end qrGet_cA;
3089
3090 protected function qrCalc_alpha
3091 "
3092 helper for QR-Decomposition based on Householder
3093 calculate -> min loss of significance
3094 author: Vitalij Ruge
3095 "
3096 input array<DAE.Exp> y;
3097 input DAE.Exp y1;
3098 input Integer m;
3099 output DAE.Exp alpha;
3100 protected
3101 DAE.Exp sgn_y1 = Expression.makeSign(y1);
3102 DAE.Exp norm_y = Expression.lenVec(y);
3103 algorithm
3104
2/2
✓ Branch 0 taken 211 times.
✓ Branch 1 taken 20 times.
231 norm_y := Expression.makeSum1(list( Expression.expPow(arrayGet(y,j), DAE.RCONST(2.0)) for j in 1:m ));
3105 20 norm_y := Expression.makePureBuiltinCall("sqrt",{norm_y},DAE.T_REAL_DEFAULT);
3106 20 alpha := Expression.expMul(sgn_y1,norm_y);
3107 end qrCalc_alpha;
3108
3109
3110 protected function qrDecomposition
3111 "
3112 author: Vitalij Ruge
3113 "
3114 input array<DAE.Exp> A;
3115 input Integer n;
3116 input array<DAE.Exp> ib;
3117 input BackendDAE.EquationArray ieqns;
3118 input BackendDAE.Variables ivars;
3119 input Integer offset;
3120 input BackendDAE.Shared ishared;
3121 output array<DAE.Exp> R = arrayCreate(n*n,DAE.RCONST(0.0));
3122 output array<DAE.Exp> b = arrayCreate(n,DAE.RCONST(0.0));
3123 output BackendDAE.EquationArray oeqns = ieqns;
3124 output BackendDAE.Variables ovars = ivars;
3125 output BackendDAE.Shared oshared;
3126 protected
3127 array<DAE.Exp> Q = arrayCreate(n*n,DAE.RCONST(0.0));
3128 array<DAE.Exp> v = arrayCreate(n,DAE.RCONST(0.0));
3129 array<DAE.Exp> u = arrayCreate(n,DAE.RCONST(0.0));
3130 array<DAE.Exp> x,y;
3131 DAE.Exp a;
3132 Integer kk = 1;
3133 Integer m = n-1;
3134 Integer nn;
3135 algorithm
3136 //Gram–Schmidt process
3137 ✗ v := qrDecomposition1(A,n,kk);
3138 ✗ (u, oeqns, ovars, oshared) := BackendEquation.normalizationVec(v,"QR$NOM$" + intString(kk), offset, oeqns, ovars, ishared);
3139
3140 ✗ for j in 1:n loop
3141 ✗ (a,_) := ExpressionSimplify.simplify(arrayGet(u,j));
3142 ✗ (a, oeqns, ovars,oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$Q$" + intString(kk + (j-1)*n), offset, oeqns, ovars,oshared);
3143 ✗ arrayUpdate(Q, kk + (j-1)*n, a);
3144 end for;
3145
3146 ✗ for k in 1:m loop
3147 ✗ v := qrDecomposition1(A,n,k+1);
3148 ✗ for j in 1:k loop
3149 ✗ u := qrDecomposition1(Q,n,j);
3150 ✗ (v, oeqns, ovars,oshared) := gramSchmidtProcessHelper(v,u,"QR$W$" + intString(kk) + "$" + intString(kk), offset, oeqns, ovars,oshared);
3151 ✗ kk := kk +1;
3152 end for;
3153 ✗ (u, oeqns, ovars, oshared) := BackendEquation.normalizationVec(v,"QR$NOM$" + intString(k+1), offset, oeqns, ovars, oshared);
3154 //qrDecomposition3(u, n, false, "u");
3155 ✗ for j in 1:n loop
3156 ✗ nn := k+1 + (j-1)*n;
3157 ✗ (a,_) := ExpressionSimplify.simplify(arrayGet(u,j));
3158 ✗ (a, oeqns, ovars, oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$Q$" + intString(nn), offset, oeqns, ovars, oshared);
3159 ✗ arrayUpdate(Q, nn, a);
3160 end for;
3161
3162 end for;
3163 //qrDecomposition3(Q, n, true, "Q");
3164
3165
3166 // R
3167 ✗ for i in 1:n loop
3168 ✗ x := qrDecomposition1(Q,n,i);
3169 ✗ m := (i-1)*n;
3170 //qrDecomposition3(x, n, false, "x" + intString(i));
3171 ✗ for j in i:n loop
3172 ✗ y := qrDecomposition1(A,n,j);
3173 ✗ a := Expression.makeScalarProduct(x,y);
3174 ✗ (a, oeqns, ovars,oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$R$" + intString(m + j), offset, oeqns, ovars,oshared);
3175 ✗ arrayUpdate(R, m+j, a);
3176 end for;
3177 end for;
3178
3179 //qrDecomposition3(R, n, true, "R");
3180 //qrDecomposition3(Q, n, true, "Q");
3181 // Q*b
3182 ✗ for i in 1:n loop
3183 ✗ x := qrDecomposition1(Q,n,i);
3184 //qrDecomposition3(x, n, false, "x" + intString(i));
3185 ✗ a := Expression.makeScalarProduct(x,ib);
3186 ✗ (a, oeqns, ovars,oshared) := BackendEquation.makeTmpEqnForExp(a, "QR$Qb$" + intString(i), offset, oeqns, ovars, oshared);
3187 ✗ arrayUpdate(b, i, a);
3188 end for;
3189 //qrDecomposition3(b, n, false, "Qb");
3190
3191 end qrDecomposition;
3192
3193 protected function qrDecomposition1
3194 "return column of A"
3195 input array<DAE.Exp> A;
3196 input Integer sizeA;
3197 input Integer i;
3198 output array<DAE.Exp> column = arrayCreate(sizeA,DAE.RCONST(0.0)) "A(:,i)";
3199 algorithm
3200 ✗ for j in 1:sizeA loop
3201 ✗ arrayUpdate(column, j, arrayGet(A,i + (j-1)*sizeA));
3202 end for;
3203 end qrDecomposition1;
3204
3205 protected function qrDecomposition2
3206 "return row of A"
3207 input array<DAE.Exp> A;
3208 input Integer sizeA;
3209 input Integer i "row";
3210 output array<DAE.Exp> row = arrayCreate(sizeA,DAE.RCONST(0.0)) "A(:,i)";
3211 protected
3212 Integer k = i - 1;
3213 algorithm
3214 ✗ for j in 1:sizeA loop
3215 ✗ arrayUpdate(row, j, arrayGet(A,j + k*sizeA));
3216 end for;
3217 end qrDecomposition2;
3218
3219 protected function qrDecomposition3
3220 "for debug"
3221 input array<DAE.Exp> A;
3222 input Integer sizeA;
3223 input Boolean isMat;
3224 input String s;
3225 protected
3226 Integer n = sizeA;
3227 Integer m = if isMat then sizeA else 1;
3228 algorithm
3229 ✗ print("\n");
3230 ✗ for i in 1:n loop
3231 ✗ print("\n");
3232 ✗ for j in 1:m loop
3233 ✗ print(s + "(" + intString(i) + "," + intString(j) + ") = " + ExpressionBasics.printExpStr(arrayGet(A, (i-1)*m + j)) + "\t");
3234 end for;
3235 end for;
3236 ✗ print("\n");
3237 end qrDecomposition3;
3238
3239 protected function gramSchmidtProcessHelper
3240 "
3241 author: Vitalij Ruge
3242 small step inside gram–schmidt process
3243 "
3244 input array<DAE.Exp> w;
3245 input array<DAE.Exp> u;
3246 input String name "var name";
3247 input Integer offset;
3248 input BackendDAE.EquationArray ieqns;
3249 input BackendDAE.Variables ivars;
3250 input BackendDAE.Shared ishared;
3251 output array<DAE.Exp> v;
3252 output BackendDAE.EquationArray oeqns;
3253 output BackendDAE.Variables ovars;
3254 output BackendDAE.Shared oshared;
3255 protected
3256 DAE.Exp h = Expression.makeScalarProduct(w,u);
3257 Integer n = arrayLength(w);
3258 algorithm
3259 ✗ (h,oeqns,ovars,oshared) := BackendEquation.makeTmpEqnForExp(h, name + "_h", offset, ieqns, ivars, ishared);
3260 ✗ v := Array.map1(u, Expression.expMul, h);
3261 ✗ v := Expression.subVec(w,v);
3262 ✗ for i in 1:n loop
3263 ✗ (h,oeqns,ovars,oshared) := BackendEquation.makeTmpEqnForExp(arrayGet(v,i), name + "_" + intString(i), offset, oeqns, ovars, oshared);
3264 ✗ arrayUpdate(v,i,h);
3265 end for;
3266
3267 end gramSchmidtProcessHelper;
3268
3269 annotation(__OpenModelica_Interface="backend");
3270 end ResolveLoops;
3271
3272