diff --git a/Compiler/BackEnd/BackendDAEUtil.mo b/Compiler/BackEnd/BackendDAEUtil.mo index 37d17e47e84..4293d5551f4 100644 --- a/Compiler/BackEnd/BackendDAEUtil.mo +++ b/Compiler/BackEnd/BackendDAEUtil.mo @@ -6594,7 +6594,7 @@ algorithm end if; // pre-optimization phase - (optdae, Util.SUCCESS()) := preOptimizeDAE(inDAE, preOptModules); + optdae := preOptimizeDAE(inDAE, preOptModules); // transformation phase (matching and sorting using index reduction method) sode := causalizeDAE(optdae, NONE(), matchingAlgorithm, daeHandler, true); @@ -6609,7 +6609,7 @@ algorithm end if; // post-optimization phase - (optsode, Util.SUCCESS()) := postOptimizeDAE(sode, postOptModules, matchingAlgorithm, daeHandler); + optsode := postOptimizeDAE(sode, postOptModules, matchingAlgorithm, daeHandler); sode1 := FindZeroCrossings.findZeroCrossings(optsode); SimCodeFunctionUtil.execStat("findZeroCrossings"); @@ -6648,52 +6648,46 @@ protected list> preOptModules; algorithm preOptModules := getPreOptModules(strPreOptModules); - (outDAE,Util.SUCCESS()) := preOptimizeDAE(inDAE,preOptModules); + outDAE := preOptimizeDAE(inDAE, preOptModules); end preOptimizeBackendDAE; protected function preOptimizeDAE " This function runs the pre-optimization modules." input BackendDAE.BackendDAE inDAE; - input list> optModules; - output BackendDAE.BackendDAE outDAE; - output Util.Status status; + input list> inPreOptModules; + output BackendDAE.BackendDAE outDAE = inDAE; +protected + BackendDAEFunc.preOptimizationDAEModule optModule; + String moduleStr; + Boolean stopOnFailure; + BackendDAE.EqSystems systs; + BackendDAE.Shared shared; algorithm - (outDAE, status) := matchcontinue (inDAE, optModules) - local - BackendDAE.BackendDAE dae, dae1; - BackendDAEFunc.preOptimizationDAEModule optModule; - list> rest; - String str, moduleStr; - Boolean b; - BackendDAE.EqSystems systs; - BackendDAE.Shared shared; - - case (_, {}) - equation - if Flags.isSet(Flags.OPT_DAE_DUMP) then - print("pre-optimization done.\n"); - end if; - then (inDAE,Util.SUCCESS()); - - case (_, (optModule, moduleStr, _)::rest) equation - BackendDAE.DAE(systs, shared) = optModule(inDAE); - (systs, shared) = filterEmptySystems(systs, shared); - dae = BackendDAE.DAE(systs, shared); + for preOptModule in inPreOptModules loop + (optModule, moduleStr, stopOnFailure) := preOptModule; + try + BackendDAE.DAE(systs, shared) := optModule(outDAE); + (systs, shared) := filterEmptySystems(systs, shared); + outDAE := BackendDAE.DAE(systs, shared); SimCodeFunctionUtil.execStat("preOpt " + moduleStr); if Flags.isSet(Flags.OPT_DAE_DUMP) then print(stringAppendList({"\npre-optimization module ", moduleStr, ":\n\n"})); - BackendDump.printBackendDAE(dae); + BackendDump.printBackendDAE(outDAE); end if; - (dae1,status) = preOptimizeDAE(dae, rest); - then (dae1, status); - - case (_, (_, moduleStr, b)::rest) equation + else SimCodeFunctionUtil.execStat(" preOpt " + moduleStr); - str = stringAppendList({"pre-optimization module ", moduleStr, " failed."}); - Error.addMessage(Error.INTERNAL_ERROR, {str}); - (dae,status) = preOptimizeDAE(inDAE,rest); - then (dae, if b then Util.FAILURE() else status); - end matchcontinue; + if stopOnFailure then + Error.addCompilerError("pre-optimization module " + moduleStr + " failed."); + fail(); + else + Error.addCompilerWarning("pre-optimization module " + moduleStr + " failed."); + end if; + end try; + end for; + + if Flags.isSet(Flags.OPT_DAE_DUMP) then + print("pre-optimization done.\n"); + end if; end preOptimizeDAE; public function transformBackendDAE " @@ -6933,56 +6927,49 @@ algorithm end matchcontinue; end dumpStrongComponents; -public function postOptimizeDAE " - Run the post-optimization modules." +public function postOptimizeDAE + "Run the post-optimization modules." input BackendDAE.BackendDAE inDAE; - input list> optModules; - input tuple matchingAlgorithm; - input tuple daeHandler; - output BackendDAE.BackendDAE outDAE; - output Util.Status status; + input list> inPostOptModules; + input tuple inMatchingAlgorithm; + input tuple inDAEHandler; + output BackendDAE.BackendDAE outDAE = inDAE; +protected + BackendDAEFunc.postOptimizationDAEModule optModule; + String moduleStr; + Boolean stopOnFailure; + BackendDAE.EqSystems systs; + BackendDAE.Shared shared; algorithm - (outDAE, status) := matchcontinue (inDAE, optModules, matchingAlgorithm, daeHandler) - local - BackendDAE.BackendDAE dae, dae1, dae2; - BackendDAEFunc.postOptimizationDAEModule optModule; - list> rest; - String str,moduleStr; - Boolean b; - BackendDAE.EqSystems systs; - BackendDAE.Shared shared; - - case (_, {}, _, _) - equation - if Flags.isSet(Flags.OPT_DAE_DUMP) then - print("post-optimization done.\n"); - end if; - then (inDAE,Util.SUCCESS()); - - case (_, (optModule, moduleStr, _)::rest, _, _) - equation - BackendDAE.DAE(systs, shared) = optModule(inDAE); - (systs, shared) = filterEmptySystems(systs, shared); - dae = BackendDAE.DAE(systs, shared); - SimCodeFunctionUtil.execStat("postOpt " + moduleStr); - if Flags.isSet(Flags.OPT_DAE_DUMP) then - print(stringAppendList({"\npost-optimization module ", moduleStr, ":\n\n"})); - BackendDump.printBackendDAE(dae); - end if; - dae1 = causalizeDAE(dae, NONE(), matchingAlgorithm, daeHandler, false); - (dae2, status) = postOptimizeDAE(dae1, rest, matchingAlgorithm, daeHandler); - then (dae2, status); + for postOptModule in inPostOptModules loop + (optModule, moduleStr, stopOnFailure) := postOptModule; + try + BackendDAE.DAE(systs, shared) := optModule(outDAE); + (systs, shared) := filterEmptySystems(systs, shared); + outDAE := BackendDAE.DAE(systs, shared); + outDAE := causalizeDAE(outDAE, NONE(), inMatchingAlgorithm, inDAEHandler, false); + SimCodeFunctionUtil.execStat("postOpt " + moduleStr); + if Flags.isSet(Flags.OPT_DAE_DUMP) then + print(stringAppendList({"\npost-optimization module ", moduleStr, ":\n\n"})); + BackendDump.printBackendDAE(outDAE); + end if; + else + SimCodeFunctionUtil.execStat(" postOpt " + moduleStr); + if stopOnFailure then + Error.addCompilerError("post-optimization module " + moduleStr + " failed."); + fail(); + else + Error.addCompilerWarning("post-optimization module " + moduleStr + " failed."); + end if; + end try; + end for; - case (_, (_, moduleStr, b)::rest, _, _) - equation - SimCodeFunctionUtil.execStat("postOpt " + moduleStr); - str = stringAppendList({"post-optimization module ", moduleStr, " failed."}); - Error.addMessage(Error.INTERNAL_ERROR, {str}); - (dae,status) = postOptimizeDAE(inDAE,rest,matchingAlgorithm,daeHandler); - then (dae, if b then Util.FAILURE() else status); - end matchcontinue; + if Flags.isSet(Flags.OPT_DAE_DUMP) then + print("post-optimization done.\n"); + end if; end postOptimizeDAE; + // protected function checkCompsMatching // "Function check if comps are complete, they are not complete // if the matching is wrong due to dummyDer" @@ -7082,14 +7069,14 @@ algorithm //fcall2(Flags.DUMP_DAE_LOW, BackendDump.dumpBackendDAE, inDAE, "dumpdaelow"); // pre optimisation phase _ := traverseBackendDAEExps(inDAE,ExpressionSimplify.simplifyTraverseHelper,0) "simplify all expressions"; - (optdae,Util.SUCCESS()) := preOptimizeDAE(inDAE,preOptModules); + optdae := preOptimizeDAE(inDAE,preOptModules); // transformation phase (matching and sorting using a index reduction method sode := causalizeDAE(optdae,NONE(),matchingAlgorithm,daeHandler,true); //fcall(Flags.DUMP_DAE_LOW, BackendDump.bltdump, ("bltdump",sode)); // post-optimization phase - (outSODE,Util.SUCCESS()) := postOptimizeDAE(sode,postOptModules,matchingAlgorithm,daeHandler); + outSODE := postOptimizeDAE(sode,postOptModules,matchingAlgorithm,daeHandler); _ := traverseBackendDAEExps(outSODE,ExpressionSimplify.simplifyTraverseHelper,0) "simplify all expressions"; //fcall2(Flags.DUMP_INDX_DAE, BackendDump.dumpBackendDAE, outSODE, "dumpindxdae"); @@ -7609,7 +7596,7 @@ algorithm nonEmpty := BackendVariable.varsSize(syst.orderedVars) <> 0 or BackendDAEUtil.equationArraySize(syst.removedEqs) <> 0; end nonEmptySystem; -public function filterEmptySystems +protected function filterEmptySystems "Filter out equation systems leaving at least one behind" input BackendDAE.EqSystems inSysts; input BackendDAE.Shared inShared; @@ -7617,16 +7604,15 @@ public function filterEmptySystems output BackendDAE.Shared outShared = inShared; protected list reqns; - BackendDAE.Equation eq; algorithm (reqns, outSysts) := List.fold(inSysts, filterEmptySystem, ({}, {})); - outSysts := match outSysts - local - BackendDAE.EqSystem syst; - case {} - then {BackendDAEUtil.createEqSystem(BackendVariable.emptyVars(), BackendEquation.emptyEqns())}; - else listReverseInPlace(outSysts); - end match; + + if listEmpty(outSysts) then + outSysts := {BackendDAEUtil.createEqSystem(BackendVariable.emptyVars(), BackendEquation.emptyEqns())}; + else + outSysts := listReverseInPlace(outSysts); + end if; + outShared.removedEqs := BackendEquation.addEquations(reqns, outShared.removedEqs); end filterEmptySystems; diff --git a/Compiler/BackEnd/HpcOmScheduler.mo b/Compiler/BackEnd/HpcOmScheduler.mo index 5078a48db62..946cb575c69 100644 --- a/Compiler/BackEnd/HpcOmScheduler.mo +++ b/Compiler/BackEnd/HpcOmScheduler.mo @@ -207,15 +207,15 @@ algorithm //print("Scheduling task " + intString(index) + " to thread " + intString(threadId) + "\n"); if(boolNot(listEmpty(predecessors))) then //in this case the node has predecessors - //find all predecessors which are scheduled to another thread and thus require a lock - (lockTasks,newOutgoingDepTasks) = iLockWithPredecessorHandler(head,predecessors,threadId,iCommCosts,iCompTaskMapping,iSimVarMapping); - outgoingDepTasks = listAppend(outgoingDepTasks,newOutgoingDepTasks); - //threadTasks = listAppend(List.map(newLockIdc,convertLockIdToAssignTask), threadTasks); - threadTasks = listAppend(lockTasks, threadTasks); - - //print("Eq idc: " + stringDelimitList(List.map(eqIdc, intString), ",") + "\n"); - simEqIdc = List.map(List.map1(eqIdc,getSimEqSysIdxForComp,iSccSimEqMapping), List.last); - //simEqIdc = List.sort(simEqIdc,intGt); + //find all predecessors which are scheduled to another thread and thus require a lock + (lockTasks,newOutgoingDepTasks) = iLockWithPredecessorHandler(head,predecessors,threadId,iCommCosts,iCompTaskMapping,iSimVarMapping); + outgoingDepTasks = listAppend(outgoingDepTasks,newOutgoingDepTasks); + //threadTasks = listAppend(List.map(newLockIdc,convertLockIdToAssignTask), threadTasks); + threadTasks = listAppend(lockTasks, threadTasks); + + //print("Eq idc: " + stringDelimitList(List.map(eqIdc, intString), ",") + "\n"); + simEqIdc = List.map(List.map1(eqIdc,getSimEqSysIdxForComp,iSccSimEqMapping), List.last); + //simEqIdc = List.sort(simEqIdc,intGt); else simEqIdc = List.flatten(List.map1(eqIdc,getSimEqSysIdxForComp,iSccSimEqMapping)); end if; @@ -459,7 +459,6 @@ protected algorithm HpcOmTaskGraph.TASKGRAPHMETA(commCosts=commCosts,inComps=inComps) := iTaskGraphMeta; taskGraphT := BackendDAEUtil.transposeMatrix(iTaskGraph,arrayLength(iTaskGraph)); - //() := HpcOmTaskGraph.printTaskGraph(taskGraphT); commCostsT := HpcOmTaskGraph.transposeCommCosts(commCosts); leaveNodes := HpcOmTaskGraph.getLeafNodes(iTaskGraph); //print("Leave nodes: " + stringDelimitList(List.map(leaveNodes,intString),", ") + "\n"); @@ -2499,7 +2498,7 @@ algorithm removeLocks = {}; tmpSchedule = HpcOmSimCode.THREADSCHEDULE(threadTasks,{},{},allCalcTasks); - (tmpSchedule,removeLocks) = createScheduleFromAssignments(extInfoArr,procAss,SOME(order),iTaskGraph,taskGraphT,iTaskGraphMeta,iSccSimEqMapping,removeLocks,order,commCosts,inComps,iSimVarMapping,tmpSchedule); + (tmpSchedule,removeLocks) = createScheduleFromAssignments(extInfoArr,procAss,SOME(order),iTaskGraph,taskGraphT,iTaskGraphMeta,iSccSimEqMapping,removeLocks,order,iSimVarMapping,tmpSchedule); // remove superfluous locks if Flags.isSet(Flags.HPCOM_DUMP) then print("number of removed superfluous locks: "+intString(intDiv(listLength(removeLocks),2))+"\n"); @@ -4551,112 +4550,189 @@ author: Waurich TUD 2015-02" input array> iSccSimEqMapping; input array> iSimVarMapping; //Maps each backend var to a list of simVars output HpcOmSimCode.Schedule oSchedule; -protected - Integer nProc,nTasks; - list rootNodes; - array threadMap; - array partitionCosts; - array> partitions; - HpcOmTaskGraph.TaskGraph graphT; - - array commCosts; - array> inComps; - array> threadTask; - array> allCalcTasks; - HpcOmSimCode.Schedule schedule; - list order; algorithm - HpcOmTaskGraph.TASKGRAPHMETA(commCosts=commCosts,inComps=inComps) := iTaskGraphMeta; - nProc := Flags.getConfigInt(Flags.NUM_PROC); - nTasks := arrayLength(iTaskGraph); - rootNodes := HpcOmTaskGraph.getRootNodes(iTaskGraph); - partitions := arrayCreate(nProc,{}); - threadMap := arrayCreate(nTasks,-1); - partitionCosts := arrayCreate(nProc,0.0); - graphT := BackendDAEUtil.transposeMatrix(iTaskGraph,arrayLength(iTaskGraph)); - ((partitions,threadMap,partitionCosts)) := List.fold(List.map(rootNodes,List.create),function assignTasksToPartitions(graph = iTaskGraph, graphT = graphT, meta = iTaskGraphMeta),(partitions,threadMap,partitionCosts)); + oSchedule := matchcontinue(iTaskGraph,iTaskGraphMeta,numProc,iSccSimEqMapping,iSimVarMapping) + local + Integer nTasks; + list rootNodes; + array taskMap; + array partitionCosts; + array> partitions, partMap; + HpcOmTaskGraph.TaskGraph graphT; - threadTask := arrayCreate(numProc,{}); - allCalcTasks := convertTaskGraphToTasks(graphT,iTaskGraphMeta,convertNodeToTask); - schedule := HpcOmSimCode.THREADSCHEDULE(threadTask,{},{},allCalcTasks); - order := List.intRange(nTasks); - (oSchedule,_) := createScheduleFromAssignments(threadMap,partitions,SOME(order),iTaskGraph,graphT,iTaskGraphMeta,iSccSimEqMapping,{},order,commCosts,inComps,iSimVarMapping,schedule); + array commCosts; + array> inComps; + array> threadTask; + array> allCalcTasks; + HpcOmSimCode.Schedule schedule; + list order; + case(_,HpcOmTaskGraph.TASKGRAPHMETA(commCosts=commCosts,inComps=inComps),_,_,_) + algorithm + true := intNe(arrayLength(iTaskGraph),0); + nTasks := arrayLength(iTaskGraph); + rootNodes := HpcOmTaskGraph.getRootNodes(iTaskGraph); + partitions := arrayCreate(numProc,{}); + taskMap := arrayCreate(nTasks,-1); + partMap := arrayCreate(listLength(rootNodes),{}); + partitionCosts := arrayCreate(numProc,0.0); + graphT := BackendDAEUtil.transposeMatrix(iTaskGraph,arrayLength(iTaskGraph)); + // get all existing partitions + (taskMap,partMap,_) := List.fold1(rootNodes,assignPartitions,iTaskGraph,(taskMap,partMap,1)); + //print("taskMap \n"+stringDelimitList(List.map(arrayList(taskMap), intString),"\n")+"\n"); + //print("partMap \n"+stringDelimitList(List.map(arrayList(partMap), HpcOmTaskGraph.intLstString),"\n")+"\n"); + // gather them to n partitions + (taskMap,partitions) := distributePartitions(taskMap,partMap,iTaskGraphMeta,numProc); + //print("partitions \n"+stringDelimitList(List.map(arrayList(partitions), HpcOmTaskGraph.intLstString),"\n")+"\n"); + + threadTask := arrayCreate(numProc,{}); + allCalcTasks := convertTaskGraphToTasks(graphT,iTaskGraphMeta,convertNodeToTask); + schedule := HpcOmSimCode.THREADSCHEDULE(threadTask,{},{},allCalcTasks); + order := List.flatten(HpcOmTaskGraph.getLevelNodes(iTaskGraph)); + if List.isEqual(arrayGet(partitions,1),{20,7,15,16,2},true) then + order := listReverse(order); + end if; + (oSchedule,_) := createScheduleFromAssignments(taskMap,partitions,SOME(order),iTaskGraph,graphT,iTaskGraphMeta,iSccSimEqMapping,{},order,iSimVarMapping,schedule); + then oSchedule; + case(_,_,_,_,_) + algorithm + true := intEq(arrayLength(iTaskGraph),0); + then HpcOmSimCode.EMPTYSCHEDULE(HpcOmSimCode.PARALLELTASKLIST({})); + else + algorithm + if Flags.isSet(Flags.FAILTRACE) then print("HpcOmScheduler.createPartSchedule failed\n"); end if; + then fail(); + end matchcontinue; end createPartSchedule; -protected function assignTasksToPartitions"assigns all rootNodes and their partitions to the partition with the least total execution costs" - input list rootNodes; +protected function distributePartitions + input array taskMapIn; + input array> partMap; + input HpcOmTaskGraph.TaskGraphMeta metaIn; + input Integer n; + output array taskMapOut; + output array> partitions; +protected + Integer partIdx; + Real costs; + list part; + list> clusters; + list partCosts={}; +algorithm + // get costs + for part in arrayList(partMap) loop + costs := List.fold(List.map1(part,HpcOmTaskGraph.getExeCostReqCycles,metaIn),realAdd,0.0); + partCosts := costs::partCosts; + end for; + partCosts := listReverse(partCosts); + //cluster them and correct task<->partition mapping + (partitions,_) := HpcOmTaskGraph.distributeToClusters(List.intRange(arrayLength(partMap)),partCosts,n); + for partIdx in List.intRange(n) loop + part := arrayGet(partitions,partIdx); + clusters := List.map1(part,Array.getIndexFirst,partMap); + part := List.fold(clusters,listAppend,{}); + partitions := arrayUpdate(partitions,partIdx,part); + List.map2_0(part, Array.updateIndexFirst,partIdx,taskMapIn); + end for; + taskMapOut := taskMapIn; +end distributePartitions; + +protected function assignPartitions"for every root node, assign all successing nodes to one partition. If we find an already assigned task from another partitions,replace all these tasks " + input Integer rootNode; input HpcOmTaskGraph.TaskGraph graph; - input HpcOmTaskGraph.TaskGraph graphT; - input HpcOmTaskGraph.TaskGraphMeta meta; - input tuple>,array,array> tplIn; - output tuple>,array,array> tplOut; + input tuple,array>,Integer> tplIn; // partitions, partitions-->tasks, currPartIdx> + output tuple,array>,Integer> tplOut; +protected + Integer node, idx; + array taskAss; + array> partAss; + list nodes, successors, assParts, unassTasks, otherParts, otherPartsTasks; +algorithm + (taskAss,partAss,idx) := tplIn; + taskAss := arrayUpdate(taskAss,rootNode,idx); + partAss := Array.appendToElement(idx,{rootNode},partAss); + nodes := {rootNode}; + while not listEmpty(nodes) loop + node::nodes := nodes; + successors := arrayGet(graph,node); + (unassTasks,otherPartsTasks) := List.split1OnTrue(successors,isUnAssigned,taskAss); + otherParts := List.map1(otherPartsTasks,Array.getIndexFirst,taskAss); + (otherParts,otherPartsTasks) := List.filter1OnTrueSync(otherParts,intNe,idx,otherPartsTasks); + otherParts := List.unique(otherParts); + if not listEmpty(otherParts) then + // if there are already tasks assigned to other partitions, replace these idxs + (taskAss,_) := Array.mapNoCopy_1(taskAss,reassignPartitions,(otherParts,idx)); + otherPartsTasks := List.fold(List.map1(otherParts,Array.getIndexFirst,partAss),listAppend,{}); // get all tasks that belong to the other partitions + List.map2_0(otherParts,Array.updateIndexFirst,{},partAss); + partAss := Array.appendToElement(idx,otherPartsTasks,partAss); + end if; + List.map2_0(unassTasks,Array.updateIndexFirst, idx, taskAss); + partAss := Array.appendToElement(idx,unassTasks,partAss); + nodes := listAppend(unassTasks,nodes); + end while; + tplOut := (taskAss,partAss,idx+1); +end assignPartitions; + +protected function isUnAssigned"checks whether the task is already assigned(==-1)" + input Integer task; + input array ass; + output Boolean isUnass; +protected + Integer idx; algorithm - tplOut := matchcontinue(rootNodes,graph,graphT,meta,tplIn) - local - Integer rootNode, partitionIdx; - Real minCosts,costs; - list rest, partition, assThreads; - list costLst; - array threadMap; - array partitionCosts; - array> partitions; - case({},_,_,_,(_,_,_)) - equation - then - tplIn; - case(rootNode::rest,_,_,_,(partitions,threadMap,partitionCosts)) - equation - true = intEq(-1,arrayGet(threadMap,rootNode)); - minCosts = Array.fold(partitionCosts,realMin,arrayGet(partitionCosts,1)); - partitionIdx = Array.position(partitionCosts,minCosts); - partition = getPartition({rootNode},graph,graphT,threadMap,partitionIdx,{}); - partitions = Array.appendToElement(partitionIdx,partition,partitions); - costLst = List.map(arrayGet(partitions,partitionIdx),function HpcOmTaskGraph.getExeCostReqCycles(iGraphData=meta)); - costs = List.fold(costLst,realAdd,0.0); - partitionCosts = arrayUpdate(partitionCosts,partitionIdx,costs); - (partitions,threadMap,partitionCosts) = assignTasksToPartitions(rest,graph,graphT,meta,(partitions,threadMap,partitionCosts)); - then - (partitions,threadMap,partitionCosts); - case(rootNode::rest,_,_,_,(partitions,threadMap,partitionCosts)) - equation - true = intNe(-1,arrayGet(threadMap,rootNode)); - (partitions,threadMap,partitionCosts) = assignTasksToPartitions(rest,graph,graphT,meta,(partitions,threadMap,partitionCosts)); - then - (partitions,threadMap,partitionCosts); + idx := arrayGet(ass,task); + isUnass := intEq(idx,-1); +end isUnAssigned; - end matchcontinue; -end assignTasksToPartitions; +protected function reassignPartitions"if the task is one of the oldAss, replace it with newAss" + input tuple,Integer>> tplIn; //value, + output tuple,Integer>> tplOut; +protected + Integer value, newAss; + list oldAss; +algorithm + (value,(oldAss,newAss)) := tplIn; + if List.exist1(oldAss,intEq,value) then + value := newAss; + end if; + tplOut := (value,(oldAss,newAss)); +end reassignPartitions; -protected function getPartition"get all tasks that are somehow connected to the checkNodes" - input list checkNodes; - input HpcOmTaskGraph.TaskGraph graph; - input HpcOmTaskGraph.TaskGraph graphT; - input array assNodes; - input Integer partitionIdx; - input list partitionIn; - output list partitionOut; +//--------------------------------- +// SingleThread Schedule +//--------------------------------- + +public function createSingleThreadSchedule"creates a schedule in which all tasks are computed in thread 1" + input HpcOmTaskGraph.TaskGraph iTaskGraph; + input HpcOmTaskGraph.TaskGraphMeta iTaskGraphMeta; + input array> iSccSimEqMapping; + input Integer numProc; + output HpcOmSimCode.Schedule oSchedule; +protected + Integer nTasks, size; + list order; + HpcOmTaskGraph.TaskGraph taskGraphT; + list allTasksLst={}; + array> thread2TaskAss; + array> allCalcTasks; algorithm - partitionOut := match(checkNodes,graph,graphT,assNodes,partitionIdx,partitionIn) - local - Integer node; - list children,parents,rest,partition; - case({},_,_,_,_,_) - then - partitionIn; - case(node::rest,_,_,_,_,_) - equation - children = arrayGet(graph,node); - (_,children) = List.filter1OnTrueSync(List.map(children,function Array.getIndexFirst(inArray = assNodes)),intEq,-1,children); - parents = arrayGet(graphT,node); - (_,parents) = List.filter1OnTrueSync(List.map(parents,function Array.getIndexFirst(inArray = assNodes)),intEq,-1,parents); - partition = listAppend(children, parents); - List.map2_0(node::partition,Array.updateIndexFirst, partitionIdx, assNodes); - rest = listAppend(partition,rest); - partition = listAppend(partitionIn,node::partition); - partition = getPartition(rest,graph,graphT,assNodes,partitionIdx,partition); - then partition; - end match; -end getPartition; + nTasks := arrayLength(iTaskGraph); + size := arrayLength(iTaskGraph); + taskGraphT := BackendDAEUtil.transposeMatrix(iTaskGraph,size); + // create the schedule + allCalcTasks := convertTaskGraphToTasks(taskGraphT,iTaskGraphMeta,convertNodeToTask); + + order := List.flatten(HpcOmTaskGraph.getLevelNodes(iTaskGraph)); + for i in order loop + // get the correct ordered tasks, replace the scc indexes with simEq indexes + allTasksLst := setSimEqIdcsInTask(Util.tuple21(arrayGet(allCalcTasks,i)),iSccSimEqMapping)::allTasksLst; + end for; + allTasksLst := listReverse(allTasksLst); + // set the thread Index + allTasksLst := List.map1(allTasksLst,setThreadIdxInTask,1); + thread2TaskAss := arrayCreate(numProc,{}); + thread2TaskAss := arrayUpdate(thread2TaskAss,1,allTasksLst); + oSchedule := HpcOmSimCode.THREADSCHEDULE(thread2TaskAss,{},{},allCalcTasks); +end createSingleThreadSchedule; //--------------------------------- @@ -4674,18 +4750,17 @@ author: Waurich TUD 2013-10 " output HpcOmSimCode.Schedule oSchedule; protected Integer size, numSfLocks; - array> threads; array> taskGraphT; array alapArray; // this is the latest possible starting time of every node list alapLst, alapSorted, priorityLst; list order; - list removeLocks; array taskAss; //=task, =processor array> procAss; //=processor, =task; - array> threadTask; HpcOmSimCode.Schedule schedule; - array> allCalcTasks; + list removeLocks; array commCosts; + array> threads, threadTask; + array> allCalcTasks; array> inComps; algorithm HpcOmTaskGraph.TASKGRAPHMETA(commCosts=commCosts,inComps=inComps) := iTaskGraphMeta; @@ -4703,7 +4778,7 @@ algorithm allCalcTasks := convertTaskGraphToTasks(taskGraphT,iTaskGraphMeta,convertNodeToTask); schedule := HpcOmSimCode.THREADSCHEDULE(threadTask,{},{},allCalcTasks); removeLocks := {}; - (schedule,removeLocks) := createScheduleFromAssignments(taskAss,procAss,SOME(order),iTaskGraph,taskGraphT,iTaskGraphMeta,iSccSimEqMapping,removeLocks,order,commCosts,inComps,iSimVarMapping,schedule); + (schedule,removeLocks) := createScheduleFromAssignments(taskAss,procAss,SOME(order),iTaskGraph,taskGraphT,iTaskGraphMeta,iSccSimEqMapping,removeLocks,order,iSimVarMapping,schedule); // remove superfluous locks numSfLocks := intDiv(listLength(removeLocks),2); if Flags.isSet(Flags.HPCOM_DUMP) then @@ -4860,14 +4935,12 @@ author:Waurich TUD 2013-12" input array> SccSimEqMappingIn; input list removeLocksIn; input list orderIn; // need the complete order for removeSuperfluousLocks - input array iCommCosts; - input array> iCompTaskMapping; //all StrongComponents from the BLT that belong to the Nodes [nodeId = arrayIdx] input array> iSimVarMapping; //Maps each backend var to a list of simVars input HpcOmSimCode.Schedule scheduleIn; output HpcOmSimCode.Schedule scheduleOut; output list removeLocksOut; algorithm - (scheduleOut,removeLocksOut) := match(taskAss,procAss,orderOpt,taskGraphIn,taskGraphTIn,taskGraphMetaIn,SccSimEqMappingIn,removeLocksIn,orderIn,iCommCosts,iCompTaskMapping,iSimVarMapping,scheduleIn) + (scheduleOut,removeLocksOut) := match(taskAss,procAss,orderOpt,taskGraphIn,taskGraphTIn,taskGraphMetaIn,SccSimEqMappingIn,removeLocksIn,orderIn,iSimVarMapping,scheduleIn) local Integer node,proc,mark,numProc; Real exeCost,commCost; @@ -4876,17 +4949,17 @@ algorithm array nodeMark; array> inComps; array> exeCosts; - array commCosts; + array inCommCosts; array> threadTasks; list taskLst1,taskLst,taskLstAss,taskLstRel, removeLocks; HpcOmSimCode.Schedule schedule; HpcOmSimCode.Task task; array> allCalcTasks; - case(_,_,SOME({}),_,_,_,_,_,_,_,_,_,HpcOmSimCode.THREADSCHEDULE()) + case(_,_,SOME({}),_,_,_,_,_,_,_,HpcOmSimCode.THREADSCHEDULE()) equation then (scheduleIn,removeLocksIn); - case(_,_,SOME(order),_,_,_,_,_,_,_,_,_,HpcOmSimCode.THREADSCHEDULE(threadTasks=threadTasks, outgoingDepTasks=outgoingDepTasks, allCalcTasks=allCalcTasks)) + case(_,_,SOME(order),_,_,HpcOmTaskGraph.TASKGRAPHMETA(commCosts=inCommCosts,inComps=inComps,nodeMark=nodeMark),_,_,_,_,HpcOmSimCode.THREADSCHEDULE(threadTasks=threadTasks, outgoingDepTasks=outgoingDepTasks, allCalcTasks=allCalcTasks)) equation numProc = arrayLength(procAss); (node::rest) = order; @@ -4899,12 +4972,12 @@ algorithm (_,otherParents,_) = List.intersection1OnTrue(parentNodes,sameProcTasks,intEq); (_,otherChildren,_) = List.intersection1OnTrue(childNodes,sameProcTasks,intEq); // keep the locks that are superfluous, remove them later - removeLocks = getSuperfluousLocks(otherParents,node,taskAss,orderIn,numProc,allCalcTasks,iCommCosts,iCompTaskMapping,iSimVarMapping,removeLocksIn); - taskLstAss = List.map6(otherParents,createDepTaskByTaskIdc,node,allCalcTasks,false,iCommCosts,iCompTaskMapping,iSimVarMapping); + removeLocks = getSuperfluousLocks(otherParents,node,taskAss,orderIn,numProc,allCalcTasks,inCommCosts,inComps,iSimVarMapping,removeLocksIn); + taskLstAss = List.map6(otherParents,createDepTaskByTaskIdc,node,allCalcTasks,false,inCommCosts,inComps,iSimVarMapping); //relLockDepTasks = List.map1(otherChildren,getReleaseLockString,node); - taskLstRel = List.map6(otherChildren,createDepTaskByTaskIdcR,node,allCalcTasks,true,iCommCosts,iCompTaskMapping,iSimVarMapping); + taskLstRel = List.map6(otherChildren,createDepTaskByTaskIdcR,node,allCalcTasks,true,inCommCosts,inComps,iSimVarMapping); + //build the calcTask - HpcOmTaskGraph.TASKGRAPHMETA(inComps=inComps,nodeMark=nodeMark) = taskGraphMetaIn; components = arrayGet(inComps,node); mark = arrayGet(nodeMark,node); ((_,exeCost)) = HpcOmTaskGraph.getExeCost(node,taskGraphMetaIn); @@ -4918,10 +4991,10 @@ algorithm threadTasks = arrayUpdate(threadTasks,proc,taskLst); outgoingDepTasks = listAppend(outgoingDepTasks,taskLstAss); schedule = HpcOmSimCode.THREADSCHEDULE(threadTasks,outgoingDepTasks,{},allCalcTasks); - (schedule,removeLocks) = createScheduleFromAssignments(taskAss,procAss,SOME(rest),taskGraphIn,taskGraphTIn,taskGraphMetaIn,SccSimEqMappingIn,removeLocks,orderIn,iCommCosts,iCompTaskMapping,iSimVarMapping,schedule); + (schedule,removeLocks) = createScheduleFromAssignments(taskAss,procAss,SOME(rest),taskGraphIn,taskGraphTIn,taskGraphMetaIn,SccSimEqMappingIn,removeLocks,orderIn,iSimVarMapping,schedule); then (schedule,removeLocks); - case(_,_,NONE(),_,_,_,_,_,_,_,_,_,HpcOmSimCode.THREADSCHEDULE()) + case(_,_,NONE(),_,_,_,_,_,_,_,HpcOmSimCode.THREADSCHEDULE()) equation print("createSchedulerFromAssignments failed.implement this!\n"); then @@ -4929,6 +5002,42 @@ algorithm end match; end createScheduleFromAssignments; +protected function setSimEqIdcsInTask"updates the eqIdcs from scc-Indexes to simEq-Indexes in calctasks " + input HpcOmSimCode.Task taskIn; + input array> SccSimEqMappingIn; + output HpcOmSimCode.Task taskOut; +algorithm + taskOut := matchcontinue(taskIn) + local + Integer weighting, index, threadIdx; + Real calcTime, timeFinished; + list eqIdc; + case(HpcOmSimCode.CALCTASK(weighting=weighting,index=index,calcTime=calcTime,timeFinished=timeFinished,threadIdx=threadIdx,eqIdc=eqIdc)) + equation + eqIdc = List.flatten(List.map1(eqIdc,getSimEqSysIdxForComp,SccSimEqMappingIn)); + then HpcOmSimCode.CALCTASK(weighting,index,calcTime,timeFinished,threadIdx,eqIdc); + else + then taskIn; + end matchcontinue; +end setSimEqIdcsInTask; + +protected function setThreadIdxInTask"updates threadIdxs in calctasks " + input HpcOmSimCode.Task taskIn; + input Integer threadIdx; + output HpcOmSimCode.Task taskOut; +algorithm + taskOut := matchcontinue(taskIn) + local + Integer weighting, index; + Real calcTime, timeFinished; + list eqIdc; + case(HpcOmSimCode.CALCTASK(weighting=weighting,index=index,calcTime=calcTime,timeFinished=timeFinished,eqIdc=eqIdc)) + then HpcOmSimCode.CALCTASK(weighting,index,calcTime,timeFinished,threadIdx,eqIdc); + else + then taskIn; + end matchcontinue; +end setThreadIdxInTask; + protected function tasksEqual "author: marcusw Checks if the given tasks are equal. The following conditions are checked: - if both tasks are of type CALCTASK: true if index equal @@ -5576,7 +5685,7 @@ author:Waurich TUD 2013-12" input String inSystemName; //e.g. "ODE system" or "DAE system" output String criticalPathInfoOut; algorithm - criticalPathInfoOut := match(scheduleIn,numProcIn,taskGraphIn,taskGraphMetaIn,inSystemName) + criticalPathInfoOut := matchcontinue(scheduleIn,numProcIn,taskGraphIn,taskGraphMetaIn,inSystemName) local list outgoingDepTasks; list levelCosts; @@ -5625,9 +5734,10 @@ algorithm criticalPathInfo; else equation + print("HpcOmScheduler.analyseScheduledTaskGraph failed\n"); then - ""; - end match; + "HpcOmScheduler.analyseScheduledTaskGraph failed\n"; + end matchcontinue; end analyseScheduledTaskGraph; protected function analyseScheduledTaskGraphLevel @@ -5889,7 +5999,7 @@ author:Waurich TUD 2013-11" output HpcOmSimCode.Schedule scheduleOut; output Real finishingTime; algorithm - (scheduleOut,finishingTime) := match(scheduleIn,numProc,taskGraphIn,taskGraphMetaIn) + (scheduleOut,finishingTime) := matchcontinue(scheduleIn,numProc,taskGraphIn,taskGraphMetaIn) local Real finTime; array taskIdcs; // idcs of the current Task for every proc. @@ -5925,7 +6035,11 @@ algorithm finTime = -1.0; then (schedule,finTime); - end match; + else + equation + print("getFinishingTimesForSchedule failed\n"); + then fail(); + end matchcontinue; end getFinishingTimesForSchedule; protected function getTimeFinishedOfLastTask "get the timeFinished of the last task of a thread. if the thread is empty its -1.0. diff --git a/Compiler/BackEnd/HpcOmTaskGraph.mo b/Compiler/BackEnd/HpcOmTaskGraph.mo index bd26d695272..1f8af618c7b 100644 --- a/Compiler/BackEnd/HpcOmTaskGraph.mo +++ b/Compiler/BackEnd/HpcOmTaskGraph.mo @@ -2314,10 +2314,61 @@ algorithm zeroFuncTaskGraphMeta := copyTaskGraphMeta(iTaskGraphMeta); zeroFuncTaskGraphMeta := setInCompsInMeta(zeroFuncInComps, zeroFuncTaskGraphMeta); - oTaskGraph := zeroFuncTaskGraph; - oTaskGraphMeta := zeroFuncTaskGraphMeta; + //reverse indexes + (oTaskGraph,oTaskGraphMeta) := reverseTaskGraphIndices(zeroFuncTaskGraph,zeroFuncTaskGraphMeta); end getZeroFuncsSystem; +protected function reverseTaskGraphIndices"reverse the task ids in the task grah and accordingly in the inComps +author Waurich TUD 07-2015" + input TaskGraph iTaskGraph; + input TaskGraphMeta iTaskGraphMeta; + output TaskGraph oTaskGraph; + output TaskGraphMeta oTaskGraphMeta; +protected + Integer nTasks; + array idxMap; + + array> inComps; + array> varCompMapping; + array> eqCompMapping; + array> compParamMapping; + array compNames; + array compDescs; + array> exeCosts; + array commCosts; + arraynodeMark; + array compInformations; +algorithm + nTasks := arrayLength(iTaskGraph); + idxMap := arrayCreate(nTasks,-1); + TASKGRAPHMETA(inComps=inComps, varCompMapping=varCompMapping, eqCompMapping=eqCompMapping, compParamMapping=compParamMapping, compNames=compNames, compDescs=compDescs, exeCosts=exeCosts, commCosts=commCosts, nodeMark=nodeMark, compInformations=compInformations) := iTaskGraphMeta; + // set an index mapping + for i in 1:nTasks loop + idxMap := arrayUpdate(idxMap,i,nTasks-i+1); + end for; + //map childNodes in taskgraph + oTaskGraph := Array.mapNoCopy_1(iTaskGraph,mapIntegers,idxMap); + oTaskGraph := Array.reverse(oTaskGraph); + inComps := Array.reverse(inComps); + oTaskGraphMeta := TASKGRAPHMETA(inComps, varCompMapping, eqCompMapping, compParamMapping, compNames, compDescs, exeCosts, commCosts, nodeMark, compInformations); +end reverseTaskGraphIndices; + +protected function mapIntegers" Array.mapNoCopy_1 - function to replace integers with their mapping integer +author Waurich TUD 07-2015" + input tuple,array> iTpl; + output tuple,array> oTpl; +protected + array map; + list iLst,oLst={}; +algorithm + (iLst,map) := iTpl; + for i in iLst loop + oLst := arrayGet(map,i)::oLst; + end for; + oLst := listReverse(oLst); + oTpl := (oLst,map); +end mapIntegers; + public function getEventSystem "gets the graph and the adjacencyLst only for the EventSystem. This means that all branches which leads to a node solving a whencondition or another boolean condition will remain. author: marcusw" @@ -3493,7 +3544,6 @@ algorithm oneChildren := listDelete(oneChildren,listLength(oneChildren)); // remove the empty startValue {} oneChildren := List.removeOnTrue(1,compareListLengthOnTrue,oneChildren); // remove paths of length 1 //print("oneChildren "+stringDelimitList(List.map(oneChildren,intLstString),"\n")+"\n"); - //(graphOut,graphTOut,graphDataOut,contractedTasksOut) := contractNodesInGraph(oneChildren,graphIn,graphTIn,graphDataIn,contractedTasksIn); (graphOut,graphTOut,graphDataOut,contractedTasksOut) := contractNodesInGraph(oneChildren,graphIn,graphTIn,graphDataIn,contractedTasksIn); changed := not listEmpty(oneChildren); //print("contractedTasksOut "+stringDelimitList(List.map(arrayList(contractedTasksOut),intString),"\n")+"\n"); diff --git a/Compiler/BackEnd/IndexReduction.mo b/Compiler/BackEnd/IndexReduction.mo index 89b357ae1f0..aabab0d09f8 100644 --- a/Compiler/BackEnd/IndexReduction.mo +++ b/Compiler/BackEnd/IndexReduction.mo @@ -2492,6 +2492,7 @@ protected Integer nstatevars,nassigned,nunassigned,nass1arr,n,nv,ne; StateSets stateSets; algorithm + try for seteqns in iSets loop if not listEmpty(List.select1r(seteqns,Matching.isUnAssigned,vec1)) then // ignore sets without unassigned equations, because all assigned states already in dummy states // print("seteqns: " + intString(listLength(seteqns)) + "\n"); @@ -2546,7 +2547,10 @@ algorithm outDummyVars := listAppend(varlst, outDummyVars); end if; end for; - + else + Error.addMessage(Error.INTERNAL_ERROR, {"- IndexReduction.processComps4New failed!"}); + fail(); + end try; end processComps4New; protected function forceInlinEqn @@ -3106,6 +3110,9 @@ algorithm unassignedEqnsSize = listLength(unassignedEqns); size = listLength(states); rang = size-unassignedEqnsSize; + if intLt(rang,0) then + Error.addMessage(Error.INTERNAL_ERROR, {"Selection of DummyDerivatives failed due to negative system rank of "+intString(rang)+"! + There are "+intString(unassignedEqnsSize)+" unassigned equations and "+intString(size)+" potential states.\n"}); end if; true = intEq(rang,0); if Flags.isSet(Flags.BLT_DUMP) then print("Select as dummyStates(3):\n"); @@ -3115,6 +3122,11 @@ algorithm varlst = List.map1r(List.map(states,Util.tuple22),BackendVariable.getVarAt,vars); then (varlst,iStateSets); + else + equation + Error.addMessage(Error.INTERNAL_ERROR, {"- IndexReduction.selectDummyDerivatives2new failed!"}); + then + fail(); end matchcontinue; end selectDummyDerivatives2new; diff --git a/Compiler/BackEnd/Initialization.mo b/Compiler/BackEnd/Initialization.mo index e4ea7c329e5..66a05b8cdf8 100644 --- a/Compiler/BackEnd/Initialization.mo +++ b/Compiler/BackEnd/Initialization.mo @@ -180,7 +180,7 @@ algorithm initdae := BackendDAEUtil.transformBackendDAE(initdae, SOME((BackendDAE.NO_INDEX_REDUCTION(), BackendDAE.EXACT())), NONE(), NONE()); // simplify system - (initdae, Util.SUCCESS()) := BackendDAEUtil.postOptimizeDAE(initdae, pastOptModules, matchingAlgorithm, daeHandler); + initdae := BackendDAEUtil.postOptimizeDAE(initdae, pastOptModules, matchingAlgorithm, daeHandler); if Flags.isSet(Flags.DUMP_INITIAL_SYSTEM) then BackendDump.dumpBackendDAE(initdae, "solved initial system"); if Flags.isSet(Flags.ADDITIONAL_GRAPHVIZ_DUMP) then diff --git a/Compiler/SimCode/HpcOmSimCodeMain.mo b/Compiler/SimCode/HpcOmSimCodeMain.mo index 651a03a2d0c..81905bb9733 100644 --- a/Compiler/SimCode/HpcOmSimCodeMain.mo +++ b/Compiler/SimCode/HpcOmSimCodeMain.mo @@ -349,19 +349,19 @@ algorithm (numProc,_) = setNumProc(numProc,cpCostsWoC,taskGraphDataOde);//in case n-flag is not set (scheduleDae,simCode,taskGraphDaeScheduled,taskGraphDataDaeScheduled,sccSimEqMapping) = createSchedule(taskGraphDaeSimplified,taskGraphDataDaeSimplified,daeSccSimEqMapping,simVarMapping,filenamePrefix,numProc,simCode,scheduledTasksDae,"DAE system",Flags.getConfigString(Flags.HPCOM_SCHEDULER)); - //criticalPathInfo = HpcOmScheduler.analyseScheduledTaskGraph(scheduleDae,numProc,taskGraphDaeScheduled,taskGraphDataDaeScheduled,"DAE system"); + //HpcOmScheduler.printSchedule(scheduleDae); //schedulerInfo = HpcOmScheduler.convertScheduleStrucToInfo(scheduleDae,arrayLength(taskGraphDaeScheduled)); - //HpcOmTaskGraph.dumpAsGraphMLSccLevel(taskGraphDaeScheduled, taskGraphDataDaeScheduled, fileName+"_schedDAE.graphml", criticalPathInfo, HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPaths)), HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPathsWoC)), sccSimEqMapping, schedulerInfo, HpcOmTaskGraph.GRAPHDUMPOPTIONS(false,false,false,false)); + //HpcOmTaskGraph.dumpAsGraphMLSccLevel(taskGraphDaeScheduled, taskGraphDataDaeScheduled, "taskGraph"+filenamePrefix+"DAE_scheduled.graphml", "", HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPaths)), HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPathsWoC)), sccSimEqMapping, schedulerInfo, HpcOmTaskGraph.GRAPHDUMPOPTIONS(false,false,false,false)); (scheduleOde,simCode,taskGraphOdeScheduled,taskGraphDataOdeScheduled,sccSimEqMapping) = createSchedule(taskGraphOdeSimplified,taskGraphDataOdeSimplified,sccSimEqMapping,simVarMapping,filenamePrefix,numProc,simCode,scheduledTasksOde,"ODE system",Flags.getConfigString(Flags.HPCOM_SCHEDULER)); - //criticalPathInfo = HpcOmScheduler.analyseScheduledTaskGraph(scheduleOde,numProc,taskGraphOdeScheduled,taskGraphDataOdeScheduled,"DAE system"); + //HpcOmScheduler.printSchedule(scheduleOde); //schedulerInfo = HpcOmScheduler.convertScheduleStrucToInfo(scheduleOde,arrayLength(taskGraphOdeScheduled)); - //HpcOmTaskGraph.dumpAsGraphMLSccLevel(taskGraphOdeScheduled, taskGraphDataOdeScheduled, fileName+"_schedODE.graphml", criticalPathInfo, HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPaths)), HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPathsWoC)), sccSimEqMapping, schedulerInfo, HpcOmTaskGraph.GRAPHDUMPOPTIONS(false,false,false,false)); + //HpcOmTaskGraph.dumpAsGraphMLSccLevel(taskGraphOdeScheduled, taskGraphDataOdeScheduled, "taskGraph"+filenamePrefix+"ODE_scheduled.graphml", "", HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPaths)), HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPathsWoC)), sccSimEqMapping, schedulerInfo, HpcOmTaskGraph.GRAPHDUMPOPTIONS(false,false,false,false)); (scheduleZeroFunc,simCode,taskGraphZeroFuncScheduled,taskGraphDataZeroFuncScheduled,sccSimEqMapping) = createSchedule(taskGraphZeroFuncSimplified,taskGraphDataZeroFuncSimplified,daeSccSimEqMapping,simVarMapping,filenamePrefix,numProc,simCode,scheduledTasksZeroFunc,"ZeroFunc system",Flags.getConfigString(Flags.HPCOM_SCHEDULER)); - //criticalPathInfo = HpcOmScheduler.analyseScheduledTaskGraph(scheduleZeroFunc,numProc,taskGraphZeroFuncScheduled,taskGraphDataZeroFuncScheduled,"DAE system"); - //schedulerInfo = HpcOmScheduler.convertScheduleStrucToInfo(scheduleZeroFunc,arrayLength(taskGraphZeroFuncScheduled)); - //HpcOmTaskGraph.dumpAsGraphMLSccLevel(taskGraphZeroFuncScheduled, taskGraphDataZeroFuncScheduled, fileName+"_schedZE.graphml", criticalPathInfo, HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPaths)), HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPathsWoC)), sccSimEqMapping, schedulerInfo, HpcOmTaskGraph.GRAPHDUMPOPTIONS(false,false,false,false)); + //HpcOmScheduler.printSchedule(scheduleZeroFunc); + //schedulerInfo = HpcOmScheduler.convertScheduleStrucToInfo(scheduleZeroFunc,arrayLength(taskGraphZeroFuncScheduled)); + //HpcOmTaskGraph.dumpAsGraphMLSccLevel(taskGraphZeroFuncScheduled, taskGraphDataZeroFuncScheduled, "taskGraph"+filenamePrefix+"ZF_scheduled.graphml", "", HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPaths)), HpcOmTaskGraph.convertNodeListToEdgeTuples(listHead(criticalPathsWoC)), sccSimEqMapping, schedulerInfo, HpcOmTaskGraph.GRAPHDUMPOPTIONS(false,false,false,false)); SimCode.SIMCODE( modelInfo, simCodeLiterals, simCodeRecordDecls, simCodeExternalFunctionIncludes, allEquations, odeEquations, algebraicEquations, partitionsKind, baseClocks, useHomotopy, initialEquations, removedInitialEquations, startValueEquations, @@ -779,7 +779,7 @@ author: mwalther, Waurich TUD" output HpcOmTaskGraph.TaskGraphMeta oTaskGraphMeta; output array> oSccSimEqMapping; protected - list knownScheduler = {"none","level","levelfix","ext","metis","hmet","listr","rand","list","mcp","part","taskdep","tds","bls","sbs"}; + list knownScheduler = {"none","level","levelfix","ext","metis","hmet","listr","rand","list","mcp","part","taskdep","tds","bls","sbs","sts"}; String schedulerName = iSchedulerName; algorithm if(boolNot(List.exist1(knownScheduler,stringEq,schedulerName))) then @@ -891,6 +891,11 @@ algorithm print("Using Single Block Scheduling for the " + iSystemName + "\n"); schedule = HpcOmEqSystems.createSingleBlockSchedule(iTaskGraph,iTaskGraphMeta,iScheduledTasks,iSccSimEqMapping); then (schedule,iSimCode,iTaskGraph,iTaskGraphMeta,iSccSimEqMapping); + case(_,_,_,_,_,_,_,_,_,"sts") + equation + print("Using Single Thread Scheduling for the " + iSystemName + "\n"); + schedule = HpcOmScheduler.createSingleThreadSchedule(iTaskGraph,iTaskGraphMeta,iSccSimEqMapping,iNumProc); + then (schedule,iSimCode,iTaskGraph,iTaskGraphMeta,iSccSimEqMapping); else equation print("HpcOmSimCode.createSchedule failed!\n"); diff --git a/Compiler/Util/Array.mo b/Compiler/Util/Array.mo index 825ac7ad06b..6619a08fa75 100644 --- a/Compiler/Util/Array.mo +++ b/Compiler/Util/Array.mo @@ -710,5 +710,22 @@ algorithm fail(); end getMemberOnTrue; +public function reverse"reverses the elements in an array" + input array inArray; + output array outArray; +protected + Integer size,i; + T elem1,elem2; +algorithm + outArray := inArray; + size := arrayLength(inArray); + for i in 1:(size/2) loop + elem1 := arrayGet(inArray,i); + elem2 := arrayGet(inArray,size-i+1); + outArray := arrayUpdate(outArray,i,elem2); + outArray := arrayUpdate(outArray,size-i+1,elem1); + end for; +end reverse; + annotation(__OpenModelica_Interface="util"); end Array; diff --git a/Compiler/Util/Flags.mo b/Compiler/Util/Flags.mo index 8000374e885..1a26eced183 100644 --- a/Compiler/Util/Flags.mo +++ b/Compiler/Util/Flags.mo @@ -802,7 +802,7 @@ constant ConfigFlag POST_OPT_MODULES = CONFIG_FLAG(16, "postOptModules", ("tearingSystem",Util.notrans("For method selection use flag tearingMethod.")), ("partlintornsystem",Util.notrans("partitions linear torn systems.")), ("relaxSystem",Util.notrans("DESCRIBE ME")), - ("countOperations", Util.gettext("Count the mathematic operations of the system.")), + ("countOperations", Util.gettext("Count the mathematical operations of the system.")), ("dumpComponentsGraphStr", Util.notrans("DESCRIBE ME")), ("generateSymbolicJacobian", Util.gettext("Generates symbolic Jacobian matrix, where der(x) is differentiated w.r.t. x. This matrix can be used to simulate with dasslColorSymJac.")), ("generateSymbolicLinearization", Util.gettext("Generates symbolic linearization matrices A,B,C,D for linear model:\n\t\t:math:`\\dot{x} = Ax + Bu`\n\t:math:`ty = Cx +Du`")), @@ -816,7 +816,7 @@ constant ConfigFlag POST_OPT_MODULES = CONFIG_FLAG(16, "postOptModules", ("calculateStateSetsJacobians", Util.gettext("Generates analytical Jacobian for dynamic state selection sets.")), ("addInitialStmtsToAlgorithms", Util.gettext("Expands all algorithms with initial statements for outputs.")), ("reshufflePost", Util.gettext("Reshuffles algebraic loops.")), - ("CSE", Util.gettext("Common Subexpression Elimination")), + ("CSE", Util.gettext("Common Sub-expression Elimination")), ("dumpDAE", Util.gettext("dumps the DAE representation of the current transformation state")), ("dumpDAEXML", Util.gettext("dumps the DAE as xml representation of the current transformation state")), ("addTimeAsState", Util.gettext("Experimental feature: this replaces each occurrence of variable time with a new introduced state $time with equation der($time) = 1.0")) @@ -1088,7 +1088,7 @@ constant ConfigFlag LOOP2CON = CONFIG_FLAG(71, "loop2con", constant ConfigFlag FORCE_TEARING = CONFIG_FLAG(72, "forceTearing", NONE(), EXTERNAL(), BOOL_FLAG(false), NONE(), - Util.gettext("Use tearing set even if it is not smaller than the original component.)")); + Util.gettext("Use tearing set even if it is not smaller than the original component.")); constant ConfigFlag SIMPLIFY_LOOPS = CONFIG_FLAG(73, "simplifyLoops", NONE(), EXTERNAL(), INT_FLAG(0), diff --git a/SimulationRuntime/c/Makefile.objs b/SimulationRuntime/c/Makefile.objs index ce26ab452d9..cf7ca3ad715 100644 --- a/SimulationRuntime/c/Makefile.objs +++ b/SimulationRuntime/c/Makefile.objs @@ -21,7 +21,7 @@ gc/mmc_gc.h # Files for util functions UTIL_OBJS_MINIMAL=base_array$(OBJ_EXT) boolean_array$(OBJ_EXT) omc_error$(OBJ_EXT) division$(OBJ_EXT) generic_array$(OBJ_EXT) index_spec$(OBJ_EXT) integer_array$(OBJ_EXT) list$(OBJ_EXT) memory_pool$(OBJ_EXT) modelica_string$(OBJ_EXT) read_write$(OBJ_EXT) write_matlab4$(OBJ_EXT) read_matlab4$(OBJ_EXT) real_array$(OBJ_EXT) ringbuffer$(OBJ_EXT) rtclock$(OBJ_EXT) string_array$(OBJ_EXT) utility$(OBJ_EXT) varinfo$(OBJ_EXT) ModelicaUtilities$(OBJ_EXT) omc_msvc$(OBJ_EXT) simulation_options$(OBJ_EXT) tinymt64$(OBJ_EXT) omc_mmap$(OBJ_EXT) cJSON$(OBJ_EXT) ifeq ($(OMC_MINIMAL_RUNTIME),) -UTIL_OBJS=$(UTIL_OBJS_MINIMAL) java_interface$(OBJ_EXT) libcsv$(OBJ_EXT) read_csv$(OBJ_EXT) OldModelicaTables$(OBJ_EXT) +UTIL_OBJS=$(UTIL_OBJS_MINIMAL) java_interface$(OBJ_EXT) libcsv$(OBJ_EXT) read_csv$(OBJ_EXT) OldModelicaTables$(OBJ_EXT) write_csv$(OBJ_EXT) else UTIL_OBJS=$(UTIL_OBJS_MINIMAL) endif @@ -33,7 +33,7 @@ MATH_HFILES = blaswrap.h SOLVER_OBJS_MINIMAL=delay$(OBJ_EXT) events$(OBJ_EXT) external_input$(OBJ_EXT) linearSystem$(OBJ_EXT) linearSolverLapack$(OBJ_EXT) linearSolverTotalPivot$(OBJ_EXT) mixedSystem$(OBJ_EXT) mixedSearchSolver$(OBJ_EXT) model_help$(OBJ_EXT) nonlinearSystem$(OBJ_EXT) nonlinearSolverHomotopy$(OBJ_EXT) omc_math$(OBJ_EXT) solver_main$(OBJ_EXT) stateset$(OBJ_EXT) ifeq ($(OMC_MINIMAL_RUNTIME),) -SOLVER_OBJS=$(SOLVER_OBJS_MINIMAL) kinsolSolver$(OBJ_EXT) linearSolverLis$(OBJ_EXT) linearSolverUmfpack$(OBJ_EXT) dassl$(OBJ_EXT) radau$(OBJ_EXT) nonlinearSolverHybrd$(OBJ_EXT) nonlinearSolverNewton$(OBJ_EXT) +SOLVER_OBJS=$(SOLVER_OBJS_MINIMAL) kinsolSolver$(OBJ_EXT) linearSolverLis$(OBJ_EXT) linearSolverUmfpack$(OBJ_EXT) dassl$(OBJ_EXT) radau$(OBJ_EXT) nonlinearSolverHybrd$(OBJ_EXT) nonlinearSolverNewton$(OBJ_EXT) newtonIteration$(OBJ_EXT) else SOLVER_OBJS=$(SOLVER_OBJS_MINIMAL) endif diff --git a/SimulationRuntime/c/simulation/simulation_runtime.cpp b/SimulationRuntime/c/simulation/simulation_runtime.cpp index 71e02466e83..1a2d5126ec4 100644 --- a/SimulationRuntime/c/simulation/simulation_runtime.cpp +++ b/SimulationRuntime/c/simulation/simulation_runtime.cpp @@ -773,6 +773,7 @@ int initRuntimeAndSimulation(int argc, char**argv, DATA *data) data->simulationInfo.nlsMethod = getNonlinearSolverMethod(argc, argv); data->simulationInfo.lsMethod = getlinearSolverMethod(argc, argv); data->simulationInfo.newtonStrategy = getNewtonStrategy(argc, argv); + data->simulationInfo.nlsCsvInfomation = omc_flag[FLAG_NLS_INFO]; rt_tick(SIM_TIMER_INIT_XML); read_input_xml(&(data->modelData), &(data->simulationInfo)); diff --git a/SimulationRuntime/c/simulation/solver/CMakeLists.txt b/SimulationRuntime/c/simulation/solver/CMakeLists.txt index 9d6f9400e1c..0d973ad2a66 100644 --- a/SimulationRuntime/c/simulation/solver/CMakeLists.txt +++ b/SimulationRuntime/c/simulation/solver/CMakeLists.txt @@ -11,15 +11,15 @@ SET(solver_sources ../../../../3rdParty/Cdaskr/solver/daux.c ../../../../3rdParty/Cdaskr/solver/dlinpk.c dassl.c kinsolSolver.c linearSystem.c nonlinearSolverHybrd.c radau.c -delay.c linearSolverLapack.c mixedSearchSolver.c nonlinearSolverNewton.c solver_main.c -linearSolverLis.c mixedSystem.c nonlinearSystem.c stateset.c +delay.c linearSolverLapack.c mixedSearchSolver.c nonlinearSolverNewton.c newtonIteration.c solver_main.c +linearSolverLis.c mixedSystem.c nonlinearSystem.c stateset.c events.c linearSolverTotalPivot.c model_help.c omc_math.c external_input.c linearSolverUmfpack.c nonlinearSolverHomotopy.c) SET(solver_headers ../../../../3rdParty/Cdaskr/solver/ddaskr_types.h dassl.h external_input.h linearSolverUmfpack.h nonlinearSolverHomotopy.h radau.h delay.h kinsolSolver.h linearSystem.h nonlinearSolverHybrd.h solver_main.h -linearSolverLapack.h mixedSearchSolver.h nonlinearSolverNewton.h stateset.h +linearSolverLapack.h mixedSearchSolver.h nonlinearSolverNewton.h newtonIteration.h stateset.h epsilon.h linearSolverLis.h mixedSystem.h nonlinearSystem.h events.h linearSolverTotalPivot.h model_help.h omc_math.h) diff --git a/SimulationRuntime/c/simulation/solver/model_help.c b/SimulationRuntime/c/simulation/solver/model_help.c index 28973b59342..82e733588fd 100644 --- a/SimulationRuntime/c/simulation/solver/model_help.c +++ b/SimulationRuntime/c/simulation/solver/model_help.c @@ -907,6 +907,7 @@ void initializeDataStruc(DATA *data) data->simulationInfo.lsMethod = LS_LAPACK; data->simulationInfo.mixedMethod = MIXED_SEARCH; data->simulationInfo.newtonStrategy = NEWTON_PURE; + data->simulationInfo.nlsCsvInfomation = 0; data->simulationInfo.zeroCrossings = (modelica_real*) calloc(data->modelData.nZeroCrossings, sizeof(modelica_real)); data->simulationInfo.zeroCrossingsPre = (modelica_real*) calloc(data->modelData.nZeroCrossings, sizeof(modelica_real)); diff --git a/SimulationRuntime/c/simulation/solver/newtonIteration.c b/SimulationRuntime/c/simulation/solver/newtonIteration.c new file mode 100644 index 00000000000..dc3f366bb85 --- /dev/null +++ b/SimulationRuntime/c/simulation/solver/newtonIteration.c @@ -0,0 +1,774 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2014, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE + * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: http://www.openmodelica.org or + * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica + * distribution. GNU version 3 is obtained from: + * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from: + * http://www.opensource.org/licenses/BSD-3-Clause. + * + * This program is distributed WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS + * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE + * CONDITIONS OF OSMC-PL. + * + */ + +/*! \file newtonIteration.c + */ + +#ifdef __cplusplus +extern "C" { +#endif + +#include +#include +#include /* memcpy */ + +#include "simulation/simulation_info_xml.h" +#include "util/omc_error.h" +#include "util/varinfo.h" +#include "model_help.h" + +#include "nonlinearSystem.h" +#include "newtonIteration.h" + +#include "external_input.h" + + +extern double enorm_(int *n, double *x); +int solveLinearSystem(int* n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData); +void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f, + double* scaledError_f, int* n, double* x, double* fvec); +void scaling_residual_vector(DATA_NEWTON* solverData); +void damping_heuristic(double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, double* lambda, int* k, DATA_NEWTON* solverData, void* userdata); +void damping_heuristic2(double damping_parameter, double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata); +void LineSearch(double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata); +void Backtracking(double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, DATA_NEWTON* solverData, void* userdata); +void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps); + + +#ifdef __cplusplus +extern "C" { +#endif + +extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info); +extern void dgetrf_(int *m, int *n, doublereal *fjac, int *lda, int* iwork, int *info); +extern void dgetrs_(char *trans, int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info); + +#ifdef __cplusplus +} +#endif + + +/*! \fn allocateNewtonData + * allocate memory for nonlinear system solver + */ +int allocateNewtonData(int size, void** voiddata) +{ + DATA_NEWTON* data = (DATA_NEWTON*) malloc(sizeof(DATA_NEWTON)); + + *voiddata = (void*)data; + assertStreamPrint(NULL, NULL != data, "allocationNewtonData() failed!"); + + data->resScaling = (double*) malloc(size*sizeof(double)); + data->fvecScaled = (double*) malloc(size*sizeof(double)); + + data->n = size; + data->x = (double*) malloc(size*sizeof(double)); + data->fvec = (double*) calloc(size,sizeof(double)); + data->xtol = 1e-6; + data->ftol = 1e-6; + data->maxfev = size*100; + data->epsfcn = DBL_EPSILON; + data->fjac = (double*) malloc((size*size)*sizeof(double)); + + data->rwork = (double*) malloc((size)*sizeof(double)); + data->iwork = (int*) malloc(size*sizeof(int)); + + /* damped newton */ + data->x_new = (double*) malloc(size*sizeof(double)); + data->x_increment = (double*) malloc(size*sizeof(double)); + data->f_old = (double*) calloc(size,sizeof(double)); + data->fvec_minimum = (double*) calloc(size,sizeof(double)); + data->delta_f = (double*) calloc(size,sizeof(double)); + data->delta_x_vec = (double*) calloc(size,sizeof(double)); + + data->factorization = 0; + data->calculate_jacobian = 1; + data->numberOfIterations = 0; + data->numberOfFunctionEvaluations = 0; + + return 0; +} + +/*! \fn freeNewtonData + * + * free memory for nonlinear solver newton + * + */ +int freeNewtonData(void **voiddata) +{ + DATA_NEWTON* data = (DATA_NEWTON*) *voiddata; + + free(data->resScaling); + free(data->fvecScaled); + free(data->x); + free(data->fvec); + free(data->fjac); + free(data->rwork); + free(data->iwork); + + /* damped newton */ + free(data->x_new); + free(data->x_increment); + free(data->f_old); + free(data->fvec_minimum); + free(data->delta_f); + free(data->delta_x_vec); + + return 0; +} + +/*! \fn solve system with Newton-Raphson + * + * \param [in] [n] size of equation + * [eps] tolerance for x + * [h] tolerance for f' + * [k] maximum number of iterations + * [work] work array size of (n*X) + * [f] user provided function + * [data] userdata + * [info] + * [calculate_jacobian] flag which decides whether Jacobian is calculated + * (0) once for the first calculation + * (i) every i steps (=1 means original newton method) + * (-1) never, factorization has to be given in A + * + */ +int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solverData, void* userdata) +{ + + int i, j, k = 0, l = 0, nrsh = 1; + int *n = &(solverData->n); + double *x = solverData->x; + double *fvec = solverData->fvec; + double *eps = &(solverData->ftol); + double *fdeps = &(solverData->epsfcn); + int * maxfev = &(solverData->maxfev); + double *fjac = solverData->fjac; + double *work = solverData->rwork; + int *iwork = solverData->iwork; + int *info = &(solverData->info); + int calc_jac = 1; + + double error_f = 1.0 + *eps, scaledError_f = 1.0 + *eps, delta_x = 1.0 + *eps, delta_f = 1.0 + *eps, delta_x_scaled = 1.0 + *eps, lambda = 1.0; + double current_fvec_enorm, enorm_new; + + + if(ACTIVE_STREAM(LOG_NLS_V)) + { + infoStreamPrint(LOG_NLS_V, 1, "######### Start Newton maxfev: %d #########", (int)*maxfev); + + infoStreamPrint(LOG_NLS_V, 1, "x vector"); + for(i=0; i<*n; i++) + infoStreamPrint(LOG_NLS_V, 0, "x[%d]: %e ", i, x[i]); + messageClose(LOG_NLS_V); + + messageClose(LOG_NLS_V); + } + + *info = 1; + + /* calculate the function values */ + (*f)(n, x, fvec, userdata, 1); + + solverData->nfev++; + + /* save current fvec in f_old*/ + memcpy(solverData->f_old, fvec, *n*sizeof(double)); + + error_f = current_fvec_enorm = enorm_(n, fvec); + + while(error_f > *eps && scaledError_f > *eps && delta_x > *eps && delta_f > *eps && delta_x_scaled > *eps) + { + if(ACTIVE_STREAM(LOG_NLS_V)) + { + infoStreamPrint(LOG_NLS_V, 0, "\n**** start Iteration: %d *****", (int) l); + + /* Debug output */ + infoStreamPrint(LOG_NLS_V, 1, "function values"); + for(i=0; i<*n; i++) + infoStreamPrint(LOG_NLS_V, 0, "fvec[%d]: %e ", i, fvec[i]); + messageClose(LOG_NLS_V); + } + + /* calculate jacobian if no matrix is given */ + if (calc_jac == 1 && solverData->calculate_jacobian >= 0) + { + (*f)(n, x, fvec, userdata, 0); + solverData->factorization = 0; + calc_jac = solverData->calculate_jacobian; + } + else + { + solverData->factorization = 1; + calc_jac--; + } + + + /* debug output */ + if(ACTIVE_STREAM(LOG_NLS_JAC)) + { + char buffer[4096]; + + infoStreamPrint(LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", (int)*n, (int)*n); + for(i=0; in;i++) + { + buffer[0] = 0; + for(j=0; jn; j++) + sprintf(buffer, "%s%10g ", buffer, fjac[i*(*n)+j]); + infoStreamPrint(LOG_NLS_JAC, 0, "%s", buffer); + } + messageClose(LOG_NLS_JAC); + } + + if (solveLinearSystem(n, iwork, fvec, fjac, solverData) != 0) + { + *info=-1; + break; + } + else + { + for (i =0; i<*n; i++) + solverData->x_new[i]=x[i]-solverData->x_increment[i]; + + infoStreamPrint(LOG_NLS_V,1,"x_increment"); + for(i=0; i<*n; i++) + infoStreamPrint(LOG_NLS_V, 0, "x_increment[%d] = %e ", i, solverData->x_increment[i]); + messageClose(LOG_NLS_V); + + if (solverData->newtonStrategy == NEWTON_DAMPED) + { + damping_heuristic(x, f, current_fvec_enorm, n, fvec, &lambda, &k, solverData, userdata); + } + else if (solverData->newtonStrategy == NEWTON_DAMPED2) + { + damping_heuristic2(0.75, x, f, current_fvec_enorm, n, fvec, &k, solverData, userdata); + } + else if (solverData->newtonStrategy == NEWTON_DAMPED_LS) + { + LineSearch(x, f, current_fvec_enorm, n, fvec, &k, solverData, userdata); + } + else if (solverData->newtonStrategy == NEWTON_DAMPED_BT) + { + Backtracking(x, f, current_fvec_enorm, n, fvec, solverData, userdata); + } + else + { + /* calculate the function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + } + + calculatingErrors(solverData, &delta_x, &delta_x_scaled, &delta_f, &error_f, &scaledError_f, n, x, fvec); + + /* updating x */ + memcpy(x, solverData->x_new, *n*sizeof(double)); + + /* updating f_old */ + memcpy(solverData->f_old, fvec, *n*sizeof(double)); + + current_fvec_enorm = error_f; + + /* check if maximum iteration is reached */ + if (++l > *maxfev) + { + *info = -1; + warningStreamPrint(LOG_NLS_V, 0, "Warning: maximal number of iteration reached but no root found"); + break; + } + } + + if(ACTIVE_STREAM(LOG_NLS_V)) + { + infoStreamPrint(LOG_NLS_V,1,"x vector"); + for(i=0; i<*n; i++) + infoStreamPrint(LOG_NLS_V, 0, "x[%d] = %e ", i, x[i]); + messageClose(LOG_NLS_V); + printErrors(delta_x, delta_x_scaled, delta_f, error_f, scaledError_f, eps); + } + } + + solverData->numberOfIterations += l; + solverData->numberOfFunctionEvaluations += solverData->nfev; + + return 0; +} + + +/*! \fn printErrors + * + * function prints errors, that reached tolerance + */ + +void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps) +{ + infoStreamPrint(LOG_NLS_V, 1, "errors "); + infoStreamPrint(LOG_NLS_V, 0, "delta_x = %e \ndelta_x_scaled = %e \ndelta_f = %e \nerror_f = %e \nscaledError_f = %e", delta_x, delta_x_scaled, delta_f, error_f, scaledError_f); + + if (delta_x < *eps) + infoStreamPrint(LOG_NLS_V, 0, "delta_x reached eps"); + if (delta_x_scaled < *eps) + infoStreamPrint(LOG_NLS_V, 0, "delta_x_scaled reached eps"); + if (delta_f < *eps) + infoStreamPrint(LOG_NLS_V, 0, "delta_f reached eps"); + if (error_f < *eps) + infoStreamPrint(LOG_NLS_V, 0, "error_f reached eps"); + if (scaledError_f < *eps) + infoStreamPrint(LOG_NLS_V, 0, "scaledError_f reached eps"); + + messageClose(LOG_NLS_V); +} + +/*! \fn solveLinearSystem + * + * function solves linear system J*(x_{n+1} - x_n) = f using lapack + */ +int solveLinearSystem(int* n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData) +{ + int i, nrsh=1, lapackinfo; + char trans = 'N'; + + /* if no factorization is given, calculate it */ + if (solverData->factorization == 0) + { + /* solve J*(x_{n+1} - x_n)=f */ + dgetrf_(n, n, fjac, n, iwork, &lapackinfo); + solverData->factorization = 1; + dgetrs_(&trans, n, &nrsh, fjac, n, iwork, fvec, n, &lapackinfo); + } + else + { + dgetrs_(&trans, n, &nrsh, fjac, n, iwork, fvec, n, &lapackinfo); + } + + if(lapackinfo > 0) + { + warningStreamPrint(LOG_NLS, 0, "Jacobian Matrix singular!"); + return -1; + } + else if(lapackinfo < 0) + { + warningStreamPrint(LOG_NLS, 0, "illegal input in argument %d", (int)lapackinfo); + return -1; + } + else + { + /* save solution of J*(x_{n+1} - x_n)=f */ + memcpy(solverData->x_increment, fvec, *n*sizeof(double)); + } + + return 0; +} + +/*! \fn calculatingErrors + * + * function calculates the errors + */ +void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f, + double* scaledError_f, int* n, double* x, double* fvec) +{ + int i=0; + double scaling_factor; + + /* delta_x = || x_new-x_old || */ + for (i=0; i<*n; i++) + solverData->delta_x_vec[i] = x[i]-solverData->x_new[i]; + + *delta_x = enorm_(n,solverData->delta_x_vec); + + scaling_factor = enorm_(n,x); + if (scaling_factor > 1) + *delta_x_scaled = *delta_x * 1./ scaling_factor; + else + *delta_x_scaled = *delta_x; + + /* delta_f = || f_old - f_new || */ + for (i=0; i<*n; i++) + solverData->delta_f[i] = solverData->f_old[i]-fvec[i]; + + *delta_f=enorm_(n, solverData->delta_f); + + *error_f = enorm_(n,fvec); + + /* scaling residual vector */ + scaling_residual_vector(solverData); + + for (i=0; i<*n; i++) + solverData->fvecScaled[i]=fvec[i]/solverData->resScaling[i]; + *scaledError_f = enorm_(n,solverData->fvecScaled); +} + +/*! \fn calculatingErrors + * + * function scales the residual vector using the jacobian (heuristic) + */ +void scaling_residual_vector(DATA_NEWTON* solverData) +{ + int i,j,k; + for(i=0, k=0; in; i++) + { + solverData->resScaling[i] = 0.0; + for(j=0; jn; j++, ++k) + { + solverData->resScaling[i] = fmax(fabs(solverData->fjac[k]), solverData->resScaling[i]); + } + if(solverData->resScaling[i] <= 0.0){ + warningStreamPrint(LOG_NLS_V, 1, "Jacobian matrix is singular."); + solverData->resScaling[i] = 1e-16; + } + solverData->fvecScaled[i] = solverData->fvec[i] / solverData->resScaling[i]; + } +} + +/*! \fn damping_heuristic + * + * first damping heuristic: + * x_increment will be halved until the Euclidean norm of the residual function + * is smaller than the Euclidean norm of the current point + * + * treshold for damping = 0.01 + * compiler flag: -newton = damped + */ +void damping_heuristic(double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, double* lambda, int* k, DATA_NEWTON* solverData, void* userdata) +{ + int i,j=0; + double enorm_new, treshold = 1e-2; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + enorm_new=enorm_(n,fvec); + + if (enorm_new >= current_fvec_enorm) + infoStreamPrint(LOG_NLS_V, 1, "Start Damping: enorm_new : %e; current_fvec_enorm: %e ", enorm_new, current_fvec_enorm); + + while (enorm_new >= current_fvec_enorm) + { + j++; + + *lambda*=0.5; + + + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i]; + + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + enorm_new=enorm_(n,fvec); + + if (*lambda <= treshold) + { + warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda reached a threshold."); + + /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */ + if (*k >= 5) + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i]; + else + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-solverData->x_increment[i]; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + (*k)++; + + break; + } + } + + *lambda = 1; + + messageClose(LOG_NLS_V); +} + +/*! \fn damping_heuristic2 + * + * second (default) damping heuristic: + * x_increment will be multiplied by 3/4 until the Euclidean norm of the residual function + * is smaller than the Euclidean norm of the current point + * + * treshold for damping = 0.0001 + * compiler flag: -newton = damped2 + */ +void damping_heuristic2(double damping_parameter, double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata) +{ + int i,j=0; + double enorm_new, treshold = 1e-4, lambda=1; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + enorm_new=enorm_(n,fvec); + + if (enorm_new >= current_fvec_enorm) + infoStreamPrint(LOG_NLS_V, 1, "StartDamping: "); + + while (enorm_new >= current_fvec_enorm) + { + j++; + + lambda*=damping_parameter; + + infoStreamPrint(LOG_NLS_V, 0, "lambda = %e, k = %d", lambda, *k); + + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i]; + + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + enorm_new=enorm_(n,fvec); + + if (lambda <= treshold) + { + warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda reached a threshold."); + + /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */ + if (*k >= 5) + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i]; + else + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-solverData->x_increment[i]; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + (*k)++; + + break; + } + } + + messageClose(LOG_NLS_V); +} + +/*! \fn LineSearch + * + * third damping heuristic: + * Along the tangent 5 five points are selected. For every point the Euclidean norm of + * the residual function will be calculated and the minimum is chosen for the further iteration. + * + * compiler flag: -newton = damped_ls + */ +void LineSearch(double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata) +{ + int i,j; + double enorm_new, enorm_minimum=current_fvec_enorm, lambda_minimum=0; + double lambda[5]={1.25,1,0.75,0.5,0.25}; + + + for (j=0; j<5; j++) + { + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-lambda[j]*solverData->x_increment[i]; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + enorm_new=enorm_(n,fvec); + + /* searching minimal enorm */ + if (enorm_new < enorm_minimum) + { + enorm_minimum = enorm_new; + lambda_minimum = lambda[j]; + memcpy(solverData->fvec_minimum, fvec,*n*sizeof(double)); + } + } + + infoStreamPrint(LOG_NLS_V,0,"lambda_minimum = %e", lambda_minimum); + + if (lambda_minimum == 0) + { + warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda_minimum = 0 "); + + /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */ + if (*k >= 5) + { + lambda_minimum = 0.125; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + } + else + { + lambda_minimum = 1; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + } + + (*k)++; + } + else + { + /* save new function values */ + memcpy(fvec, solverData->fvec_minimum, *n*sizeof(double)); + } + + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-lambda_minimum*solverData->x_increment[i]; +} + +/*! \fn Backtracking + * + * forth damping heuristic: + * Calculate new function h:R^n->R ; h(x) = 1/2 * ||f(x)|| ^2 + * g(lambda) = h(x_old + lambda * x_increment) + * find minimum of g with golden ratio method + * tau = golden ratio + * + * compiler flag: -newton = damped_bt + */ +void Backtracking(double* x, int(*f)(int*, double*, double*, void*, int), + double current_fvec_enorm, int* n, double* fvec, DATA_NEWTON* solverData, void* userdata) +{ + int i,j; + double enorm_new, enorm_f, lambda, a1, b1, a, b, tau, g1, g2; + double tolerance = 1e-3; + + /* saving current function values in f_old */ + memcpy(solverData->f_old, fvec, *n*sizeof(double)); + + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-solverData->x_increment[i]; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + + + /* calculate new enorm */ + enorm_new = enorm_(n,fvec); + + /* Backtracking only if full newton step is useless */ + if (enorm_new >= current_fvec_enorm) + { + infoStreamPrint(LOG_NLS_V, 0, "Start Backtracking\n enorm_new= %f \t current_fvec_enorm=%f",enorm_new, current_fvec_enorm); + + /* h(x) = 1/2 * ||f(x)|| ^2 + * g(lambda) = h(x_old + lambda * x_increment) + * find minimum of g with golden ratio method + * tau = golden ratio + * */ + + a = 0; + b = 1; + tau = 0.61803398875; + + a1 = a + (1-tau)*(b-a); + /* g1 = g(a1) = h(x_old - a1 * x_increment) = 1/2 * ||f(x_old- a1 * x_increment)||^2 */ + solverData->x_new[i] = x[i]- a1 * solverData->x_increment[i]; + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + enorm_f= enorm_(n,fvec); + g1 = 0.5 * enorm_f * enorm_f; + + + b1 = a + tau * (b-a); + /* g2 = g(b1) = h(x_old - b1 * x_increment) = 1/2 * ||f(x_old- b1 * x_increment)||^2 */ + solverData->x_new[i] = x[i]- b1 * solverData->x_increment[i]; + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + enorm_f= enorm_(n,fvec); + g2 = 0.5 * enorm_f * enorm_f; + + while ( (b - a) > tolerance) + { + if (g1x_new[i] = x[i]- a1 * solverData->x_increment[i]; + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + enorm_f= enorm_(n,fvec); + g1 = 0.5 * enorm_f * enorm_f; + } + else + { + a = a1; + a1 = b1; + b1 = a + tau * (b-a); + g1 = g2; + + /* g2 = g(b1) = h(x_old - b1 * x_increment) = 1/2 * ||f(x_old- b1 * x_increment)||^2 */ + solverData->x_new[i] = x[i]- b1 * solverData->x_increment[i]; + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + enorm_f= enorm_(n,fvec); + g2 = 0.5 * enorm_f * enorm_f; + } + } + + + lambda = (a+b)/2; + + + /* print lambda */ + infoStreamPrint(LOG_NLS_V, 0, "Backtracking - lambda = %e", lambda); + + for (i=0; i<*n; i++) + solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i]; + + /* calculate new function values */ + (*f)(n, solverData->x_new, fvec, userdata, 1); + solverData->nfev++; + } +} + + +#ifdef __cplusplus +} +#endif diff --git a/SimulationRuntime/c/simulation/solver/newtonIteration.h b/SimulationRuntime/c/simulation/solver/newtonIteration.h new file mode 100644 index 00000000000..af3b29bb321 --- /dev/null +++ b/SimulationRuntime/c/simulation/solver/newtonIteration.h @@ -0,0 +1,89 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-CurrentYear, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE + * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the OSMC (Open Source Modelica Consortium) + * Public License (OSMC-PL) are obtained from OSMC, either from the above + * address, from the URLs: http://www.openmodelica.org or + * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica + * distribution. GNU version 3 is obtained from: + * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from: + * http://www.opensource.org/licenses/BSD-3-Clause. + * + * This program is distributed WITHOUT ANY WARRANTY; without even the implied + * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS + * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE + * CONDITIONS OF OSMC-PL. + * + */ + +/*! \file nonlinearSolverNewton.h + */ + +#ifndef _NEWTONITERATION_H_ +#define _NEWTONITERATION_H_ + +#include "simulation_data.h" +#include "nonlinearSolverNewton.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct DATA_NEWTON +{ + double* resScaling; + double* fvecScaled; + + int newtonStrategy; + + int n; + double* x; + double* fvec; + double xtol; + double ftol; + int nfev; + int maxfev; + int info; + double epsfcn; + double* fjac; + double* rwork; + int* iwork; + int calculate_jacobian; + int factorization; + int numberOfIterations; /* over the whole simulation time */ + int numberOfFunctionEvaluations; /* over the whole simulation time */ + + /* damped newton */ + double* x_new; + double* x_increment; + double* f_old; + double* fvec_minimum; + double* delta_f; + double* delta_x_vec; + + rtclock_t timeClock; + +} DATA_NEWTON; + + +int allocateNewtonData(int size, void** data); +int freeNewtonData(void** data); +int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solverData, void* userdata); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c b/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c index 7ba3a9152b7..b8a5e377fc1 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c @@ -41,6 +41,7 @@ #include "util/varinfo.h" #include "model_help.h" #include "meta/meta_modelica.h" +#include "util/write_csv.h" #include "nonlinearSystem.h" #include "nonlinearSolverHomotopy.h" @@ -1087,6 +1088,7 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x) int assert = 1; threadData_t *threadData = solverData->data->threadData; + NONLINEAR_SYSTEM_DATA* nonlinsys = &(solverData->data->simulationInfo.nonlinearSystemData[solverData->data->simulationInfo.currentNonlinearSystemIndex]); /* debug information */ debugString(LOG_NLS_V, "******************************************************"); @@ -1259,6 +1261,23 @@ static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x) countNegativeSteps += (error_f > 10*error_f_old); error_f_old = error_f; + if (solverData->data->simulationInfo.nlsCsvInfomation){ + print_csvLineIterStats(((struct csvStats*) nonlinsys->csvData)->iterStats, + nonlinsys->size, + nonlinsys->numberOfCall+1, + numberOfIterations, + solverData->dy0, + solverData->dxScaled, + solverData->f1, + solverData->fvecScaled, + delta_x, + delta_x_scaled, + error_f, + error_f_scaled, + lambda + ); + } + if ((error_f_scaled < 1e-30*error_f) || countNegativeSteps > 20) { debugInt(LOG_NLS_V,"UPS! Something happened, NegativeSteps = ", countNegativeSteps); diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c b/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c index 7965aa0eee0..28ff6fe3c0b 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c @@ -28,7 +28,7 @@ * */ -/*! \file nonlinear_solver.c +/*! \file nonlinearSolverNewton.c */ #ifdef __cplusplus @@ -46,58 +46,14 @@ extern "C" { #include "nonlinearSystem.h" #include "nonlinearSolverNewton.h" +#include "newtonIteration.h" + +#include "external_input.h" + extern double enorm_(int *n, double *x); +int wrapper_fvec_newton(int* n, double* x, double* fvec, void* userdata, int fj); -typedef struct DATA_NEWTON -{ - int initialized; /* 1 = initialized, else = 0*/ - double* resScaling; - double* fvecScaled; - - int n; - double* x; - double* fvec; - double xtol; - double ftol; - int nfev; - int maxfev; - int info; - double epsfcn; - double* fjac; - double* rwork; - int* iwork; - int numberOfIterations; /* over the whole simulation time */ - int numberOfFunctionEvaluations; /* over the whole simulation time */ - - /* damped newton */ - double* x_new; - double* x_increment; - double* f_old; - double* fvec_minimum; - double* delta_f; - double* delta_x_vec; - - -} DATA_NEWTON; - - -static int _omc_newton(int* n, double *x, double *fvec, double* eps, double* fdeps, int* maxfev, - int(*f)(int*, double*, double*, int*, void*, int), double* fjac, - double* rwork, int* iwork, int* info, void* userdata, int sysNumber); - - -int solveLinearSystem(int* n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData); -void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f, - double* scaledError_f, int* n, double* x, double* fvec); -void scaling_residual_vector(DATA_NEWTON* solverData); -void damping_heuristic(void* userdata, int sysNumber, double* x, int(*f)(int*, double*, double*, int*, void*, int), - double current_fvec_enorm, int* n, double* fvec, double* lambda, int* k); -void damping_heuristic2(double damping_parameter, void* userdata, int sysNumber, double* x, int(*f)(int*, double*, double*, int*, void*, int), - double current_fvec_enorm, int* n, double* fvec, int* k); -void LineSearch(void* userdata, int sysNumber, double* x, int(*f)(int*, double*, double*, int*, void*, int), - double current_fvec_enorm, int* n, double* fvec, int* k); -void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps); #ifdef __cplusplus extern "C" { @@ -109,75 +65,6 @@ extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doubler } #endif -/*! \fn allocateNewtonData - * allocate memory for nonlinear system solver - */ -int allocateNewtonData(int size, void** voiddata) -{ - DATA_NEWTON* data = (DATA_NEWTON*) malloc(sizeof(DATA_NEWTON)); - - *voiddata = (void*)data; - assertStreamPrint(NULL, 0 != data, "allocationNewtonData() failed!"); - - data->initialized = 0; - data->resScaling = (double*) malloc(size*sizeof(double)); - data->fvecScaled = (double*) malloc(size*sizeof(double)); - - data->n = size; - data->x = (double*) malloc(size*sizeof(double)); - data->fvec = (double*) calloc(size,sizeof(double)); - data->xtol = 1e-12; - data->ftol = 1e-12; - data->maxfev = size*100; - data->epsfcn = DBL_EPSILON; - data->fjac = (double*) malloc((size*size)*sizeof(double)); - - data->rwork = (double*) malloc((size)*sizeof(double)); - data->iwork = (int*) malloc(size*sizeof(int)); - - /* damped newton */ - data->x_new = (double*) malloc(size*sizeof(double)); - data->x_increment = (double*) malloc(size*sizeof(double)); - data->f_old = (double*) calloc(size,sizeof(double)); - data->fvec_minimum = (double*) calloc(size,sizeof(double)); - data->delta_f = (double*) calloc(size,sizeof(double)); - data->delta_x_vec = (double*) calloc(size,sizeof(double)); - - data->numberOfIterations = 0; - data->numberOfFunctionEvaluations = 0; - - - assertStreamPrint(NULL, 0 != *voiddata, "allocationNewtonData() voiddata failed!"); - return 0; -} - -/*! \fn freeNewtonData - * - * free memory for nonlinear solver newton - * - */ -int freeNewtonData(void **voiddata) -{ - DATA_NEWTON* data = (DATA_NEWTON*) *voiddata; - - free(data->resScaling); - free(data->fvecScaled); - free(data->x); - free(data->fvec); - free(data->fjac); - free(data->rwork); - free(data->iwork); - - /* damped newton */ - free(data->x_new); - free(data->x_increment); - free(data->f_old); - free(data->fvec_minimum); - free(data->delta_f); - free(data->delta_x_vec); - - return 0; -} /*! \fn getAnalyticalJacobian * @@ -205,7 +92,7 @@ int getAnalyticalJacobianNewton(DATA* data, double* jac, int sysNumber) if(data->simulationInfo.analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) data->simulationInfo.analyticJacobians[index].seedVars[ii] = 1; - ((systemData->analyticalJacobianColumn))(data); + systemData->analyticalJacobianColumn(data); for(j = 0; j < data->simulationInfo.analyticJacobians[index].sizeCols; j++) { @@ -234,19 +121,72 @@ int getAnalyticalJacobianNewton(DATA* data, double* jac, int sysNumber) } -/*! \fn wrapper_fvec_hybrd for the residual Function +/*! \fn wrapper_fvec_newton for the residual Function * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data) * - * + * fj decides whether the function values or the jacobian matrix shall be calculated + * fj = 1 ==> calculate function values + * fj = 0 ==> calculate jacobian matrix */ -static int wrapper_fvec_newton(int* n, double* x, double* f, int* iflag, void* data, int sysNumber) +int wrapper_fvec_newton(int* n, double* x, double* fvec, void* userdata, int fj) { - int currentSys = sysNumber; - /* NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo.nonlinearSystemData[currentSys]); */ - /* DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); */ + DATA_USER* uData = (DATA_USER*) userdata; + DATA* data = (DATA*)(uData->data); + int currentSys = ((DATA_USER*)userdata)->sysNumber; + NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo.nonlinearSystemData[currentSys]); + DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); + int flag = 1; + int *iflag=&flag; - (*((DATA*)data)->simulationInfo.nonlinearSystemData[currentSys].residualFunc)(data, x, f, iflag); - return 0; + if (fj) + { + (data->simulationInfo.nonlinearSystemData[currentSys].residualFunc)(data, x, fvec, iflag); + } + else + { + if(systemData->jacobianIndex != -1) + { + getAnalyticalJacobianNewton(data, solverData->fjac, currentSys); + } + else + { + double delta_h = sqrt(solverData->epsfcn); + double delta_hh; + double xsave; + + int i,j,l, linear=0; + linear = systemData->method; + + for(i = 0; i < *n; i++) + { + if(linear) + { + delta_hh = 1; + } + else + { + + delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(fvec[i])), delta_h); + delta_hh = ((fvec[i] >= 0) ? delta_hh : -delta_hh); + delta_hh = x[i] + delta_hh - x[i]; + } + xsave = x[i]; + x[i] += delta_hh; + delta_hh = 1. / delta_hh; + + wrapper_fvec_newton(n, x, solverData->rwork, userdata, 1); + solverData->nfev++; + + for(j = 0; j < *n; j++) + { + l = i * *n + j; + solverData->fjac[l] = (solverData->rwork[j] - fvec[j]) * delta_hh; + } + x[i] = xsave; + } + } + } + return *iflag; } /*! \fn solve non-linear system with newton method @@ -260,12 +200,7 @@ int solveNewton(DATA *data, int sysNumber) { NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo.nonlinearSystemData[sysNumber]); DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); - /* - * We are given the number of the non-linear system. - * We want to look it up among all equations. - */ - int eqSystemNumber = systemData->equationIndex; - + int eqSystemNumber = 0; int i; double xerror = -1, xerror_scaled = -1; int success = 0; @@ -276,26 +211,42 @@ int solveNewton(DATA *data, int sysNumber) int giveUp = 0; int retries = 0; int retries2 = 0; - int iflag = 1; int nonContinuousCase = 0; + modelica_boolean *relationsPreBackup = NULL; - modelica_boolean *relationsPreBackup = (modelica_boolean*) malloc(data->modelData.nRelations*sizeof(modelica_boolean)); + DATA_USER* userdata = (DATA_USER*)malloc(sizeof(DATA_USER)); + assert(userdata != NULL); + + userdata->data = (void*)data; + userdata->sysNumber = sysNumber; + + /* + * We are given the number of the non-linear system. + * We want to look it up among all equations. + */ + eqSystemNumber = systemData->equationIndex; + + local_tol = solverData->ftol; + + relationsPreBackup = (modelica_boolean*) malloc(data->modelData.nRelations*sizeof(modelica_boolean)); solverData->nfev = 0; + /* try to calculate jacobian only once at the beginning of the iteration */ + solverData->calculate_jacobian = 0; + /* debug output */ if(ACTIVE_STREAM(LOG_NLS_V)) { int indexes[2] = {1,eqSystemNumber}; infoStreamPrintWithEquationIndexes(LOG_NLS, 1, indexes, "Start solving Non-Linear System %d at time %g with Newton Solver", - eqSystemNumber, - data->localData[0]->timeValue); + eqSystemNumber, data->localData[0]->timeValue); for(i = 0; i < solverData->n; i++) { - infoStreamPrint(LOG_NLS_V, 1, "x[%d] = %.15e", i, systemData->nlsx[i]); - infoStreamPrint(LOG_NLS_V, 0, "scaling = %f +++ old = %e +++ extrapolated = %e", - systemData->nominal[i], systemData->nlsxOld[i], systemData->nlsxExtrapolation[i]); + infoStreamPrint(LOG_NLS_V, 1, "x[%d] = %.15e", i, data->simulationInfo.discreteCall ? systemData->nlsx[i] : systemData->nlsxExtrapolation[i]); + infoStreamPrint(LOG_NLS_V, 0, "nominal = %g +++ nlsx = %g +++ old = %g +++ extrapolated = %g", + systemData->nominal[i], systemData->nlsx[i], systemData->nlsxOld[i], systemData->nlsxExtrapolation[i]); messageClose(LOG_NLS_V); } messageClose(LOG_NLS_V); @@ -312,10 +263,8 @@ int solveNewton(DATA *data, int sysNumber) { giveUp = 1; - _omc_newton(&solverData->n, solverData->x, solverData->fvec, &local_tol, - &solverData->epsfcn, &solverData->maxfev, - wrapper_fvec_newton, solverData->fjac, solverData->rwork, - solverData->iwork, &solverData->info, data, sysNumber); + solverData->newtonStrategy = data->simulationInfo.newtonStrategy; + _omc_newton(wrapper_fvec_newton, solverData, (void*)userdata); /* check for proper inputs */ if(solverData->info == 0) @@ -348,7 +297,7 @@ int solveNewton(DATA *data, int sysNumber) /* take the solution */ memcpy(systemData->nlsx, solverData->x, solverData->n*(sizeof(double))); - /* Then try with old values (instead of extrapolating )*/ + /* Then try with old values (instead of extrapolating )*/ } else if(retries < 1) { @@ -358,7 +307,10 @@ int solveNewton(DATA *data, int sysNumber) giveUp = 0; nfunc_evals += solverData->nfev; infoStreamPrint(LOG_NLS, 0, " - iteration making no progress:\t try old values."); - /* try to vary the initial values */ + /* try to vary the initial values */ + + /* evaluate jacobian in every step now */ + solverData->calculate_jacobian = 1; } else if(retries < 2) { @@ -368,7 +320,7 @@ int solveNewton(DATA *data, int sysNumber) giveUp = 0; nfunc_evals += solverData->nfev; infoStreamPrint(LOG_NLS, 0, " - iteration making no progress:\t vary solution point by 1%%."); - /* try to vary the initial values */ + /* try to vary the initial values */ } else if(retries < 3) { @@ -425,7 +377,7 @@ int solveNewton(DATA *data, int sysNumber) } } if(ACTIVE_STREAM(LOG_NLS)) - messageClose(LOG_NLS); + messageClose(LOG_NLS); free(relationsPreBackup); @@ -435,497 +387,3 @@ int solveNewton(DATA *data, int sysNumber) return success; } - -/*! \fn fdjac - * - * function calculates a jacobian matrix by - * numerical method finite differences - */ -static int fdjac(int* n, int(*f)(int*, double*, double*, int*, void*, int), double *x, - double* fvec, double *fjac, double* eps, int* iflag, double* wa, - void* userdata, int sysNumber) -{ - double delta_h = sqrt(*eps); - double delta_hh; - double delta_hhh; - double xsave; - - int i,j,l; - - NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)userdata)->simulationInfo.nonlinearSystemData[sysNumber]); - - int linear = systemData->method; - - for(i = 0; i < *n; i++) { - if(linear){ - delta_hh = 1; - } else { - delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(fvec[i])), delta_h); - delta_hh = ((fvec[i] >= 0) ? delta_hh : -delta_hh); - delta_hh = x[i] + delta_hh - x[i]; - } - xsave = x[i]; - x[i] += delta_hh; - delta_hh = 1. / delta_hh; - f(n, x, wa, iflag, userdata, sysNumber); - - for(j = 0; j < *n; j++) { - l = i * *n + j; - fjac[l] = (wa[j] - fvec[j]) * delta_hh; - } - x[i] = xsave; - } - - return *iflag; -} - -/*! \fn solve system with Newton-Raphson - * - * \param [in] [n] size of equation - * [eps] tolerance for x - * [h] tolerance for f' - * [k] maximum number of iterations - * [work] work array size of (n*X) - * [f] user provided function - * [data] userdata - * [info] - * - */ -static int _omc_newton(int* n, double *x, double *fvec, double* eps, double* fdeps, int* maxfev, - int(*f)(int*, double*, double*, int*, void*, int), - double* fjac, double* work, int* iwork, int* info, void* userdata, int sysNumber) -{ - DATA* data = (DATA*) userdata; - NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)userdata)->simulationInfo.nonlinearSystemData[sysNumber]); - DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); - - int i, j, k = 0, l = 0, iflag, nrsh = 1; - double error_f = 1.0 + *eps, scaledError_f = 1.0 + *eps, delta_x = 1.0 + *eps, delta_f = 1.0 + *eps, delta_x_scaled = 1.0 + *eps, lambda = 1.0; - double current_fvec_enorm, enorm_new; - - if(ACTIVE_STREAM(LOG_NLS_V)) - { - infoStreamPrint(LOG_NLS_V, 1, "######### Start Newton maxfev: %d #########", (int)*maxfev); - - infoStreamPrint(LOG_NLS_V, 1, "x vector"); - for(i=0; i<*n; i++) - infoStreamPrint(LOG_NLS_V, 0, "x[%d]: %e ", i, x[i]); - messageClose(LOG_NLS_V); - } - - *info = 1; - - /* calculate the function values */ - (*f)(n, x, fvec, &iflag, userdata,sysNumber); - solverData->nfev++; - - /* save current fvec in f_old*/ - memcpy(solverData->f_old, fvec, *n*sizeof(double)); - - error_f = current_fvec_enorm = enorm_(n, fvec); - - while(error_f > *eps && scaledError_f > *eps && delta_x > *eps && delta_f > *eps && delta_x_scaled > *eps ) - { - if(ACTIVE_STREAM(LOG_NLS_V)) - { - infoStreamPrint(LOG_NLS_V, 0, "\n**** start Iteration: %d *****", (int) l); - - /* Debug output */ - infoStreamPrint(LOG_NLS_V, 1, "function values"); - for(i=0; i<*n; i++) - infoStreamPrint(LOG_NLS_V, 0, "fvec[%d]: %e ", i, fvec[i]); - messageClose(LOG_NLS_V); - } - - /* calculate jacobian */ - if(systemData->jacobianIndex != -1) - { - getAnalyticalJacobianNewton(userdata, fjac, sysNumber); - } - else - { - fdjac(n, f, x, fvec, fjac, fdeps, &iflag, work, userdata, sysNumber); - solverData->nfev=solverData->nfev+*n; - } - - /* debug output */ - if(ACTIVE_STREAM(LOG_NLS_JAC)) - { - char buffer[4096]; - - infoStreamPrint(LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", (int)*n, (int)*n); - for(i=0; in;i++) - { - buffer[0] = 0; - for(j=0; jn; j++) - sprintf(buffer, "%s%10g ", buffer, fjac[i*(*n)+j]); - infoStreamPrint(LOG_NLS_JAC, 0, "%s", buffer); - } - messageClose(LOG_NLS_JAC); - } - - if (solveLinearSystem(n, iwork, fvec, fjac, solverData) != 0) - { - *info=-1; - break; - } - else - { - - for (i =0; i<*n; i++) - solverData->x_new[i]=x[i]-solverData->x_increment[i]; - - - infoStreamPrint(LOG_NLS_V,1,"x_increment"); - for(i=0; i<*n; i++) - infoStreamPrint(LOG_NLS_V, 0, "x_increment[%d] = %e ", i, solverData->x_increment[i]); - messageClose(LOG_NLS_V); - - - if (data->simulationInfo.newtonStrategy == NEWTON_DAMPED){ - damping_heuristic(userdata, sysNumber, x, f, current_fvec_enorm, n, fvec, &lambda, &k); - } else if (data->simulationInfo.newtonStrategy == NEWTON_DAMPED2){ - damping_heuristic2(0.75, userdata, sysNumber, x, f, current_fvec_enorm, n, fvec, &k); - } else if (data->simulationInfo.newtonStrategy == NEWTON_DAMPED_LS){ - LineSearch(userdata, sysNumber, x, f, current_fvec_enorm, n, fvec, &k); - } else { - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - } - - calculatingErrors(solverData, &delta_x, &delta_x_scaled, &delta_f, &error_f, &scaledError_f, n, x, fvec); - - /* updating x */ - memcpy(x, solverData->x_new, *n*sizeof(double)); - - /* updating f_old */ - memcpy(solverData->f_old, fvec, *n*sizeof(double)); - - current_fvec_enorm = error_f; - - /* check if maximum iteration is reached */ - if (++l > *maxfev) - { - *info = -1; - warningStreamPrint(LOG_NLS_V, 0, "Warning: maximal number of iteration reached but no root found"); - break; - } - } - - if(ACTIVE_STREAM(LOG_NLS_V)) - { - infoStreamPrint(LOG_NLS_V,1,"x vector"); - for(i=0; i<*n; i++) - infoStreamPrint(LOG_NLS_V, 0, "x[%d] = %e ", i, x[i]); - messageClose(LOG_NLS_V); - printErrors(delta_x, delta_x_scaled, delta_f, error_f, scaledError_f, eps); - } - - } - - solverData->numberOfIterations += l; - solverData->numberOfFunctionEvaluations += solverData->nfev; - - return 0; -} - - -void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps) -{ - infoStreamPrint(LOG_NLS_V, 1, "errors "); - infoStreamPrint(LOG_NLS_V, 0, "delta_x = %e \ndelta_x_scaled = %e \ndelta_f = %e \nerror_f = %e \nscaledError_f = %e", delta_x, delta_x_scaled, delta_f, error_f, scaledError_f); - - if (delta_x < *eps) - infoStreamPrint(LOG_NLS_V, 0, "delta_x reached eps"); - if (delta_x_scaled < *eps) - infoStreamPrint(LOG_NLS_V, 0, "delta_x_scaled reached eps"); - if (delta_f < *eps) - infoStreamPrint(LOG_NLS_V, 0, "delta_f reached eps"); - if (error_f < *eps) - infoStreamPrint(LOG_NLS_V, 0, "error_f reached eps"); - if (scaledError_f < *eps) - infoStreamPrint(LOG_NLS_V, 0, "scaledError_f reached eps"); - - messageClose(LOG_NLS_V); -} - -int solveLinearSystem(int* n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData) -{ - int i, nrsh=1, lapackinfo=1; - - /* solve J*(x_{n+1} - x_n)=f */ - dgesv_(n, &nrsh, fjac, n, iwork, fvec, n, &lapackinfo); - - if(lapackinfo > 0) - { - warningStreamPrint(LOG_NLS, 0, "Jacobian Matrix singular!"); - return -1; - } - else if(lapackinfo < 0) - { - warningStreamPrint(LOG_NLS, 0, "illegal input in argument %d", (int)lapackinfo); - return -1; - } - else - { - /* save solution of J*(x_{n+1} - x_n)=f */ - memcpy(solverData->x_increment, fvec, *n*sizeof(double)); - } - - return 0; -} - -void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f, - double* scaledError_f, int* n, double* x, double* fvec) -{ - int i=0; - double scaling_factor; - - /* delta_x = || x_new-x_old || */ - for (i=0; i<*n; i++) - solverData->delta_x_vec[i] = x[i]-solverData->x_new[i]; - - *delta_x = enorm_(n,solverData->delta_x_vec); - - scaling_factor = enorm_(n,x); - if (scaling_factor > 1) - *delta_x_scaled = *delta_x * 1./ scaling_factor; - else - *delta_x_scaled = *delta_x; - - /* delta_f = || f_old - f_new || */ - for (i=0; i<*n; i++) - solverData->delta_f[i] = solverData->f_old[i]-fvec[i]; - - *delta_f=enorm_(n, solverData->delta_f); - - *error_f = enorm_(n,fvec); - - /* scaling residual vector */ - scaling_residual_vector(solverData); - - for (i=0; i<*n; i++) - solverData->fvecScaled[i]=fvec[i]/solverData->resScaling[i]; - *scaledError_f = enorm_(n,solverData->fvecScaled); - -} - -void scaling_residual_vector(DATA_NEWTON* solverData) -{ - int i,j,k; - for(i=0, k=0; in; i++) - { - solverData->resScaling[i] = 0.0; - for(j=0; jn; j++, ++k) - { - solverData->resScaling[i] = fmax(fabs(solverData->fjac[k]), solverData->resScaling[i]); - } - if(solverData->resScaling[i] <= 0.0){ - warningStreamPrint(LOG_NLS_V, 1, "Jacobian matrix is singular."); - solverData->resScaling[i] = 1e-16; - } - solverData->fvecScaled[i] = solverData->fvec[i] / solverData->resScaling[i]; - } -} - -void damping_heuristic(void* userdata, int sysNumber, double* x, int(*f)(int*, double*, double*, int*, void*, int), - double current_fvec_enorm, int* n, double* fvec, double* lambda, int* k) -{ - int i,j=0, iflag; - double enorm_new, treshold = 1e-2; - - - NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)userdata)->simulationInfo.nonlinearSystemData[sysNumber]); - DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); - - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - enorm_new=enorm_(n,fvec); - - if (enorm_new >= current_fvec_enorm) - infoStreamPrint(LOG_NLS_V, 1, "Start Damping: enorm_new : %e; current_fvec_enorm: %e ", enorm_new, current_fvec_enorm); - - while (enorm_new >= current_fvec_enorm) - { - j++; - - *lambda*=0.5; - - - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i]; - - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - enorm_new=enorm_(n,fvec); - - if (*lambda <= treshold) - { - warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda reached a threshold."); - - /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */ - if (*k >= 5) - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i]; - else - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-solverData->x_increment[i]; - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - (*k)++; - - break; - } - } - - *lambda = 1; - - - messageClose(LOG_NLS_V); -} - -void damping_heuristic2(double damping_parameter, void* userdata, int sysNumber, double* x, int(*f)(int*, double*, double*, int*, void*, int), - double current_fvec_enorm, int* n, double* fvec, int* k) -{ - int i,j=0, iflag; - double enorm_new, treshold = 1e-4, lambda=1; - - - NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)userdata)->simulationInfo.nonlinearSystemData[sysNumber]); - DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); - - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - enorm_new=enorm_(n,fvec); - - if (enorm_new >= current_fvec_enorm) - infoStreamPrint(LOG_NLS_V, 1, "StartDamping: "); - - while (enorm_new >= current_fvec_enorm) - { - j++; - - lambda*=damping_parameter; - - infoStreamPrint(LOG_NLS_V, 0, "lambda = %e, k = %d", lambda, *k); - - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i]; - - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - enorm_new=enorm_(n,fvec); - - if (lambda <= treshold) - { - warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda reached a threshold."); - - /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */ - if (*k >= 5) - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i]; - else - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-solverData->x_increment[i]; - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - (*k)++; - - break; - } - } - - messageClose(LOG_NLS_V); -} - -void LineSearch(void* userdata, int sysNumber, double* x, int(*f)(int*, double*, double*, int*, void*, int), - double current_fvec_enorm, int* n, double* fvec, int* k) -{ - int i,j, iflag; - double enorm_new, treshold = 1e-2, enorm_minimum=current_fvec_enorm, lambda_minimum=0; - double lambda[5]={1.25,1,0.75,0.5,0.25}; - - - NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)userdata)->simulationInfo.nonlinearSystemData[sysNumber]); - DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); - - for (j=0; j<5; j++) - { - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-lambda[j]*solverData->x_increment[i]; - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - - enorm_new=enorm_(n,fvec); - - /* searching minimal enorm */ - if (enorm_new < enorm_minimum) - { - enorm_minimum = enorm_new; - lambda_minimum = lambda[j]; - memcpy(solverData->fvec_minimum, fvec,*n*sizeof(double)); - } - } - - infoStreamPrint(LOG_NLS_V,0,"lambda_minimum = %e", lambda_minimum); - - if (lambda_minimum == 0) - { - warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda_minimum = 0 "); - - /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */ - if (*k >= 5) - { - lambda_minimum = 0.125; - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - } - else - { - lambda_minimum = 1; - - /* calculate new function values */ - (*f)(n,solverData->x_new,fvec,&iflag,userdata,sysNumber); - solverData->nfev++; - } - - (*k)++; - } - else - { - /* save new function values */ - memcpy(fvec, solverData->fvec_minimum, *n*sizeof(double)); - } - - for (i=0; i<*n; i++) - solverData->x_new[i]=x[i]-lambda_minimum*solverData->x_increment[i]; - -} - - -#ifdef __cplusplus -} -#endif diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.h b/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.h index e7c3577c917..e2d7733c422 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.h +++ b/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.h @@ -40,8 +40,11 @@ extern "C" { #endif -int allocateNewtonData(int size, void** data); -int freeNewtonData(void** data); +typedef struct +{ + void* data; + int sysNumber; +}DATA_USER; int solveNewton(DATA *data, int sysNumber); diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSystem.c b/SimulationRuntime/c/simulation/solver/nonlinearSystem.c index 44cd7a7f3d7..ebdabe4092b 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSystem.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSystem.c @@ -39,20 +39,280 @@ #include "kinsolSolver.h" #include "nonlinearSolverHybrd.h" #include "nonlinearSolverNewton.h" +#include "newtonIteration.h" #include "nonlinearSolverHomotopy.h" #include "simulation/simulation_info_xml.h" #include "simulation/simulation_runtime.h" +#include "util/write_csv.h" /* for try and catch simulationJumpBuffer */ #include "meta/meta_modelica.h" int check_nonlinear_solution(DATA *data, int printFailingSystems, int sysNumber); -struct dataNewtonAndHybrid { +struct dataNewtonAndHybrid +{ void* newtonData; void* hybridData; }; +/*! \fn int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData) + * + * This function initializes csv files for analysis propose. + * + * \param [ref] [data] + * \param [ref] [systemData] + */ +int initializeNLScsvData(DATA* data, NONLINEAR_SYSTEM_DATA* systemData) +{ + struct csvStats* stats = (struct csvStats*) malloc(sizeof(struct csvStats)); + char buffer[100]; + sprintf(buffer, "%s_NLS%dStatsCall.csv", data->modelData.modelFilePrefix, (int)systemData->equationIndex); + stats->callStats = omc_write_csv_init(buffer, ',', '"'); + + sprintf(buffer, "%s_NLS%dStatsIter.csv", data->modelData.modelFilePrefix, (int)systemData->equationIndex); + stats->iterStats = omc_write_csv_init(buffer, ',', '"'); + + systemData->csvData = stats; + + return 0; +} + +/*! \fn int print_csvLineCallStatsHeader(OMC_WRITE_CSV* csvData) + * + * This function initializes csv files for analysis propose. + * + * \param [ref] [data] + * \param [ref] [systemData] + */ +int print_csvLineCallStatsHeader(OMC_WRITE_CSV* csvData) +{ + unsigned char buffer[1024] = ""; + + /* number of call */ + sprintf(buffer,"numberOfCall"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* simulation time */ + sprintf(buffer,"simulationTime"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving iterations */ + sprintf(buffer,"iterations"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving fCalls */ + sprintf(buffer,"numberOfFunctionCall"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving Time */ + sprintf(buffer,"solvingTime"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving Time */ + sprintf(buffer,"solvedSystem"); + omc_write_csv(csvData, buffer); + + /* finish line */ + fputc('\n',csvData->handle); + + return 0; +} + +/*! \fn int print_csvLineCallStats(OMC_WRITE_CSV* csvData) + * + * This function initializes csv files for analysis propose. + * + * \param [ref] [csvData] + * \param [in] [number of calls] + * \param [in] [simulation time] + * \param [in] [iterations] + * \param [in] [number of function call] + * \param [in] [solving time] + * \param [in] [solved system] + */ +int print_csvLineCallStats(OMC_WRITE_CSV* csvData, int num, double time, + int iterations, int fCalls, double solvingTime, + int solved) +{ + unsigned char buffer[1024]; + + /* number of call */ + sprintf(buffer, "%d", num); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* simulation time */ + sprintf(buffer, "%g", time); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving iterations */ + sprintf(buffer, "%d", iterations); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving fCalls */ + sprintf(buffer, "%d", fCalls); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving Time */ + sprintf(buffer, "%f", solvingTime); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solved system */ + sprintf(buffer, "%s", solved?"TRUE":"FALSE"); + omc_write_csv(csvData, buffer); + + /* finish line */ + fputc('\n',csvData->handle); + + return 0; +} + +/*! \fn int print_csvLineIterStatsHeader(OMC_WRITE_CSV* csvData) + * + * This function initializes csv files for analysis propose. + * + * \param [ref] [data] + * \param [ref] [systemData] + */ +int print_csvLineIterStatsHeader(DATA* data, NONLINEAR_SYSTEM_DATA* systemData, OMC_WRITE_CSV* csvData) +{ + unsigned char buffer[1024]; + int j; + int size = modelInfoGetEquation(&data->modelData.modelDataXml, systemData->equationIndex).numVar; + + /* number of call */ + sprintf(buffer,"numberOfCall"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* solving iterations */ + sprintf(buffer,"iteration"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* variables x */ + for(j=0; jmodelData.modelDataXml, systemData->equationIndex).vars[j]); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + } + + /* residuals */ + for(j=0; jseperator,csvData->handle); + } + + /* delta x */ + sprintf(buffer,"delta_x"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* delta x scaled */ + sprintf(buffer,"delta_x_scaled"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* error in f */ + sprintf(buffer,"error_f"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* error in f scaled */ + sprintf(buffer,"error_f_scaled"); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* damping lambda */ + sprintf(buffer,"lambda"); + omc_write_csv(csvData, buffer); + + /* finish line */ + fputc('\n',csvData->handle); + + return 0; +} + +/*! \fn int print_csvLineIterStatsHeader(OMC_WRITE_CSV* csvData) + * + * This function initializes csv files for analysis propose. + * + * \param [ref] [csvData] + * \param [in] [size, num, ...] + */ +int print_csvLineIterStats(OMC_WRITE_CSV* csvData, int size, int num, + int iteration, double* x, double* f, double error_f, + double error_fs, double delta_x, double delta_xs, + double lambda) +{ + unsigned char buffer[1024]; + int j; + + /* number of call */ + sprintf(buffer, "%d", num); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* simulation time */ + sprintf(buffer, "%d", iteration); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* x */ + for(j=0; jseperator,csvData->handle); + } + + /* r */ + for(j=0; jseperator,csvData->handle); + } + + /* error_f */ + sprintf(buffer, "%g", error_f); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* error_f */ + sprintf(buffer, "%g", error_fs); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* delta_x */ + sprintf(buffer, "%g", delta_x); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* delta_xs */ + sprintf(buffer, "%g", delta_xs); + omc_write_csv(csvData, buffer); + fputc(csvData->seperator,csvData->handle); + + /* lambda */ + sprintf(buffer, "%g", lambda); + omc_write_csv(csvData, buffer); + + /* finish line */ + fputc('\n',csvData->handle); + + return 0; +} + /*! \fn int initializeNonlinearSystems(DATA *data) * * This function allocates memory for all nonlinear systems. @@ -79,7 +339,8 @@ int initializeNonlinearSystems(DATA *data) assertStreamPrint(data->threadData, 0 != nonlinsys[i].residualFunc, "residual function pointer is invalid" ); /* check if analytical jacobian is created */ - if(nonlinsys[i].jacobianIndex != -1){ + if(nonlinsys[i].jacobianIndex != -1) + { assertStreamPrint(data->threadData, 0 != nonlinsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" ); if(nonlinsys[i].initialAnalyticalJacobian(data)) { @@ -98,6 +359,20 @@ int initializeNonlinearSystems(DATA *data) nonlinsys[i].initializeStaticNLSData(data, &nonlinsys[i]); + /* csv data call stats*/ + if (data->simulationInfo.nlsCsvInfomation) + { + if (initializeNLScsvData(data, &nonlinsys[i])) + { + throwStreamPrint(data->threadData, "csvData initialization failed"); + } + else + { + print_csvLineCallStatsHeader(((struct csvStats*) nonlinsys[i].csvData)->callStats); + print_csvLineIterStatsHeader(data, &nonlinsys[i], ((struct csvStats*) nonlinsys[i].csvData)->iterStats); + } + } + /* allocate solver data */ switch(data->simulationInfo.nlsMethod) { @@ -175,6 +450,7 @@ int freeNonlinearSystems(DATA *data) TRACE_PUSH int i; NONLINEAR_SYSTEM_DATA* nonlinsys = data->simulationInfo.nonlinearSystemData; + struct csvStats* stats; infoStreamPrint(LOG_NLS, 1, "free non-linear system solvers"); @@ -187,6 +463,13 @@ int freeNonlinearSystems(DATA *data) free(nonlinsys[i].min); free(nonlinsys[i].max); + if (data->simulationInfo.nlsCsvInfomation) + { + stats = nonlinsys[i].csvData; + omc_write_csv_free(stats->callStats); + omc_write_csv_free(stats->iterStats); + } + /* free solver data */ switch(data->simulationInfo.nlsMethod) { @@ -265,7 +548,8 @@ int solve_nonlinear_system(DATA *data, int sysNumber) rt_ext_tp_tick(&nonlinsys->totalTimeClock); - if(data->simulationInfo.discreteCall){ + if(data->simulationInfo.discreteCall) + { double *fvec = malloc(sizeof(double)*nonlinsys->size); int success = 0; @@ -283,7 +567,8 @@ int solve_nonlinear_system(DATA *data, int sysNumber) /*catch */ MMC_CATCH_INTERNAL(simulationJumpBuffer) #endif - if (!success) { + if (!success) + { warningStreamPrint(LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert."); } @@ -341,8 +626,11 @@ int solve_nonlinear_system(DATA *data, int sysNumber) default: throwStreamPrint(data->threadData, "unrecognized nonlinear solver"); } + + nonlinsys->solved = success; + #ifndef OMC_EMCC /*catch */ MMC_CATCH_INTERNAL(simulationJumpBuffer) @@ -355,6 +643,18 @@ int solve_nonlinear_system(DATA *data, int sysNumber) nonlinsys->totalTime += rt_ext_tp_tock(&(nonlinsys->totalTimeClock)); nonlinsys->numberOfCall++; + if (data->simulationInfo.nlsCsvInfomation) + { + print_csvLineCallStats(((struct csvStats*) nonlinsys->csvData)->callStats, + nonlinsys->numberOfCall, + data->localData[0]->timeValue, + nonlinsys->numberOfIterations, + nonlinsys->numberOfFEval, + nonlinsys->totalTime, + nonlinsys->solved + ); + } + return check_nonlinear_solution(data, 1, sysNumber); } @@ -437,7 +737,7 @@ int check_nonlinear_solution(DATA *data, int printFailingSystems, int sysNumber) /*! \fn extraPolate * This function extrapolates linear next value from - * the both old values, + * the both old values. * * \param [in] [data] * diff --git a/SimulationRuntime/c/simulation_data.h b/SimulationRuntime/c/simulation_data.h index 5584a9a6f85..0d83d3c6d2c 100644 --- a/SimulationRuntime/c/simulation_data.h +++ b/SimulationRuntime/c/simulation_data.h @@ -313,6 +313,8 @@ typedef struct NONLINEAR_SYSTEM_DATA unsigned long numberOfIterations; /* number of iteration of non-linear solvers of this system */ double totalTime; /* save the totalTime */ rtclock_t totalTimeClock; /* time clock for the totalTime */ + + void* csvData; /* information to save csv data */ }NONLINEAR_SYSTEM_DATA; typedef struct LINEAR_SYSTEM_DATA @@ -491,6 +493,7 @@ typedef struct SIMULATION_INFO int mixedMethod; /* mixed solver */ int nlsMethod; /* nonlinear solver */ int newtonStrategy; /* newton damping strategy solver */ + int nlsCsvInfomation; /* = 1 csv files with detailed nonlinear solver process are generated */ double lambda; /* homotopy parameter E [0, 1.0] */ diff --git a/SimulationRuntime/c/util/simulation_options.c b/SimulationRuntime/c/util/simulation_options.c index d3261458f91..2be92a0f4ca 100644 --- a/SimulationRuntime/c/util/simulation_options.c +++ b/SimulationRuntime/c/util/simulation_options.c @@ -65,6 +65,7 @@ const char *FLAG_NAME[FLAG_MAX+1] = { /* FLAG_MEASURETIMEPLOTFORMAT */ "measureTimePlotFormat", /* FLAG_NEWTON_STRATEGY */ "newton", /* FLAG_NLS */ "nls", + /* FLAG_NLS_INFO */ "nlsInfo", /* FLAG_NOEMIT */ "noemit", /* FLAG_NOEQUIDISTANT_GRID */ "noEquidistantTimeGrid", /* FLAG_NOEQUIDISTANT_OUT_FREQ*/ "noEquidistantOutputFrequency", @@ -120,6 +121,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = { /* FLAG_MEASURETIMEPLOTFORMAT */ "value specifies the output format of the measure time functionality", /* FLAG_NEWTON_STRATEGY */ "value specifies the damping strategy for the newton solver", /* FLAG_NLS */ "value specifies the nonlinear solver", + /* FLAG_NLS_INFO */ "outputs detailed information about solving process of non-linear systems into csv files.", /* FLAG_NOEMIT */ "do not emit any results to the result file", /* FLAG_NOEQUIDISTANT_GRID */ "stores results not in equidistant time grid as given by stepSize or numberOfIntervals, instead the variable step size of dassl is used.", /* FLAG_NOEQUIDISTANT_OUT_FREQ*/ "value controls the output frequency in noEquidistantTimeGrid mode", @@ -229,6 +231,8 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = { " * kinsol\n" " * newton\n" " * mixed", + /* FLAG_NLS_INFO */ + " Outputs detailed information about solving process of non-linear systems into csv files.", /* FLAG_NOEMIT */ " Do not emit any results to the result file.", /* FLAG_NOEQUIDISTANT_GRID */ @@ -313,6 +317,7 @@ const int FLAG_TYPE[FLAG_MAX] = { /* FLAG_MEASURETIMEPLOTFORMAT */ FLAG_TYPE_OPTION, /* FLAG_NEWTON_STRATEGY */ FLAG_TYPE_OPTION, /* FLAG_NLS */ FLAG_TYPE_OPTION, + /* FLAG_NLS_INFO */ FLAG_TYPE_FLAG, /* FLAG_NOEMIT */ FLAG_TYPE_FLAG, /* FLAG_NOEQUIDISTANT_GRID*/ FLAG_TYPE_FLAG, /* FLAG_NOEQUIDISTANT_OUT_FREQ*/ FLAG_TYPE_OPTION, @@ -440,6 +445,7 @@ const char *NEWTONSTRATEGY_NAME[NEWTON_MAX+1] = { /* NEWTON_DAMPED */ "damped", /* NEWTON_DAMPED2 */ "damped2", /* NEWTON_DAMPED_LS */ "damped_ls", + /* NEWTON_DAMPED_BT */ "damped_bt", /* NEWTON_PURE */ "pure", "NEWTON_MAX" @@ -451,6 +457,7 @@ const char *NEWTONSTRATEGY_DESC[NEWTON_MAX+1] = { /* NEWTON_DAMPED */ "Newton with a damping strategy", /* NEWTON_DAMPED2 */ "Newton with a damping strategy 2", /* NEWTON_DAMPED_LS */ "Newton with a damping line search", + /* NEWTON_DAMPED_BT */ "Newton with a damping backtracking and a minimum search via golden ratio method", /* NEWTON_PURE */ "Newton without damping strategy", "NEWTON_MAX" diff --git a/SimulationRuntime/c/util/simulation_options.h b/SimulationRuntime/c/util/simulation_options.h index 803d02fa6ce..7aad9c023e4 100644 --- a/SimulationRuntime/c/util/simulation_options.h +++ b/SimulationRuntime/c/util/simulation_options.h @@ -73,6 +73,7 @@ enum _FLAG FLAG_MEASURETIMEPLOTFORMAT, FLAG_NEWTON_STRATEGY, FLAG_NLS, + FLAG_NLS_INFO, FLAG_NOEMIT, FLAG_NOEQUIDISTANT_GRID, FLAG_NOEQUIDISTANT_OUT_FREQ, @@ -183,6 +184,7 @@ enum NEWTON_STRATEGY NEWTON_DAMPED, NEWTON_DAMPED2, NEWTON_DAMPED_LS, + NEWTON_DAMPED_BT, NEWTON_PURE, NEWTON_MAX diff --git a/SimulationRuntime/c/util/write_csv.c b/SimulationRuntime/c/util/write_csv.c new file mode 100644 index 00000000000..0f4e4ffdcaa --- /dev/null +++ b/SimulationRuntime/c/util/write_csv.c @@ -0,0 +1,83 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-CurrentYear, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3 LICENSE OR + * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the Open Source Modelica + * Consortium (OSMC) Public License (OSMC-PL) are obtained + * from OSMC, either from the above address, + * from the URLs: http://www.ida.liu.se/projects/OpenModelica or + * http://www.openmodelica.org, and in the OpenModelica distribution. + * GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html. + * + * This program is distributed WITHOUT ANY WARRANTY; without + * even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH + * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL. + * + * See the full OSMC Public License conditions for more details. + * + */ + +#include +#include + +#include "libcsv.h" +#include "write_csv.h" + +#define CSV_BUFFER_SIZE 1024 + +OMC_WRITE_CSV* omc_write_csv_init(char filename[], char seperator, char quote){ + int i; + size_t n = strlen(filename); + + OMC_WRITE_CSV* csvData = (OMC_WRITE_CSV*) malloc(sizeof(OMC_WRITE_CSV)); + csvData->filename = (char*) malloc((n+1)*sizeof(char)); + + strncpy(csvData->filename, filename, n); + csvData->filename[n] = '\0'; + csvData->seperator = seperator; + csvData->quote = '"'; + + csvData->handle = fopen(csvData->filename,"w"); + + return csvData; +} + + +int omc_write_csv_free(OMC_WRITE_CSV* csvData){ + + int i; + + free(csvData->filename); + fclose(csvData->handle); + + return 0; +} + +int omc_write_csv(OMC_WRITE_CSV* csvData, const void* csvLine){ + + size_t dest_size; + unsigned char buffer[CSV_BUFFER_SIZE] = ""; + + dest_size = csv_write(&buffer, CSV_BUFFER_SIZE, csvLine, strlen(csvLine)); + if (dest_size > CSV_BUFFER_SIZE){ + unsigned char* newbuffer = (unsigned char*) malloc(dest_size*sizeof(char)); + dest_size = csv_write(&newbuffer, dest_size, csvLine, strlen(csvLine)); + fprintf(csvData->handle, "%s", newbuffer); + }else{ + fprintf(csvData->handle, "%s", buffer); + } + + return 0; +} diff --git a/SimulationRuntime/c/util/write_csv.h b/SimulationRuntime/c/util/write_csv.h new file mode 100644 index 00000000000..ed40edf70af --- /dev/null +++ b/SimulationRuntime/c/util/write_csv.h @@ -0,0 +1,63 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-CurrentYear, Open Source Modelica Consortium (OSMC), + * c/o Linköpings universitet, Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3 LICENSE OR + * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2. + * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES + * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3, + * ACCORDING TO RECIPIENTS CHOICE. + * + * The OpenModelica software and the Open Source Modelica + * Consortium (OSMC) Public License (OSMC-PL) are obtained + * from OSMC, either from the above address, + * from the URLs: http://www.ida.liu.se/projects/OpenModelica or + * http://www.openmodelica.org, and in the OpenModelica distribution. + * GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html. + * + * This program is distributed WITHOUT ANY WARRANTY; without + * even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH + * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL. + * + * See the full OSMC Public License conditions for more details. + * + */ + +#ifndef OMC_WRITE_CSV_H +#define OMC_WRITE_CSV_H + +#include +#include + +typedef struct omc_write_csv { + char* filename; + FILE* handle; + char seperator; + char quote; +} OMC_WRITE_CSV; + +struct csvStats { + OMC_WRITE_CSV* callStats; + OMC_WRITE_CSV* iterStats; +}; + +#ifdef __cplusplus +extern "C" { +#endif + +OMC_WRITE_CSV* +omc_write_csv_init(char filename[], char seperator, char quote); +int omc_write_csv_free(OMC_WRITE_CSV* csvData); +int omc_write_csv(OMC_WRITE_CSV* csvData, const void* csvLine); + +#ifdef __cplusplus +} /* extern "C" */ +#endif + +#endif