many changes

-->
now works with Maple 7
correctly generates C code for AH LHS function


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@511 f88db872-0e4f-0410-b76b-b9085cfa78c5

4 files changed, 666 insertions, 12 deletions

diff --git a/src/maple/Diff.maple b/src/maple/Diff.maple
index 8755611..3740bb4 100644
--- a/src/maple/Diff.maple
+++ b/src/maple/Diff.maple
@@ -235,6 +235,9 @@ end proc;
 # It currently knows about the following simplifications:
 # - Diff(X_ud[u,i], x_xyz[j]) --> X_udd[u,i,j]
 #
+# Anything else is returned unchanged.  (To avoid infinite recursion,
+# such a return is *unevaluated*.)
+#
 # Arguments:
 # operand = (in) The thing to be differentiated.
 # var_seq = (in) (varargs) An expression sequence of the variables to
@@ -253,8 +256,8 @@ var_list := [args[2..nargs]];
 if ( type(operand, indexed) and (op(0,operand) = 'X_ud')
      and (nops(var_list) = 1) and member(var_list[1],x_xyz_list,'posn') )
    then return X_udd[op(operand), posn];
-   else return 'Diff'(operand, op(var_list));	# return is *unevaluated*
-						# to avoid infinite recursion!
 end if;
 
+# unevaluated return to avoid infinite recursion
+return 'Diff'(operand, op(var_list));
 end proc;
