reorganize common to Lagrange and Hermite interpolators

(i.e. code depending only on molecule size/shape)
--> now there are separate *functions* rather than just code fragments
--> ../template.c should hopefully compile in finite time.....


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@109 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 294 insertions, 365 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/common/1d.log b/src/GeneralizedPolynomial-Uniform/common/1d.log
index 4c1a387..845a8c1 100644
--- a/src/GeneralizedPolynomial-Uniform/common/1d.log
+++ b/src/GeneralizedPolynomial-Uniform/common/1d.log
@@ -10,7 +10,7 @@
 # fix_rationals - convert numbers to RATIONAL() calls
 # nonmatching_names - find names in a list which *don't* have a specified prefix
 # sprint_numeric_list - convert a numeric list to a valid C identifier suffix
-# print_name_list_dcl - print a C declaration for a list of names
+# print_name_list_dcl - print C declarations for a list of names
 #
 # hypercube_points - compute all (integer) points in an N-dimensional hypercube
 #
@@ -297,49 +297,42 @@ end proc
 ################################################################################
 > 
 #
-# This function prints a C declaration for a list of names.
+# This function prints a sequence of C declarations for a list of names.
 #
 # Argument:
 # name_list = A list of the names.
-# name_type = The C type of the names, eg. "double".
+# type_name = 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,
+>       type_name::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, %);
+> map(
+>        proc(var::{name,string})
+>        fprintf(file_name,
+> 	       "%s %s;\n", 
+> 	       type_name, var);
+>        end proc
+>      ,
+>        name_list
+>    );
 > 
 > fclose(file_name);
 > NULL;
 > end proc;
 print_name_list_dcl := proc(
-name_list::list({name, string}), name_type::string, file_name::string)
+name_list::list({name, string}), type_name::string, file_name::string)
 local blanks, separator_string;
     ftruncate(file_name);
-    seq(" ", i = 1 .. length(name_type));
-    separator_string := cat(",\n", %, " ");
-    map(convert, name_list, string);
-    ListTools[Join](%, separator_string);
-    cat(op(%));
-    fprintf(file_name, "%s %s;\n", name_type, %);
+    map(proc(var::{name, string})
+            fprintf(file_name, "%s %s;\n", type_name, var)
+        end proc, name_list);
     fclose(file_name);
     NULL
 end proc
@@ -445,7 +438,7 @@ ftruncate :=
     proc(file_name::string) fopen(file_name, 'WRITE'); fclose(%); NULL end proc
 
 # interpolate.maple -- compute interpolation formulas/coefficients
-# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/interpolate.maple,v 1.8 2002/08/20 16:46:05 jthorn Exp $
+# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/interpolate.maple,v 1.10 2002/08/28 11:31:09 jthorn Exp $
 > 
 #
 # <<<representation of numbers, data values, etc>>>
@@ -454,8 +447,8 @@ ftruncate :=
 # coeff_as_lc_of_data - coefficients of ... (linear combination of data)
 #
 # print_coeff__lc_of_data - print C code to compute coefficients
-# print_data_var_assign - print C code to assign data-value variables
-# print_interp_coeff_var_store - print C code to store coeff vars "somewhere"
+# print_fetch_data - print C code to fetch input array chunk into struct data
+# print_store_coeffs - print C code to store struct coeffs "somewhere"
 # print_interp_cmpt__lc_of_data - print C code for computation of interpolant
 #
 # coeff_name - name of coefficient of data at a given [m] coordinate
@@ -543,7 +536,7 @@ end proc
 #
 # This function computes a Hermite polynomial interpolant in any
 # number of dimensions.  This is a polynomial which
-#  has values which match the given data DATA() at a specified set of
+# * has values which match the given data DATA() at a specified set of
 #   points, and
 # * has derivatives which match the specified finite-difference derivatives
 #   of the given data DATA() at a specified set of points
@@ -562,50 +555,45 @@ end proc
 #		+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y
 #		+ c00     + c10*x     + c20*x^2     + c30*x^3
 #		end proc;
-# coeff_list = A set of the interpolation coefficients (coefficients in
+# coeff_set = A set of the interpolation coefficients (coefficients in
 #	       the interpolation function), for example
-#			[
+#			{
 #			c03, c13, c23, c33,
 #			c02, c12, c22, c32,
 #			c01, c11, c21, c31,
 #			c00, c10, c20, c30
-#			]
+#			}
 # coord_list = A list of the coordinates (independent variables in the
 #	       interpolation function), for example [x,y].
