Make duneevaluate work with row vectors of input

src/duneevaluate.cc

@@ -11,12 +11,11 @@
 
 #include <cassert>
 
-DEFUN_DLD(duneevaluate, args, nargout, "duneevaluate(y)\n\
+DEFUN_DLD(duneevaluate, args, nargout, "duneevaluate(A)\n\
 \n\
-Evaluate x -> 1/2<Ax,x> - <b,x> + H(|x|) at y using DUNE\n") {
+Evaluate x -> 1/2<Ax,x> - <b,x> + H(|x|) at each point y that is a column vector of A using DUNE\n") {
   assert(args.length() == 1);
   assert(nargout <= 1);
-  ColumnVector point(args(0).vector_value());
 
   int const dim = 2;
 
@@ -35,7 +34,13 @@ Evaluate x -> 1/2<Ax,x> - <b,x> + H(|x|) at y using DUNE\n") {
   Functional J(A, b, Dune::MyNonlinearity<dim>(f));
   // }}}
 
+  Matrix points(args(0).matrix_value());
+  size_t const columns = points.columns();
+  RowVector ret(columns);
   Functional::SmallVector our_point;
-  Dune::octaveToDune<dim>(point, our_point);
-  return octave_value(J(our_point));
+  for (size_t i = 0; i < columns; ++i) {
+    Dune::octaveToDune<dim>(points.column(i), our_point);
+    ret(i) = J(our_point);
+  }
+  return octave_value(ret);
 }

src/foo.m

@@ -11,12 +11,16 @@ y = -125:1:50;
 xlabel('x')
 ylabel('y')
 
+tic
 for i=1:length(y)
+  a = zeros(2,length(x));
   for j=1:length(x)
-    f(i,j) = duneevaluate([ X(i,j); Y(i,j) ]);
+    a(:,j) = [ X(i,j); Y(i,j) ];
   end
+  f(i,:) = duneevaluate(a);
 end
 clear X Y;
+toc
 
 surf(x, y, f)
 hold on;