diff --git a/src/maple/codegen2.maple b/src/maple/codegen2.maple
new file mode 100644
index 0000000..a3a64a2
--- /dev/null
+++ b/src/maple/codegen2.maple
@@ -0,0 +1,589 @@
+# codegen2.maple -- generate C code from Maple expressions
+# $Id$
+#
+# codegen2 - generate C code from Maple expressions
+#
+# cvt_to_eqnlist - convert an expression to an equation list
+# fix_Diff - convert Diff() calls to PARTIAL_*() calls
+#   fix_Diff/remap_table - remapping table for fix_Diff()
+# temps_in_eqnlist - find temporaries used by an equation list
+# is_result - is a name equal to a "result" or an array/table component of it?
+# deindex_names - remove indices from indexed names
+# unindex_names - convert indexed names to scalars
+# fix_rationals - convert numbers to RATIONAL() calls
+# print_name_list_dcl - print a C declaration for a list of names
+#
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function is a high-level driver to generate C code for a Maple
+# expression.  It does the following:
+# - convert the input into a list of equations of the form  name = expression
+# - rewrite Diff(...) calls to DIFF_RHO_RHO(...) etc calls
+# - generate an "optimized" form of the computation sequence with
+#   codegen[optimize, tryhard]
+# - determine which temporary variables have been introduced by the
+#   codegen[optimize] library
+# - "unindex" all array accesses, eg R_dd[2,3] becomes R_dd_23.
+# - convert rational numbers to RATIONAL(p,q) calls
+# - print C declarations for the temporary variables
+# - print C code for the the optimized computation sequence
+#
+# Note that "codegen" is a Maple package; this function is called "codegen2".
+#
+# Sample usage:
+#	codegen2([R_dd__fnd, R__fnd], ['R_dd', 'R'], "Ricci.c");
+#
+# Arguments:
+# expr = (in) The expression or list of expressions for which code is to
+#	      be generated.  This will typically be the name of a gridfn
+#	      array functional-dependence form, or a list of such names,
+#	      but this isn't required.
+# lhs_name = (in) The name or list of names to be used for the results in
+#		  the generated C code, eg R_dd.
+# output_file_name = (in) The name of the file to which the generated
+#			  code is to be written.
+#
+# Arguments (as global variables)
+# `saveit/level`
+#	= (in) (optional)
+#	       If this global variable is assigned, intermediate results
+#	       are saved for debugging purposes.  If it's assigned an
+#	       integral value, making this value larger may increase the
+#	       level of debugging output.  (10 is a good number for typical
+#	       debugging.)
+#
+codegen2 :=
+proc(expr_in::{algebraic, list(algebraic)},
+     lhs_name::{name, list(name)},
+     output_file_name::string)
+local expr;
+
+printf("codegen(%a)\n", lhs_name);
+
+expr := expr_in;
+
+printf("   convert --> equation list\n");
+expr := cvt_to_eqnlist(expr, lhs_name);
+saveit(10, procname, "eqnlist", expr);
+
+printf("   optimizing computation sequence\n");
+expr := [codegen[optimize](expr)];
+##expr := [codegen[optimize](expr, tryhard)];
+saveit(10, procname, "optimize", expr);
+
+printf("   find temporary variables\n");
+temps := temps_in_eqnlist(expr, lhs_name);
+saveit(10, procname, "temps", expr);
+
+printf("   convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc\n");
+expr := fix_Diff(expr);
+saveit(10, procname, "fix_Diff", expr);
+
+printf("   convert R_dd[2,3] --> R_dd_23 etc\n");
+expr := unindex_names(expr);
+saveit(10, procname, "unindex", expr);
+
+printf("   convert p/q --> RATIONAL(p/q)\n");
+expr := fix_rationals(expr);
+saveit(10, procname, "fix_rationals", expr);
+
+printf("   generating C declarations for temporary variables\n");
+print_name_list_dcl(temps, "fp", output_file_name);
+printf("   generating C code\n");
+codegen[C](expr, filename=output_file_name);
+
+NULL;
+end proc;
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function converts an expression or list of expressions into an
+# equation list.  That is, given an expression  expr , this function
+# computes an equation list of the form
+#
+# if type(expr, algebraic)
+#	[ lhs_name = expr ]
+#
+# if type(expr, array)			# illustrated here for a rank-1 array
+#	[
+#	lhs_name[1] = expr[1],		# note equations are in lexicographic
+#	lhs_name[2] = expr[2],		# order of the indices, and include
+#	lhs_name[3] = expr[3],		# only those array elements that are
+#	    ...				# explicitly stored (as reported by
+#	lhs_name[N] = expr[N]		#  indices() )
+#	]
+#
+# if type(expr, list)
+#	then concatenate the equations lists from expr's elements
+#
+# Arguments:
+# expr = (in) The expression to be converted.
+# lhs_name = (in) The unevaluated name or list of names to use for the
+#		  left hand side(s) in the equation list.
+#
+# Results:
+# The equation list is returned as the function result.
+#
+cvt_to_eqnlist :=
+proc(expr::{algebraic, array, list({algebraic, array})},
+     lhs_name::{name, list(name)})
+
+# ... test for array first since otherwise expr itself is a "name",
+#     which would match type "algebraic" as well
+if   ( type(expr, array) and type(lhs_name, name) )
+   then return map(
+		      proc(ii)
+		      return lhs_name[op(ii)] = expr[op(ii)];
+		      end
+		    ,
+		      indices_in_order(expr)
+		  );
+
+elif ( type(expr, algebraic) and type(lhs_name, name) )
+   then return [lhs_name = expr];
+
+elif ( type(expr, list({algebraic, array})) and type(lhs_name, list(name)) )
+   then return zip(op @ cvt_to_eqnlist, expr, lhs_name);
+
+else error "unknown type for expression!\n"
+	   "   expr=%1\n"
+	   "   whattype(expr)=%2\n"
+	   ,
+	   expr, whattype(expr);
+fi;
+end;
+
+################################################################################
+
+#
+# This function converts  Diff()  calls into  PARTIAL_*()  calls, eg
+# Diff(src, rho, sigma) --> PARTIAL_RHO_SIGMA(src).
+#
+fix_Diff :=
+proc(expr::{algebraic, name = algebraic, list({algebraic, name = algebraic})})
+local nn, k, base, power, fn, fn_args_list, Darg, Dvars;
+global `fix_Diff/remap_table`;
+
+# recurse over lists
+if (type(expr, list))
+   then return map(fix_Diff, expr);
+fi;
+
+# recurse over equation right hand sides
+if (type(expr, name = algebraic))
+   then return lhs(expr) = fix_Diff(rhs(expr));
+fi;
+
+nn := nops(expr);
+
+# recurse over sums
+if (type(expr, `+`))
+   then return sum('fix_Diff(op(k,expr))', 'k'=1..nn);
+fi;
+
+# recurse over products
+if (type(expr, `*`))
+   then return product('fix_Diff(op(k,expr))', 'k'=1..nn);
+fi;
+
+# recurse over powers
+if (type(expr, `^`))
+   then 
+	base := op(1, expr);
+	power := op(2, expr);
+
+	return fix_Diff(base) ^ power;
+	fi;
+
+# recurse over non-Diff functions
+if type(expr, function) and (op(0, expr) <> 'Diff')
+   then 
+	fn := op(0, expr);
+	fn_args_list := [op(expr)];
+	
+	fn; return '%'( op(map(fix_Diff, fn_args_list)) );
+fi;
+
+# remap derivatives
+if type(expr, function) and (op(0, expr) = 'Diff')
+   then
+	Darg := op(1, expr);
+	Dvars := [op(2..nn, expr)];
+	if (assigned(`fix_Diff/remap_table`[op(Dvars)]))
+	   then `fix_Diff/remap_table`[op(Dvars)]; return '%'(Darg);
+	   else error "don't know how to remap Diff() call!\n"
+		      "   Darg = %1\n"
+		      "   Dvars = %2\n"
+		      ,
+		      Darg, Dvars;
+	fi;
+fi;
+
+# otherwise, the identity function
+return expr;
+end;
+
+########################################
+
+#
+# this table defines the remapping of Diff() calls for  fix_Diff()  (above)
+# n.b. Diff() should already have canonicalized the order of variables
+#
+`fix_Diff/remap_table`[rho  ] := 'PARTIAL_RHO';
+`fix_Diff/remap_table`[sigma] := 'PARTIAL_SIGMA';
+`fix_Diff/remap_table`[rho  , rho  ] := 'PARTIAL_RHO_RHO';
+`fix_Diff/remap_table`[rho  , sigma] := 'PARTIAL_RHO_SIGMA';
+`fix_Diff/remap_table`[sigma, sigma] := 'PARTIAL_SIGMA_SIGMA';
+
+`fix_Diff/remap_table`[xx] := 'PARTIAL_X';
+`fix_Diff/remap_table`[yy] := 'PARTIAL_Y';
+`fix_Diff/remap_table`[zz] := 'PARTIAL_Z';
+`fix_Diff/remap_table`[xx,xx] := 'PARTIAL_XX';
+`fix_Diff/remap_table`[xx,yy] := 'PARTIAL_XY';
+`fix_Diff/remap_table`[xx,zz] := 'PARTIAL_XZ';
+`fix_Diff/remap_table`[yy,yy] := 'PARTIAL_YY';
+`fix_Diff/remap_table`[yy,zz] := 'PARTIAL_YZ';
+`fix_Diff/remap_table`[zz,zz] := 'PARTIAL_ZZ';
+
+################################################################################
+
+#
+# Given an equation list, this function finds all the temporaries
+# assigned by it.  A "temporary" is defined here to be a name on the
+# left hand side of an equation, which isn't the result or a component
+# of it.
+#
+# Arguments:
+# expr = (in) The equation list to operate on.
+# result_name = (in) The result name or list/set of result names.
+#
+# Results:
+# The function returns the list of temporaries assigned.
+#
+temps_in_eqnlist :=
+proc(expr::list(name = algebraic),
+     result_name::{name, list(name), set(name)})
+
+# "temporary" = lhs name which isn't a result
+return remove(is_result, map(lhs,expr), result_name);
+end;
+
+################################################################################
+
+#
+# This function tests whether or not a name is a "result" name or
+# an array/table component of it.  Either a single result name, or a
+# list/set of these, may be specified; in the latter case the function
+# tests whether or not a name matches *any* of the result names.
+#
+# Arguments:
+# rname = (in) The name to test.
+# result_name = (in) The name or list/set of names of the result to test
+#		     against.
+#
+# Results:
+# The function returns  true  if the name is equal to the result or an
+# array/table component of it,  false  otherwise.
+#
+is_result :=
+proc(rname::name,
+     result_name_in::{name, list(name), set(name)})
+local result_name, rn;
+
+if type(in_result_name, name)
+   then result_name := { result_name_in };
+   else result_name := result_name_in;
+fi;
+
+	for rn in result_name
+	do
+	if (rname = rn)
+	   then return true;
+	elif (type(rname, indexed) and (op(0,rname) = rn))
+	   then return true;
+	fi;
+	od;
+
+return false;
+end;
+
+################################################################################
+
+#
+# This function removes all indices from indexed names, eg
+#	A[1,2,3] --> A .
+#
+# Arguments:
+# expr = (in) The expression to be converted.
+#
+# Results:
+# The converted expression is returned as the function result.
+#
+deindex_names :=
+proc(expr::{name, list(name), set(name)})
+
+# recurse over lists and sets
+if (type(expr, {list, set}))
+   then return map(deindex_names, expr);
+fi;
+
+# return non-indexed names and numbers unchanged
+if (type(expr, {string, numeric}))
+   then return expr;
+fi;
+
+# convert indexed names
+if (type(expr, indexed))
+   then return op(0, expr);
+fi;
+
+# unknown type
+error "expr has unknown type!\n"
+      "whattype(expr)=%1\n"
+      ,
+      whattype(expr);
+end;
+
+################################################################################
+
+#
+# This function converts all occurence of indexed names in an expression
+# to new non-indexed "scalar names" of the form
+#	A[1,2,3] --> A_123 .
+#
+# Arguments:
+# expr = (in) The expression to be converted.
+#
+# Results:
+# The converted expression is returned as the function result.
+#
+unindex_names :=
+proc(expr::{
+	   algebraic, name = algebraic,
+	   list({algebraic, name = algebraic}),
+	   set({algebraic, name = algebraic})
+	   })
+local nn, k,
+      base, power,
+      fn_name, fn_args_list, converted_fn_args_list,
+      base_name, index_seq;
+
+# recurse over lists and sets
+if (type(expr, {list, set}))
+   then return map(unindex_names, expr);
+fi;
+
+# recurse over equations (both lhs and rhs)
+if (type(expr, `=`))
+   then return unindex_names(lhs(expr))  =  unindex_names(rhs(expr));
+fi;
+
+nn := nops(expr);
+
+# recurse over sums
+if (type(expr, `+`))
+   then return sum('unindex_names(op(k,expr))', 'k'=1..nn);
+fi;
+
+# recurse over products
+if (type(expr, `*`))
+   then return product('unindex_names(op(k,expr))', 'k'=1..nn);
+fi;
+
+# recurse over powers
+if (type(expr, `^`))
+   then
+	base := op(1, expr);
+	power := op(2, expr);
+	return unindex_names(base) ^ power;
+fi;
+
+# recurse over function calls
+if (type(expr, function))
+   then
+	fn_name := op(0, expr);
+	fn_args_list := [op(expr)];
+	converted_fn_args_list := map(unindex_names, fn_args_list);
+	return fn_name(op(converted_fn_args_list));
+fi;
+
+# convert indexed names
+if (type(expr, indexed))
+   then
+	base_name := op(0, expr);
+	index_seq := op(expr);
+	return cat(base_name,"_",index_seq);
+fi;
+
+# return numbers and non-indexed names
+if (type(expr, {numeric, name}))
+   then return expr;
+fi;
+
+# unknown type
+error "expr has unknown type!\n"
+      "whattype(expr)=%1\n"
+      ,
+      whattype(expr);
+end;
+
+################################################################################
+
+#
+# This function converts all {integer, rational} subexpressions of its
+# input except integer exponents and -1 factors in products, into function
+# calls  RATIONAL(num,den)  with  num  and  den  integers.
+#
+# This is useful in conjunction with the  C() library function, since
+#
+#	C( (1/3) * foo * bar )
+#		t0 = foo*bar/3;
+#
+# generates a (slow) division (and runs the risk of mixed-mode-arithmetic
+# problems).  In contrast, with this function
+#
+#	fix_rationals((1/3) * foo * bar);
+#	     RATIONAL(1,3) foo bar
+#	codegen[C](%);
+#	     t0 = RATIONAL(1.0,3.0)*foo*bar;
+#
+# which a C preprocessor macro can easily convert to the desired
+#
+#	     t0 = (1.0/3.0)*foo*bar;
+#
+# Arguments:
+# expr = (in) The expression to be converted.
+#
+fix_rationals :=
+proc(expr::{algebraic, name = algebraic, list({algebraic, name = algebraic})})
+local nn, k,
+      base, power, fbase, fpower,
+      fn, fn_args_list, fixed_fn_args_list,
+      num, den, mult;
+
+# recurse over lists
+if (type(expr, list))
+   then return map(fix_rationals, expr);
+fi;
+
+# recurse over equation right hand sides
+if (type(expr, name = algebraic))
+   then return lhs(expr) = fix_rationals(rhs(expr));
+fi;
+
+# recurse over functions other than  RATIONAL()
+if (type(expr, function))
+   then
+	fn := op(0, expr);
+	if (fn <> 'RATIONAL')
+	   then
+		fn_args_list := [op(expr)];
+		fixed_fn_args_list := map(fix_rationals, fn_args_list);
+		fn; return '%'(op(fixed_fn_args_list));
+	fi;
+fi;
+
+nn := nops(expr);
+
+# recurse over sums
+if (type(expr, `+`))
+   then return sum('fix_rationals(op(k,expr))', 'k'=1..nn);
+fi;
+
+# recurse over products
+# ... leaving -1 factors intact
+if (type(expr, `*`))
+   then return product('fix_rationals(op(k,expr))', 'k'=1..nn);
+fi;
+
+# recurse over powers
+# ... leaving integer exponents intact
+if (type(expr, `^`))
+   then
+	base := op(1, expr);
+	power := op(2, expr);
+
+	fbase := fix_rationals(base);
+	if (type(power, integer))
+	   then fpower := power;
+	   else fpower := fix_rationals(power);
+	fi;
+	return fbase ^ fpower;;
+fi;
+
+# fix integers and fractions
+# ... leaving explicit -1 intact
+if (type(expr, integer))
+   then return 'RATIONAL'(expr, 1);
+fi;
+if (type(expr, fraction))
+   then
+	num := op(1, expr);
+	den := op(2, expr);
+	return 'RATIONAL'(num, den);
+fi;
+
+# turn Maple floating-point into integer fraction, then recursively fix that
+if (type(expr, float))
+   then
+	mult := op(1, expr);
+	power := op(2, expr);
+	return fix_rationals(mult * 10^power);
+fi;
+
+
+
+# identity op on names
+if (type(expr, name))
+   then return expr;
+fi;
+
+# unknown type
+error "expr has unknown type!\n"
+      "whattype(expr)=%1\n"
+      "expr=%2\n"
+      ,
+      whattype(expr), expr;
+end;
+
+################################################################################
+
+#
+# This function prints a C declaration for a list of names.
+#
+# Argument:
+# name_list = A list of the names.
+# name_type = The C type of the names, eg. "double".
+# file_name = The file name to write the declaration to.  This is
+#	      truncated before writing.
+#
+print_name_list_dcl :=
+proc( name_list::list({name,string}),
+      name_type::string,
+      file_name::string )
+local blanks, separator_string;
+
+ftruncate(file_name);
+
+# a sequence of blanks with the same length as name_type
+seq(" ", i=1..length(name_type));
+
+# string to separate names
+separator_string := cat(",\n", %, " ");
+
+map(convert, name_list, string);
+ListTools[Join](%, separator_string);
+cat(op(%));
+
+fprintf(file_name,
+	"%s %s;\n",
+	name_type, %);
+NULL;
+end proc;
diff --git a/src/maple/setup.maple b/src/maple/setup.maple
index c23e547..6835cb6 100644
--- a/src/maple/setup.maple
+++ b/src/maple/setup.maple
@@ -24,3 +24,5 @@ read "../maple/gfa.mm";
 
 read "../maple/coeffs.mm";
 setup_coeff_gfas();