-# deriv_coord_list = A list of lists of coordinates specifying which
-#		     derivatives are computed by the deriv_proc_list[]
-#		     procedures, for example
-#			[[x], [y], [x,y]]
-# deriv_proc_list = A list of procedures for computing finite-difference
-#		    derivatives.  Each procedure should take N_dims integer
-#		    arguments specifying an evaluation point, and return
-#		    a suitable linear combination of the DATA() for the
-#		    derivative at that point.  For example
-#		    example,
-#			[
-#			  proc(i::integer, j::integer)
-#			  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
-#			  end proc
+# deriv_set = A set of equations of the form
+#		{coords} = proc
+#	      giving the derivatives which are to be matched, and the
+#	      procedures to compute their finite-difference approximations.
+#	      Each procedure should take N_dims integer arguments specifying
+#	      an evaluation point, and return a suitable linear combination
+#	      of the DATA() for the derivative at that point.  For example
+#			{
+#			  {x}   = proc(i::integer, j::integer)
+#				  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
+#				  end proc
 #			,
-#			  proc(i::integer, j::integer)
-#			  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
-#			  end proc
+#			  {y}   = proc(i::integer, j::integer)
+#				  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
+#				  end proc
 #			,
-#			  proc(i::integer, j::integer)
-#			  - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
-#			  + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
-#			  end proc
-#			]
-# fn_posn_list = A list of positions (each a list of numeric values)
-#		 where the interpolant is to match the given data DATA(),
-#		 for example
-#			[[0,0], [0,1], [1,0], [1,1]]
-# deriv_posn_list = A list of positions (each a list of numeric values)
-#		    where the interpolant is to match each possible product
-#		    of 1st derivatives"1st"
-#		    derivatives (actually all derivatives the given data DATA(),
-#		    for example
-#			[[0,0], [0,1], [1,0], [1,1]]
+#			  {x,y} = proc(i::integer, j::integer)
+#				  - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
+#				  + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
+#				  end proc
+#			}
+# fn_posn_set = A set of positions (each a list of numeric values)
+#		where the interpolant is to match the given data DATA(),
+#		for example
+#			{[0,0], [0,1], [1,0], [1,1]}
+# deriv_posn_set = A list of positions (each a list of numeric values)
+#		   where the interpolant is to match the derivatives
+#		   specified by  deriv_set , for example
+#			{[0,0], [0,1], [1,0], [1,1]}
 #
 # Results:
 # This function returns the interpolating polynomial, in the form of
@@ -613,93 +601,123 @@ end proc
 #
 > Hermite_polynomial_interpolant :=
 > proc(
->       fn::procedure, coeff_list::list(name), coord_list::list(name),
->       deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure),
->       fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric))
+>       fn::procedure,
+>       coeff_set::set(name),
+>       coord_list::list(name),
+>       deriv_set::set(set(name) = procedure),
+>       fn_posn_set::set(list(numeric)),
+>       deriv_posn_set::set(list(numeric))
 >     )
-> local fn_eset,
->       fn_expr, deriv_eset;
-> 
-# set of equations {fn(posn) = DATA(posn)}
-> fn_eset := map(
-> 		  # return equation that fn(this point) = DATA(this point)
-> 		  proc(posn::list(integer))
-> 		  fn(op(posn)) = 'DATA'(op(posn));
-> 		  end proc
-> 		,
-> 		  {op(fn_posn_list)}
-> 	      );
-> 
-# set of sets of equations
-#	{
-#	  { deriv1(posn1) = fd_deriv1(posn1),
-#	    deriv2(posn1) = fd_deriv2(posn1),
-#	    ... },
-#	  { deriv1(posn2) = fd_deriv1(posn2),
-#	    deriv2(posn2) = fd_deriv2(posn2),
-#	    ... },
-#	  ...
-#	}
-> fn_expr := fn(op(coord_list));
+> local fn_eqnset, deriv_eqnset, coeff_eqns, subs_eqnset;
+> 
+> 
+#
+# compute a set of equations
+#	{fn(posn) = DATA(posn)}
+# giving the function values to be matched
+#
+> fn_eqnset := map(
+> 		    # return equation that fn(posn) = DATA(posn)
+> 		    proc(posn::list(integer))
+> 		    fn(op(posn)) = 'DATA'(op(posn));
+> 		    end proc
+> 		  ,
+> 		    fn_posn_set
+> 		);
+> 
+> 
+#
+# compute a set of equations
+#	{ diff(fn,coords)(posn) = DERIV(coords)(posn) }
+# giving the derivative values to be matched, where DERIV(coords)
+# is a placeholder for the appropriate derivative
+#
 > map(
->        # return set of equations
->        #	{deriv(posn) = fd_deriv(posn)}
->        # for this point
->        proc(posn::list(integer))
->        {
-> 	op(
-> 	 zip(
-> 		# return equation that
-> 		#	deriv(posn) = fd_deriv(posn)
-> 		# for this deriv and this point
-> 		proc(deriv_coords::list(name), fd_deriv_proc::procedure)
-> 		local fn_deriv_proc;
-> 		fn_deriv_proc := unapply( diff(fn_expr,deriv_coords),
-> 					  op(coord_list) );
-> 		fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn));
-> 		end proc
-> 	      ,
-> 		[op(deriv_coord_list)]
-> 	      ,
-> 		[op(fd_deriv_proc_list)]
-> 	    )
-> 	  )
->        }
+>        # return set of equations for this particular derivative
+>        proc(deriv_coords::set(name))
+>        local deriv_fn;
+>        fn(op(coord_list));
+>        diff(%, op(deriv_coords));
+>        deriv_fn := unapply(%, op(coord_list));
+>        map(
+> 	      proc(posn::list(integer))
+> 	      deriv_fn(op(posn)) = 'DERIV'(op(deriv_coords))(op(posn));
+> 	      end proc
+> 	    ,
+> 	      deriv_posn_set
+> 	  );
 >        end proc
 >      ,
