Fixed the stuff removed in luc's merge.... :(

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@1732 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

modeq/modpar/Codegen.cpp

@@ -1,4 +1,7 @@
 #include "Codegen.hpp"
+#include "boost/tokenizer.hpp"
+using namespace std;
+using namespace boost;
 
 const string TAB = string("        ");
 
@@ -344,6 +347,8 @@ void Codegen::generateSubTaskCode(VertexID task)
     break;
   case NonLinSys:
     m_cstream << TAB << "/* NonLinear system*/" << endl;
+    generateNonLinearResidualFunc(task);
+    m_cstream << TAB << "/* end NonLinear system*/" << endl;
     break;
   case Copy:
     m_cstream << TAB << getResultName(task,m_tg) 
@@ -356,6 +361,59 @@ void Codegen::generateSubTaskCode(VertexID task)
   };
 }
 
+  // Generates the residual function for solving nonlinear equation systems.
+void Codegen::generateNonLinearResidualFunc(VertexID task) 
+{
+  int tasknumber = getTaskID(task,m_tg); // Task no. used to make unique funcion name
+  m_cstreamFunc << "void residualFunc" << tasknumber 
+		<< "( int n, double *xloc, double *res, int iflag) {" << endl;
+  m_cstreamFunc << getVertexName(task,m_tg) << endl;
+  m_cstreamFunc << "}" << endl;
+
+  int n = out_degree(task,*m_tg);
+  int lr = (n*(n+1))/2;
+  m_cstream << TAB <<"{ double nls_x[" << n <<"];" << endl;
+  m_cstream << TAB <<"double nls_diag[" << n <<"];" << endl;
+  m_cstream << TAB <<"double nls_qtf[" << n << "];" << endl;
+  m_cstream << TAB <<"double nls_wa1[" << n << "];" << endl;
+  m_cstream << TAB <<"double nls_wa2[" << n << "];" << endl;
+  m_cstream << TAB <<"double nls_wa3[" << n << "];" << endl;
+  m_cstream << TAB << "double nls_wa4[" << n << "];" << endl;
+  m_cstream << TAB << "int info,nfev;" << endl;
+  m_cstream << TAB << "double nls_fjac[" << n*n << "];" << endl;
+
+  InEdgeIterator e,e_end;		// Iterate over parents and set values 
+  int i;
+  for (i=0,tie(e,e_end) = in_edges(task,*m_tg); e != e_end; e++,i++) {
+    ResultSet &s = getResultSet(*e,m_tg);
+    m_cstream << TAB << "nls_x[" << i << "] = " << s.top() << ";\n" << endl;
+  }
+
+  m_cstream << TAB << "hybrd(residualFunc" << tasknumber << "," << n << ", &nls_x, 1e-6," 
+	    << "2000, " << n-1 << ", " << n-1 << ", 1e-6, &nls_diag, 1, 100.0,"
+	    << " -1, &info, &nfev, &nls_fjac, "<< n <<", " << lr
+	    <<", &nls_qtf, &nls_wa1, &nls_wa2, &nls_wa3, &nls_wa4);" << endl;
+  m_cstream << TAB << "if (info == 0) { printf(\"improper input parameters to nonlinear system nuber  " << tasknumber << "\");" << endl;
+  m_cstream << TAB << "exit(-3);" << endl;
+  m_cstream << TAB << " }"<< endl;
+
+  m_cstream << TAB << "if (info >= 2 && info <= 5) { printf(\"error solving nonlinear system number " << tasknumber << "\");" << endl;
+  m_cstream << TAB << "exit(-2);" << endl;
+  m_cstream << TAB << " }"<< endl;
+
+
+  string & results = getResultName(task,m_tg);
+  char_separator<char> sep(",");  
+  tokenizer<char_separator<char> >  tokens(results,sep);
+  i=0;
+  for (tokenizer<char_separator<char> >::iterator beg=tokens.begin(); beg != tokens.end();++beg,++i) {
+    m_cstream << TAB << *beg <<  " =  nls_x[" << i << "];\n" << endl;
+
+  }
+
+  m_cstream << "}" << endl; // Closing scope for nonlinear system.
+}
+
 void Codegen::generateParallelFunctionHeader(int procno) 
 {
   m_cstreamFunc << "void proc" << procno << "(";
@@ -473,6 +531,7 @@ void  Codegen::generateParallelMPIHeaders()
 void  Codegen::generateParallelMPIGlobals()
 {
   m_cstreamFunc  << "#include <mpi.h>" << endl;
+  m_cstreamFunc  << "#include <stdio.h>" << endl;
   m_cstreamFunc  << "/* Declaration of MPI Global variables */" << endl;
   m_cstreamFunc  << "extern MPI_Status status;" << endl;
   m_cstreamFunc  << "extern MPI_Request request;" << endl;

modeq/modpar/Codegen.hpp

@@ -62,6 +62,9 @@ class Codegen
   void generateTemporaries();
   void generateParallelFunctionLocals(TaskList *tasks);
 
+  // Generates the residual function for solving nonlinear equation systems.
+  void generateNonLinearResidualFunc(VertexID task);
+
   void generateKillCommand();
   void generateTmpDeclarations();
 

modeq/taskgraph.rml

@@ -136,6 +136,7 @@ with "dae.rml"
 with "ceval.rml"
 with "values.rml"
 with "print.rml"
+with "vartransform.rml"
 
 relation build_taskgraph: (DAELow.DAELow,  int vector, int vector, int list list) => () =
 
@@ -345,27 +346,214 @@ relation build_equation:(DAELow.DAELow, int vector, int vector, int) => () =
 	--------------------------------
 	build_equation(DAELow.DAELOW(vars,_,eqns,_,_,_),ass1,ass2,e)	
 
-  rule	int_sub(e,1) => e' &
+  rule	(* non-state nonlinear *)
+	int_sub(e,1) => e' &
 	DAELow.equation_nth(eqns,e') => DAELow.EQUATION(e1,e2) &
 	vector_nth(ass2,e') => v & (* v==variable no solved in this equation *)
 	int_sub(v,1) => v' &
-	DAELow.vararray_nth(vararr,v') => DAELow.VAR(cr,_,_,_,_,_,_,_,_,origname,_,dae_var_attr,comment) &  
+	DAELow.vararray_nth(vararr,v') => DAELow.VAR(cr,kind,_,_,_,_,_,_,_,origname,_,dae_var_attr,_) &  
+	is_non_state kind &
 	let varexp = Exp.CREF(cr,Exp.REAL) &
 	not Exp.solve(e1,e2,varexp) => _ &
-	print "nonlinear equation not implemented yet\n" 
+	build_nonlinear_equations([varexp],[Exp.BINARY(e1,Exp.SUB(Exp.REAL),e2)]) 
 	--------------------------------
-	build_equation(DAELow.DAELOW(DAELow.VARIABLES(_,vararr,_,_),_,eqns,_,_,_),ass1,ass2,e) => fail
+	build_equation(DAELow.DAELOW(DAELow.VARIABLES(_,vararr,_,_),_,eqns,_,_,_),ass1,ass2,e) => ()
+
+  rule	(* state nonlinear *)
+	int_sub(e,1) => e' &
+	DAELow.equation_nth(eqns,e') => DAELow.EQUATION(e1,e2) &
+	vector_nth(ass2,e') => v &
+	(* v == variable no solved in this equation *)
+	int_sub(v,1) => v' &
+	DAELow.var_list(vars) => varlst &
+	list_nth(varlst,v') => DAELow.VAR(cr,DAELow.STATE,_,_,_,_,_,_,indx,origname,_,dae_var_attr,_) &
+	int_string(indx) => indxs &
+	Util.string_append_list(["xd[",indxs,"]"]) => id &
+	let cr' = Exp.CREF_IDENT(id,[]) &
+	let varexp = Exp.CREF(cr',Exp.REAL) &
+	not Exp.solve(e1,e2,varexp) => _ &
+	build_nonlinear_equations([varexp],[Exp.BINARY(e1,Exp.SUB(Exp.REAL),e2)]) 
+	--------------------------------
+	build_equation(DAELow.DAELOW(vars,_,eqns,_,_,_),ass1,ass2,e)	
 
   rule	print "-build_equation failed\n" 
 	--------------------------------
 	build_equation(_,_,_,_) => fail
 end
 
+(** relation: build_nonlinear_equations
+ ** builds task graph for solving non-linear equations
+ **)
+
+relation build_nonlinear_equations: (Exp.Exp list, (* variables *)
+				     Exp.Exp list) (* residuals *)
+	=> () =
+
+  rule	list_length(vars) => size &
+	int_string(size) => size_str &
+	build_residual_code(vars,residuals) => taskname &
+	TaskGraphExt.new_task(taskname) => tid &
+	TaskGraphExt.set_tasktype(tid,3) & (* See TaskType in TaskGraph.hpp *)
+	build_nonlinear_equations2(tid,vars,residuals) &
+	Util.list_map(vars,Exp.print_exp_str) => varnames &
+	store_multiple_results(varnames,tid)
+	---------------------------------------------
+	build_nonlinear_equations(vars,residuals) => ()
+
+  rule	print "build_nonlinear_equatins failed\n" 
+	-----------------------------------------
+	build_nonlinear_equations(vars,residuals) => fail
+end
+
+(** relation: build_residual_code
+ ** This relation takes a list of expressions and builds code for 
+ ** calculating the residuals as a string. Used for e.g. solving non-linear equations.
+ **)
+relation build_residual_code:(Exp.Exp list, (* vars *)
+			      Exp.Exp list) (* residuals *)
+	  => string =
+
+  rule	make_residual_replacements(vars) => repl &
+	build_residual_code2(es,0,repl) => res
+	-------------------
+	build_residual_code(vars,es) => res
+
+  rule	print "build_residual_code failed\n" 
+	--------------------------
+	build_residual_code(_,_) => fail
+end
+
+(** relation: make_residual_replacements
+ ** This relation makes replacement rules for variables occuring in a 
+ ** nonlinear equation system. They should be replaced by x[index], i.e.
+ ** an unique index in the x vector.
+ **)
+relation make_residual_replacements:(Exp.Exp list) => 
+	VarTransform.VariableReplacements =
+
+  rule	VarTransform.empty_replacements() => repl &
+	make_residual_replacements2(repl,expl,0) => repl'
+	---------------------
+	make_residual_replacements(expl) => repl'
+end
+
+relation make_residual_replacements2:(VarTransform.VariableReplacements,
+				      Exp.Exp list,
+				      int) => 
+	VarTransform.VariableReplacements =
+  axiom	make_residual_replacements2(repl,[],_) => repl
+
+  rule	int_string(pos) => pstr &
+	Util.string_append_list(["xloc[",pstr,"]"]) => str &
+	VarTransform.add_replacement(repl,cr,Exp.CREF_IDENT(str,[])) => repl' &
+	pos + 1 => pos' &
+	make_residual_replacements2(repl',es,pos') => repl''
+	------------------------------
+	make_residual_replacements2(repl,Exp.CREF(cr,_)::es,pos) => repl''
+end
+
+
+relation build_residual_code2:(Exp.Exp list, 
+			       int, 
+			       VarTransform.VariableReplacements) 
+	  => string =
+
+  axiom	build_residual_code2([],_,_) => "" 
+
+  rule	VarTransform.replace_exp(e,repl) => e' &
+	Exp.print_exp_str(e') => s1 &
+	pos + 1 => pos'  &
+	build_residual_code2(es,pos',repl) => s2 &
+	int_string(pos) => pstr &
+	Util.string_append_list(["res[",pstr,"]=",s1,";\n",s2]) => res
+	-------------------------
+	build_residual_code2(e::es,pos,repl) => res
+
+  rule	print "build_residual_code2 failed\n" 
+	-------------------------
+	build_residual_code2(_,_,_) => fail
+end
+
+(** relation store_multiple_results
+ ** When a task calculates several values, this relation is used.
+ ** It collects the names of the values into one string, separated by semicolons
+ ** and uses that as the resultstring.
+ **)
+relation store_multiple_results: (string list, (* var names*)
+				  int (* task id*)
+				  )
+	=> () =
+
+  rule	Util.string_delimit_list(varnames,";") => result_str &
+	TaskGraphExt.store_result(result_str,tid,true,result_str)
+	-----------------------------------------
+	store_multiple_results(varnames,tid) => ()
+
+  rule	print "store_multiple_results failed\n" 
+	-----------------------
+	store_multiple_results(_,_) => fail
+end
+
+
+relation build_nonlinear_equations2: (int, (*task id *)
+				      Exp.Exp list, (* vars *)
+				      Exp.Exp list) (* residuals *)
+	=> () =
+
+  axiom	build_nonlinear_equations2(tid,_,[]) => ()
+
+  rule	(* Collect all variables and construct
+	 a string for the residual, that can be directly used in codegen.*)
+	Exp.get_cref_from_exp(res) => vars1 &
+	Util.list_map(vars,Exp.get_cref_from_exp) => vars' &
+	Util.list_flatten(vars') => vars2 &
+	(* No duplicate elements *)
+	Util.list_union_p(vars1,vars2,Exp.cref_equal) => vars1' &
+	Util.list_setdifference_p(vars1',vars2,Exp.cref_equal) => vars &
+	add_edges_from_vars(vars,tid,0) 
+	--------------------------------
+	build_nonlinear_equations2(tid,vars,res::residuals) => ()
+
+  rule	print "build_nonlinear_equations2 failed\n" &
+	Exp.print_exp_str e => es & print "first residual :" &
+	print es & print "\n" 
+	----------------------------
+	build_nonlinear_equations2(_,_,e::_) => fail
+end
+
+(** relation: add_edges_from_vars
+ ** Adds an edge between the tasks where the variables are defined and the tasks
+ ** given as second argument.
+ **)
+relation add_edges_from_vars: (Exp.ComponentRef list, 
+			       int (* task *),
+			       int (* priority*)) => () =
+
+  axiom	add_edges_from_vars([],_,_) => ()
+
+  rule	Exp.cref_str(v) => v_str &
+	TaskGraphExt.get_task(v_str) => predt &
+	TaskGraphExt.add_edge(predt,tid,v_str,prio) &
+	prio + 1 => prio' &
+	add_edges_from_vars(vs,tid,prio') 
+	------------------------------
+	add_edges_from_vars(v::vs,tid,prio)
+
+  rule	Exp.cref_str(v) => v_str &
+	not TaskGraphExt.get_task(v_str) => _ &
+	print "task " & print v_str & print " not found\n" 
+	------------------------------------------------
+	add_edges_from_vars(v::vs,_,_) => fail
+
+  rule	print "add_edges_from_vars failed\n" 
+	----------------------
+	add_edges_from_vars(_,_,_) => fail
+end
+
 (* Build task graph for a system of equations *)
 relation build_system:(DAELow.DAELow, int vector, int vector, int list) => () =
 
-  rule	print "build system\n" &
-	TaskGraphExt.new_task("equation system") => tid &
+  rule	TaskGraphExt.new_task("equation system") => tid &
 	build_system2(dae,ass1,ass2,system,tid) => predtasks &
 	Util.list_map(predtasks, TaskGraphExt.get_task) => predtaskids &
 	add_predecessors(tid,predtaskids,predtasks,0) 
@@ -540,6 +728,14 @@ relation build_expression:(Exp.Exp) => (int, string) =
 	build_call_str(funcstr,numargs) => ts &
 	Util.list_map__2(expl,build_expression) => (tasks,strs) &
 	TaskGraphExt.new_task(ts) => t &
+	print "building call " & print ts & 
+	int_string t => t_str & print " for task " & print  t_str & print "\n" &
+	Util.string_delimit_list(strs,", ") => s &
+	print " predecessor strs :" & print s & print "\n" &
+	Util.list_map(tasks,int_string) => tasks_str &
+	Util.string_delimit_list(tasks_str,", ") => t_str &
+	print " predecessor tasks :" & print t_str & print "\n" &
+
 	add_predecessors(t,tasks,strs,0) 
 	------------------------------
 	build_expression(Exp.CALL(func,expl,_,_)) => (t,"") 
@@ -603,9 +799,9 @@ relation add_predecessors: (int, int list,string list,int (* prio*)) => () =
 
 axiom add_predecessors(_,[],[],_) 
 
-  rule	TaskGraphExt.add_edge(t,t1,s,prio) &
+  rule	TaskGraphExt.add_edge(t1,t,s,prio) &
 	int_add(prio,1) => prio' &
 	add_predecessors(t,ts,strs,prio')
 	---------------------
 	add_predecessors(t,t1::ts,s::strs,prio)
-end
+end