+
+read "../maple/codegen2.mm";
diff --git a/src/maple/util.maple b/src/maple/util.maple
index fe7e0ea..2cf6562 100644
--- a/src/maple/util.maple
+++ b/src/maple/util.maple
@@ -6,6 +6,8 @@
 # arctan_xy - 2-argument arctan() function
 # ssqrt - symbolic sqrt()
 #
+# indices_in_order - list of indices of table/array, in lexicographic order
+# lexorder_integer_list - lexorder() for lists
 # sort_var_list - sort a list of variables into a canonical order
 # lexorder_vars - lexorder() for "interesting" variables
 #
@@ -95,6 +97,63 @@ end proc;
 ################################################################################
 
 #
+# This function is identical to Maple's  indices()  built-in function,
+# except that (a) it returns a list, not an expression sequence, and
+# (b) the indices are returned in lexicographic order rather than in
+# Maple's internal hash ordering.
+#
+# Arguments:
+# T = (in) A table or array with (only) integer-sequence indices.
+#
+# Results:
+# The lexicographic-ordered expression sequence of indices is returned
+# as the function result.
+#
+# Bugs:
+# - The function fails if the indices aren't a list of numbers.
+#
+lindices_in_order :=
+proc(T::table)
+return sort([indices(T)], lexorder_integer_list);
+end;
+
+###############################################################################
+
+#
+# This function lexicographically orders lists of integers.
+#
+# Arguments:
+# list1 = list of integers to be sorted
+# list2 = list of integers to be sorted
+#
+# Results:
+# The function returns  true  iff  list1 < list2  lexicographically,
+# with the [1] list elements being most significant.
+#
+lexorder_integer_list :=
+proc(list1::list(numeric), list2::list(numeric))
+local len1, len2, k;
+
+len1 := nops(list1);
+len2 := nops(list2);
+
+	for k from 1 to min(len1,len2)
+	do
+	if   (list1[k] < list2[k])
+	   then return true;
+	elif (list1[k] > list2[k])
+	   then return false;
+	fi;
+	od;
+
+# get to here ==> the shorter list is an exact prefix of the longer one
+#             ==> order the shorter one < the longer one
+return evalb(len1 < len2);
+end;
+
+################################################################################
+
+#
 # This function sorts a list of names (in practice, a list of coordinates)
 # into a canonical order.
 #