->        {op(deriv_posn_list)}
+>        map(lhs, deriv_set)
 >    );
+> deriv_eqnset := `union`(op(%));
+> 
+> 
+#
+# solve overall set of equations for coefficients
+# in terms of DATA() and DERIV() values
+#
+> coeff_eqns := solve[linear](fn_eqnset union deriv_eqnset, coeff_set);
+> if (indets(map(rhs,%)) <> {})
+>    then error "no unique solution for coefficients -- %1 eqns for %2 coeffs",
+> 	      nops(fn_eqnset union deriv_eqnset),
+> 	      nops(coeff_set);
+> fi;
 > 
-# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)}
-> deriv_eset := `union`(op(%));
 > 
-# solve equations for coeffs
-> solve(fn_eset union deriv_eset, {op(coeff_list)});
+#
+# compute a set of substitution equations
+#	{'DERIV'(coords) = procedure}
+#
+> subs_eqnset := map(
+> 		      proc(eqn::set(name) = procedure)
+> 		      'DERIV'(op(lhs(eqn))) = rhs(eqn);
+> 		      end proc
+> 		    ,
+> 		      deriv_set
+> 		  );
 > 
-# interpolant as a polynomial in the coordinates
-> return subs(%, eval(fn))(op(coord_list));
+> 
+#
+# compute the coefficients in terms of the DATA() values
+#
+> subs(subs_eqnset, coeff_eqns);
+> eval(%);
+> 
+#
+# compute the interpolant as a polynomial in the coordinates
+#
+> subs(%, fn(op(coord_list)));
 > end proc;
-Hermite_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
-coord_list::list(name), deriv_coord_list::list(list(name)),
-fd_deriv_proc_list::list(procedure), fn_posn_list::list(list(numeric)),
-deriv_posn_list::list(list(numeric)))
-local fn_eset, fn_expr, deriv_eset;
-    fn_eset := map(
+Hermite_polynomial_interpolant := proc(fn::procedure, coeff_set::set(name),
+coord_list::list(name), deriv_set::set(set(name) = procedure),
+fn_posn_set::set(list(numeric)), deriv_posn_set::set(list(numeric)))
+local fn_eqnset, deriv_eqnset, coeff_eqns, subs_eqnset;
+    fn_eqnset := map(
         proc(posn::list(integer)) fn(op(posn)) = 'DATA'(op(posn)) end proc,
-        {op(fn_posn_list)});
-    fn_expr := fn(op(coord_list));
-    map(proc(posn::list(integer))
-            {op(zip(proc(deriv_coords::list(name), fd_deriv_proc::procedure
-            )
-            local fn_deriv_proc;
-                fn_deriv_proc :=
-                    unapply(diff(fn_expr, deriv_coords), op(coord_list));
-                fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn))
-            end proc, [op(deriv_coord_list)], [op(fd_deriv_proc_list)]))}
-        end proc, {op(deriv_posn_list)});
-    deriv_eset := `union`(op(%));
-    solve(fn_eset union deriv_eset, {op(coeff_list)});
-    return subs(%, eval(fn))(op(coord_list))
+        fn_posn_set);
+    map(proc(deriv_coords::set(name))
+        local deriv_fn;
+            fn(op(coord_list));
+            diff(%, op(deriv_coords));
+            deriv_fn := unapply(%, op(coord_list));
+            map(proc(posn::list(integer))
+                    deriv_fn(op(posn)) =
+                    'DERIV'(op(deriv_coords))(op(posn))
+                end proc, deriv_posn_set)
+        end proc, map(lhs, deriv_set));
+    deriv_eqnset := `union`(op(%));
+    coeff_eqns := solve[linear](fn_eqnset union deriv_eqnset, coeff_set);
+    if indets(map(rhs, %)) <> {} then error
+        "no unique solution for coefficients -- %1 eqns for %2 coeffs",
+        nops(fn_eqnset union deriv_eqnset), nops(coeff_set)
+    end if;
+    subs_eqnset := map(proc(eqn::(set(name) = procedure))
+            'DERIV'(op(lhs(eqn))) = rhs(eqn)
+        end proc, deriv_set);
+    subs(subs_eqnset, coeff_eqns);
+    eval(%);
+    subs(%, fn(op(coord_list)))
 end proc
 
 > 
@@ -858,7 +876,7 @@ end proc
 #
 # This function prints a sequence of C expression to assign the data-value
 # variables, eg
-#	data_m1_p1 = DATA(-1,1);
+#	data->data_m1_p1 = DATA(-1,1);
 #
 # Arguments:
 # posn_list = The same list of positions as was used to compute the
@@ -867,7 +885,7 @@ end proc
 # file_name = The file name to write the coefficients to.  This is
 #	      truncated before writing.
 #
-> print_data_var_assign :=
+> print_fetch_data :=
 > proc(
 >       posn_list::list(list(numeric)),
 >       data_var_name_prefix::string,
@@ -889,7 +907,7 @@ end proc
 > 
 > NULL;
 > end proc;
-print_data_var_assign := proc(posn_list::list(list(numeric)),
+print_fetch_data := proc(posn_list::list(list(numeric)),
 data_var_name_prefix::string, file_name::string)
     ftruncate(file_name);
     map(proc(posn::list(numeric))
@@ -906,25 +924,18 @@ end proc
 #
 # This function prints a sequence of C expression to store the interpolation
 # coefficients in  COEFF(...)  expressions, eg
-#	COEFF(1,-1) = factor * coeff_dx_p1_m1;
+#	COEFF(1,-1) = factor * coeffs->coeff_p1_m1;
 #
 # Arguments:
-# posn_list = The same list of positions as was used to compute the
-#	      interpolating polynomial.
-# RHS_factor_name = If this string is non-empty, then the coefficient is
-#		    multiplied by this factor before being stored, eg
-#		    setting this to "factor" would give the example above.
-#		    If this string is empty (""), the multiplication is
-#		    omitted, eg
-#			COEFF(1,-1) = coeff_dx_p1_m1;
-# coeff_name_prefix = A prefix string for the coefficient names.
+# posn_list = The list of positions in the molecule.
+# coeff_name_prefix = A prefix string for the coefficient names,
+#		      eg "factor * coeffs->coeff_"
 # file_name = The file name to write the coefficients to.  This is
 #	      truncated before writing.
 #
-> print_interp_coeff_var_store :=
+> print_store_coeffs :=
 > proc(
 >       posn_list::list(list(numeric)),
->       RHS_factor_name::string,
 >       coeff_name_prefix::string,
 >       file_name::string
 >     )
@@ -932,17 +943,10 @@ end proc
 > ftruncate(file_name);
 > map(
 >        proc(posn::list(numeric))
->        if (length(RHS_factor_name) > 0)
-> 	  then fprintf(file_name,
-> 		       "%a = %s * %s;\n",
-> 		       'COEFF'(op(posn)),
-> 		       RHS_factor_name,
-> 		       coeff_name(posn,coeff_name_prefix));
-> 	  else fprintf(file_name,
-> 		       "%a = %s;\n",
-> 		       'COEFF'(op(posn)),
-> 		       coeff_name(posn,coeff_name_prefix));
->        end if;
+>        fprintf(file_name,
+> 	       "%a = %s;\n",
+> 	       'COEFF'(op(posn)),
+> 	       coeff_name(posn,coeff_name_prefix));
 >        end proc
 >      ,
 >        posn_list
@@ -951,16 +955,12 @@ end proc
 > 
 > NULL;
 > end proc;
-print_interp_coeff_var_store := proc(posn_list::list(list(numeric)),
-RHS_factor_name::string, coeff_name_prefix::string, file_name::string)
+print_store_coeffs := proc(posn_list::list(list(numeric)),
+coeff_name_prefix::string, file_name::string)
     ftruncate(file_name);
     map(proc(posn::list(numeric))
-            if 0 < length(RHS_factor_name) then fprintf(file_name,
-                "%a = %s * %s;\n", 'COEFF'(op(posn)), RHS_factor_name,
-                coeff_name(posn, coeff_name_prefix))
-            else fprintf(file_name, "%a = %s;\n", 'COEFF'(op(posn)),
-                coeff_name(posn, coeff_name_prefix))
-            end if
+            fprintf(file_name, "%a = %s;\n", 'COEFF'(op(posn)),
+            coeff_name(posn, coeff_name_prefix))
         end proc, posn_list);
     fclose(file_name);
     NULL
@@ -970,25 +970,19 @@ end proc
 ################################################################################
 > 
 #
-# This function prints a C expression to compute the interpolant,
-# using the coefficients computed by  print_coeff__lc_of_data()
-# (i.e. expressing the interpolant as a linear combination of the
-# data values).
+# This function prints a C expression to evaluate a molecule, i.e.
+# to compute the molecule as a linear combination of the data values.
 #
 # Arguments:
-# posn_list = The same list of positions as was used to compute the
-#	      interpolating polynomial.
-# result_var_name = The (string) name of the variable to which the
-#		    result is to be assigned.
+# posn_list = The list of positions in the molecule.
 # coeff_name_prefix = A prefix string for the coefficient names.
 # data_var_name_prefix = A prefix string for the data variable names.
 # file_name = The file name to write the coefficients to.  This is
 #	      truncated before writing.
 #
-> print_interp_cmpt__lc_of_data :=
+> print_evaluate_molecule :=
 > proc(
 >       posn_list::list(list(numeric)),
->       result_var_name::string,
 >       coeff_name_prefix::string,
 >       data_var_name_prefix::string,
 >       file_name::string
@@ -996,8 +990,6 @@ end proc
 > 
 > ftruncate(file_name);
 > 
-> fprintf(file_name, "%s =\n", result_var_name);
-> 
 # list of "coeff*data_var" terms
 > map(
 >        proc(posn::list(numeric))
@@ -1009,26 +1001,24 @@ end proc
 >        posn_list
 >    );
 > 
-> ListTools[Join](%, "\n\t+ ");
+> ListTools[Join](%, "\n  + ");
 > cat(op(%));
-> fprintf(file_name, "\t%s;\n", %);
+> fprintf(file_name, "    %s;\n", %);
 > 
 > fclose(file_name);
 > 
 > NULL;
 > end proc;
-print_interp_cmpt__lc_of_data := proc(posn_list::list(list(numeric)),
-result_var_name::string, coeff_name_prefix::string,
-data_var_name_prefix::string, file_name::string)
+print_evaluate_molecule := proc(posn_list::list(list(numeric)),
+coeff_name_prefix::string, data_var_name_prefix::string, file_name::string)
     ftruncate(file_name);
-    fprintf(file_name, "%s =\n", result_var_name);
     map(proc(posn::list(numeric))
             sprintf("%s*%s", coeff_name(posn, coeff_name_prefix),
             data_var_name(posn, data_var_name_prefix))
         end proc, posn_list);
-    ListTools[Join](%, "\n\t+ ");
+    ListTools[Join](%, "\n  + ");
     cat(op(%));
-    fprintf(file_name, "\t%s;\n", %);
+    fprintf(file_name, "    %s;\n", %);
     fclose(file_name);
     NULL
 end proc
@@ -1331,31 +1321,26 @@ posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 # generic stuff for 1d, cube, size=2
 #
 > 
-> data_var_list_1d_size2 := map(data_var_name, posn_list_1d_size2, "data_");
-                data_var_list_1d_size2 := ["data_0", "data_p1"]
+> data_list_1d_size2 := map(data_var_name, posn_list_1d_size2, "data_");
+                  data_list_1d_size2 := ["data_0", "data_p1"]
+
+> coeffs_list_1d_size2 := map(coeff_name, posn_list_1d_size2, "coeff_");
+                coeffs_list_1d_size2 := ["coeff_0", "coeff_p1"]
 
 > 
-> print_name_list_dcl(data_var_list_1d_size2, "fp",
-> 		    "1d.cube.size2/data-var.dcl.c");
-> print_data_var_assign(posn_list_1d_size2, "data_",
-> 		      "1d.cube.size2/data-var.assign.c");
-> 
-> print_interp_coeff_var_store(posn_list_1d_size2, "", "coeff_I_",
-> 			     "1d.cube.size2/coeff-I.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size2, "factor", "coeff_dx_",
-> 			     "1d.cube.size2/coeff-dx.store.c");
-> 
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size2, "coeff_I_"), "fp",
-> 		    "1d.cube.size2/coeff-I.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size2, "coeff_dx_"), "fp",
-> 		    "1d.cube.size2/coeff-dx.dcl.c");
+> print_name_list_dcl(data_list_1d_size2, "fp",
+> 		    "1d.cube.size2/data-dcl.h");
+> print_name_list_dcl(coeffs_list_1d_size2, "fp",
+> 		    "1d.cube.size2/coeffs-dcl.h");
 > 
-> print_interp_cmpt__lc_of_data(posn_list_1d_size2,
-> 			      "result", "coeff_I_", "data_",
-> 			      "1d.cube.size2/interp-I.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size2,
-> 			      "result", "coeff_dx_", "data_",
-> 			      "1d.cube.size2/interp-dx.compute.c");
+> print_fetch_data(posn_list_1d_size2, "data->data_",
+> 		 "1d.cube.size2/fetch-data.c");
+> print_evaluate_molecule(posn_list_1d_size2,
+> 			"coeffs->coeff_", "data->data_",
+> 			"1d.cube.size2/evaluate-molecule.c");
+> print_store_coeffs(posn_list_1d_size2,
+> 		   "factor * coeffs->coeff_",
+> 		   "1d.cube.size2/store-coeffs.c");
 > 
 ################################################################################
 > 
@@ -1363,38 +1348,26 @@ posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 # generic stuff for 1d, cube, size=3
 #
 > 
-> data_var_list_1d_size3 := map(data_var_name, posn_list_1d_size3, "data_");
-           data_var_list_1d_size3 := ["data_m1", "data_0", "data_p1"]
+> data_list_1d_size3 := map(data_var_name, posn_list_1d_size3, "data_");
+             data_list_1d_size3 := ["data_m1", "data_0", "data_p1"]
+
+> coeffs_list_1d_size3 := map(coeff_name, posn_list_1d_size3, "coeff_");
+          coeffs_list_1d_size3 := ["coeff_m1", "coeff_0", "coeff_p1"]
 
 > 
-> print_name_list_dcl(data_var_list_1d_size3, "fp",
-> 		    "1d.cube.size3/data-var.dcl.c");
-> print_data_var_assign(posn_list_1d_size3, "data_",
-> 		      "1d.cube.size3/data-var.assign.c");
-> 
-> print_interp_coeff_var_store(posn_list_1d_size3, "", "coeff_I_",
-> 			     "1d.cube.size3/coeff-I.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size3, "factor", "coeff_dx_",
-> 			     "1d.cube.size3/coeff-dx.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size3, "factor", "coeff_dxx_",
-> 			     "1d.cube.size3/coeff-dxx.store.c");
-> 
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size3, "coeff_I_"), "fp",
-> 		    "1d.cube.size3/coeff-I.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size3, "coeff_dx_"), "fp",
-> 		    "1d.cube.size3/coeff-dx.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size3, "coeff_dxx_"), "fp",
-> 		    "1d.cube.size3/coeff-dxx.dcl.c");
+> print_name_list_dcl(data_list_1d_size3, "fp",
+> 		    "1d.cube.size3/data-dcl.h");
+> print_name_list_dcl(coeffs_list_1d_size3, "fp",
+> 		    "1d.cube.size3/coeffs-dcl.h");
 > 
-> print_interp_cmpt__lc_of_data(posn_list_1d_size3,
-> 			      "result", "coeff_I_", "data_",
-> 			      "1d.cube.size3/interp-I.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size3,
-> 			      "result", "coeff_dx_", "data_",
-> 			      "1d.cube.size3/interp-dx.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size3,
-> 			      "result", "coeff_dxx_", "data_",
-> 			      "1d.cube.size3/interp-dxx.compute.c");
+> print_fetch_data(posn_list_1d_size3, "data->data_",
+> 		 "1d.cube.size3/fetch-data.c");
+> print_evaluate_molecule(posn_list_1d_size3,
+> 			"coeffs->coeff_", "data->data_",
+> 			"1d.cube.size3/evaluate-molecule.c");
+> print_store_coeffs(posn_list_1d_size3,
+> 		   "factor * coeffs->coeff_",
+> 		   "1d.cube.size3/store-coeffs.c");
 > 
 ################################################################################
 > 
@@ -1402,38 +1375,26 @@ posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 # generic stuff for 1d, cube, size=4
 #
 > 
-> data_var_list_1d_size4 := map(data_var_name, posn_list_1d_size4, "data_");
-     data_var_list_1d_size4 := ["data_m1", "data_0", "data_p1", "data_p2"]
+> data_list_1d_size4 := map(data_var_name, posn_list_1d_size4, "data_");
+       data_list_1d_size4 := ["data_m1", "data_0", "data_p1", "data_p2"]
+
+> coeffs_list_1d_size4 := map(coeff_name, posn_list_1d_size4, "coeff_");
+    coeffs_list_1d_size4 := ["coeff_m1", "coeff_0", "coeff_p1", "coeff_p2"]
 
 > 
-> print_name_list_dcl(data_var_list_1d_size4, "fp",
-> 		    "1d.cube.size4/data-var.dcl.c");
-> print_data_var_assign(posn_list_1d_size4, "data_",
-> 		      "1d.cube.size4/data-var.assign.c");
-> 
-> print_interp_coeff_var_store(posn_list_1d_size4, "", "coeff_I_",
-> 			     "1d.cube.size4/coeff-I.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size4, "factor", "coeff_dx_",
-> 			     "1d.cube.size4/coeff-dx.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size4, "factor", "coeff_dxx_",
-> 			     "1d.cube.size4/coeff-dxx.store.c");
-> 
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size4, "coeff_I_"), "fp",
-> 		    "1d.cube.size4/coeff-I.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size4, "coeff_dx_"), "fp",
-> 		    "1d.cube.size4/coeff-dx.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size4, "coeff_dxx_"), "fp",
-> 		    "1d.cube.size4/coeff-dxx.dcl.c");
+> print_name_list_dcl(data_list_1d_size4, "fp",
+> 		    "1d.cube.size4/data-dcl.h");
+> print_name_list_dcl(coeffs_list_1d_size4, "fp",
+> 		    "1d.cube.size4/coeffs-dcl.h");
 > 