@@ -105,7 +164,7 @@ global lexorder_vars;
 
 # only get to here the first time we see a given variable list
 # (i.e. if it's not in the remember table)
-RETURN( sort(var_list, lexorder_vars) );
+return sort(var_list, lexorder_vars);
 end;
 
 ################################################################################
@@ -217,12 +276,13 @@ end proc;
 # Maple doesn't have a proper debugger :-( .
 #
 # Arguments (explicit):
-# fn_name = (in) The calling function or subsystem's name.
+# n = (in) An integer controlling how important this saveit() call is;
+#	   the saving will only be done if `saveit/level` >= n, i.e.
+#	   smaller values of n make the saving more likely.
+# fn = (in) The calling function or subsystem's name.
+# label = (in) A string identifying this saveit() call among all those
+#	       within the calling function.
 # expr = (in) The temporary result to save.
-# name_postfix... = (in) (varargs)
-#			 One or more arguments which this function
-#			 concatenates to form the name postfix for
-#			 the global variable.
 #
 # Arguments (implicit, as global variables)
 # `saveit/level` = (in) (optional) If this name is assigned (as determined
@@ -231,7 +291,7 @@ end proc;
 #				  this function is a nop.
 #
 saveit :=
-proc(n::integer, fn_name::string, label::string, expr::anything)
+proc(n::integer, fn::{procedure,string}, label::string, expr::anything)
 global `saveit/level`;
 local save_name;
 
@@ -239,9 +299,9 @@ if ( assigned(`saveit/level`)
      and type(`saveit/level`, integer)
      and (`saveit/level` >= n) )
    then
-	save_name := cat(fn_name,`/`,label);
-	lprint(printf(`      --> \"%s\"`, save_name));
-	assign(  eval(save_name,1) = expr  );
+	save_name := cat(convert(fn,string),"/",label);
+	printf("      --> `%s`\n", save_name);
+	assign(  convert(eval(save_name,1),name) = expr  );
 end if;
 
 NULL;