-> print_interp_cmpt__lc_of_data(posn_list_1d_size4,
-> 			      "result", "coeff_I_", "data_",
-> 			      "1d.cube.size4/interp-I.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size4,
-> 			      "result", "coeff_dx_", "data_",
-> 			      "1d.cube.size4/interp-dx.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size4,
-> 			      "result", "coeff_dxx_", "data_",
-> 			      "1d.cube.size4/interp-dxx.compute.c");
+> print_fetch_data(posn_list_1d_size4, "data->data_",
+> 		 "1d.cube.size4/fetch-data.c");
+> print_evaluate_molecule(posn_list_1d_size4,
+> 			"coeffs->coeff_", "data->data_",
+> 			"1d.cube.size4/evaluate-molecule.c");
+> print_store_coeffs(posn_list_1d_size4,
+> 		   "factor * coeffs->coeff_",
+> 		   "1d.cube.size4/store-coeffs.c");
 > 
 ################################################################################
 > 
@@ -1441,40 +1402,28 @@ posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 # generic stuff for 1d, cube, size=5
 #
 > 
-> data_var_list_1d_size5 := map(data_var_name, posn_list_1d_size5, "data_");
-data_var_list_1d_size5 :=
+> data_list_1d_size5 := map(data_var_name, posn_list_1d_size5, "data_");
+  data_list_1d_size5 := ["data_m2", "data_m1", "data_0", "data_p1", "data_p2"]
 
-    ["data_m2", "data_m1", "data_0", "data_p1", "data_p2"]
+> coeffs_list_1d_size5 := map(coeff_name, posn_list_1d_size5, "coeff_");
+coeffs_list_1d_size5 :=
+
+    ["coeff_m2", "coeff_m1", "coeff_0", "coeff_p1", "coeff_p2"]
 
 > 
-> print_name_list_dcl(data_var_list_1d_size5, "fp",
-> 		    "1d.cube.size5/data-var.dcl.c");
-> print_data_var_assign(posn_list_1d_size5, "data_",
-> 		      "1d.cube.size5/data-var.assign.c");
-> 
-> print_interp_coeff_var_store(posn_list_1d_size5, "", "coeff_I_",
-> 			     "1d.cube.size5/coeff-I.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size5, "factor", "coeff_dx_",
-> 			     "1d.cube.size5/coeff-dx.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size5, "factor", "coeff_dxx_",
-> 			     "1d.cube.size5/coeff-dxx.store.c");
-> 
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size5, "coeff_I_"), "fp",
-> 		    "1d.cube.size5/coeff-I.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size5, "coeff_dx_"), "fp",
-> 		    "1d.cube.size5/coeff-dx.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size5, "coeff_dxx_"), "fp",
-> 		    "1d.cube.size5/coeff-dxx.dcl.c");
+> print_name_list_dcl(data_list_1d_size5, "fp",
+> 		    "1d.cube.size5/data-dcl.h");
+> print_name_list_dcl(coeffs_list_1d_size5, "fp",
+> 		    "1d.cube.size5/coeffs-dcl.h");
 > 
-> print_interp_cmpt__lc_of_data(posn_list_1d_size5,
-> 			      "result", "coeff_I_", "data_",
-> 			      "1d.cube.size5/interp-I.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size5,
-> 			      "result", "coeff_dx_", "data_",
-> 			      "1d.cube.size5/interp-dx.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size5,
-> 			      "result", "coeff_dxx_", "data_",
-> 			      "1d.cube.size5/interp-dxx.compute.c");
+> print_fetch_data(posn_list_1d_size5, "data->data_",
+> 		 "1d.cube.size5/fetch-data.c");
+> print_evaluate_molecule(posn_list_1d_size5,
+> 			"coeffs->coeff_", "data->data_",
+> 			"1d.cube.size5/evaluate-molecule.c");
+> print_store_coeffs(posn_list_1d_size5,
+> 		   "factor * coeffs->coeff_",
+> 		   "1d.cube.size5/store-coeffs.c");
 > 
 ################################################################################
 > 
@@ -1482,40 +1431,30 @@ data_var_list_1d_size5 :=
 # generic stuff for 1d, cube, size=6
 #
 > 
-> data_var_list_1d_size6 := map(data_var_name, posn_list_1d_size6, "data_");
-data_var_list_1d_size6 :=
+> data_list_1d_size6 := map(data_var_name, posn_list_1d_size6, "data_");
+data_list_1d_size6 :=
 
     ["data_m2", "data_m1", "data_0", "data_p1", "data_p2", "data_p3"]
 
+> coeffs_list_1d_size6 := map(coeff_name, posn_list_1d_size6, "coeff_");
+coeffs_list_1d_size6 :=
+
+    ["coeff_m2", "coeff_m1", "coeff_0", "coeff_p1", "coeff_p2", "coeff_p3"]
+
 > 
-> print_name_list_dcl(data_var_list_1d_size6, "fp",
-> 		    "1d.cube.size6/data-var.dcl.c");
-> print_data_var_assign(posn_list_1d_size6, "data_",
-> 		      "1d.cube.size6/data-var.assign.c");
-> 
-> print_interp_coeff_var_store(posn_list_1d_size6, "", "coeff_I_",
-> 			     "1d.cube.size6/coeff-I.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size6, "factor", "coeff_dx_",
-> 			     "1d.cube.size6/coeff-dx.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size6, "factor", "coeff_dxx_",
-> 			     "1d.cube.size6/coeff-dxx.store.c");
-> 
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size6, "coeff_I_"), "fp",
-> 		    "1d.cube.size6/coeff-I.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size6, "coeff_dx_"), "fp",
-> 		    "1d.cube.size6/coeff-dx.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size6, "coeff_dxx_"), "fp",
-> 		    "1d.cube.size6/coeff-dxx.dcl.c");
-> 
-> print_interp_cmpt__lc_of_data(posn_list_1d_size6,
-> 			      "result", "coeff_I_", "data_",
-> 			      "1d.cube.size6/interp-I.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size6,
-> 			      "result", "coeff_dx_", "data_",
-> 			      "1d.cube.size6/interp-dx.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size6,
-> 			      "result", "coeff_dxx_", "data_",
-> 			      "1d.cube.size6/interp-dxx.compute.c");
+> print_name_list_dcl(data_list_1d_size6, "fp",
+> 		    "1d.cube.size6/data-dcl.h");
+> print_name_list_dcl(coeffs_list_1d_size6, "fp",
+> 		    "1d.cube.size6/coeffs-dcl.h");
+> 
+> print_fetch_data(posn_list_1d_size6, "data->data_",
+> 		 "1d.cube.size6/fetch-data.c");
+> print_evaluate_molecule(posn_list_1d_size6,
+> 			"coeffs->coeff_", "data->data_",
+> 			"1d.cube.size6/evaluate-molecule.c");
+> print_store_coeffs(posn_list_1d_size6,
+> 		   "factor * coeffs->coeff_",
+> 		   "1d.cube.size6/store-coeffs.c");
 > 
 ################################################################################
 > 
@@ -1523,41 +1462,31 @@ data_var_list_1d_size6 :=
 # generic stuff for 1d, cube, size=7
 #
 > 
-> data_var_list_1d_size7 := map(data_var_name, posn_list_1d_size7, "data_");
-data_var_list_1d_size7 := [
+> data_list_1d_size7 := map(data_var_name, posn_list_1d_size7, "data_");
+data_list_1d_size7 := [
 
     "data_m3", "data_m2", "data_m1", "data_0", "data_p1", "data_p2", "data_p3"]
 
+> coeffs_list_1d_size7 := map(coeff_name, posn_list_1d_size7, "coeff_");
+coeffs_list_1d_size7 := ["coeff_m3", "coeff_m2", "coeff_m1", "coeff_0",
+
+    "coeff_p1", "coeff_p2", "coeff_p3"]
+
+> 
+> print_name_list_dcl(data_list_1d_size7, "fp",
+> 		    "1d.cube.size7/data-dcl.h");
+> print_name_list_dcl(coeffs_list_1d_size7, "fp",
+> 		    "1d.cube.size7/coeffs-dcl.h");
 > 
-> print_name_list_dcl(data_var_list_1d_size7, "fp",
-> 		    "1d.cube.size7/data-var.dcl.c");
-> print_data_var_assign(posn_list_1d_size7, "data_",
-> 		      "1d.cube.size7/data-var.assign.c");
-> 
-> print_interp_coeff_var_store(posn_list_1d_size7, "", "coeff_I_",
-> 			     "1d.cube.size7/coeff-I.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size7, "factor", "coeff_dx_",
-> 			     "1d.cube.size7/coeff-dx.store.c");
-> print_interp_coeff_var_store(posn_list_1d_size7, "factor", "coeff_dxx_",
-> 			     "1d.cube.size7/coeff-dxx.store.c");
-> 
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size7, "coeff_I_"), "fp",
-> 		    "1d.cube.size7/coeff-I.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size7, "coeff_dx_"), "fp",
-> 		    "1d.cube.size7/coeff-dx.dcl.c");
-> print_name_list_dcl(map(coeff_name, posn_list_1d_size7, "coeff_dxx_"), "fp",
-> 		    "1d.cube.size7/coeff-dxx.dcl.c");
-> 
-> print_interp_cmpt__lc_of_data(posn_list_1d_size7,
-> 			      "result", "coeff_I_", "data_",
-> 			      "1d.cube.size7/interp-I.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size7,
-> 			      "result", "coeff_dx_", "data_",
-> 			      "1d.cube.size7/interp-dx.compute.c");
-> print_interp_cmpt__lc_of_data(posn_list_1d_size7,
-> 			      "result", "coeff_dxx_", "data_",
-> 			      "1d.cube.size7/interp-dxx.compute.c");
+> print_fetch_data(posn_list_1d_size7, "data->data_",
+> 		 "1d.cube.size7/fetch-data.c");
+> print_evaluate_molecule(posn_list_1d_size7,
+> 			"coeffs->coeff_", "data->data_",
+> 			"1d.cube.size7/evaluate-molecule.c");
+> print_store_coeffs(posn_list_1d_size7,
+> 		   "factor * coeffs->coeff_",
+> 		   "1d.cube.size7/store-coeffs.c");
 > 
 ################################################################################
 > quit
-bytes used=797948, alloc=720764, time=0.13
+bytes used=584788, alloc=589716, time=0.10