changeset 9756:b134960cea23

implement built-in tril/triu
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 23 Oct 2009 10:10:37 +0200
parents 4f4873f6f873
children 95ad9c2a27e2
files scripts/ChangeLog scripts/general/Makefile.in scripts/general/tril.m scripts/general/triu.m src/ChangeLog src/DLD-FUNCTIONS/tril.cc src/Makefile.in
diffstat 7 files changed, 436 insertions(+), 185 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Tue Oct 20 20:22:10 2009 +0200
+++ b/scripts/ChangeLog	Fri Oct 23 10:10:37 2009 +0200
@@ -1,3 +1,8 @@
+2009-10-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* general/tril.m, general/triu.m: Remove sources.
+	* general/Makefile.in: Update.
+
 2009-10-20  Soren Hauberg  <hauberg@gmail.com>
 
 	* general/interp2.m: improved error checking and support for bicubic
--- a/scripts/general/Makefile.in	Tue Oct 20 20:22:10 2009 +0200
+++ b/scripts/general/Makefile.in	Fri Oct 23 10:10:37 2009 +0200
@@ -46,7 +46,7 @@
   polyarea.m postpad.m prepad.m quadgk.m quadl.m quadv.m rat.m \
   rem.m repmat.m rot90.m rotdim.m runlength.m saveobj.m shift.m shiftdim.m \
   sortrows.m sph2cart.m strerror.m structfun.m subsindex.m \
-  triplequad.m trapz.m tril.m triu.m
+  triplequad.m trapz.m
 
 DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
 
--- a/scripts/general/tril.m	Tue Oct 20 20:22:10 2009 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,112 +0,0 @@
-## Copyright (C) 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2004, 2005,
-##               2006, 2007, 2008 John W. Eaton
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {} tril (@var{a}, @var{k})
-## @deftypefnx {Function File} {} triu (@var{a}, @var{k})
-## Return a new matrix formed by extracting the lower (@code{tril})
-## or upper (@code{triu}) triangular part of the matrix @var{a}, and
-## setting all other elements to zero.  The second argument is optional,
-## and specifies how many diagonals above or below the main diagonal should
-## also be set to zero.
-##
-## The default value of @var{k} is zero, so that @code{triu} and
-## @code{tril} normally include the main diagonal as part of the result
-## matrix.
-##
-## If the value of @var{k} is negative, additional elements above (for
-## @code{tril}) or below (for @code{triu}) the main diagonal are also
-## selected.
-##
-## The absolute value of @var{k} must not be greater than the number of
-## sub- or super-diagonals.
-##
-## For example,
-##
-## @example
-## @group
-## tril (ones (3), -1)
-##      @result{}  0  0  0
-##          1  0  0
-##          1  1  0
-## @end group
-## @end example
-##
-## @noindent
-## and
-##
-## @example
-## @group
-## tril (ones (3), 1)
-##      @result{}  1  1  0
-##          1  1  1
-##          1  1  1
-## @end group
-## @end example
-## @seealso{triu, diag}
-## @end deftypefn
-
-## Author: jwe
-
-function retval = tril (x, k)
-
-  if (nargin > 0)
-    if (isstruct (x))
-       error ("tril: structure arrays not supported");
-     endif 
-    [nr, nc] = size (x);
-  endif
-
-  if (nargin == 1)
-    k = 0;
-  elseif (nargin == 2)
-    if ((k > 0 && k > nc) || (k < 0 && k < -nr))
-      error ("tril: requested diagonal out of range");
-    endif
-  else
-    print_usage ();
-  endif
-
-  retval = resize (resize (x, 0), nr, nc);
-  for j = 1 : min (nc, nr+k)
-    nr_limit = max (1, j-k);
-    retval (nr_limit:nr, j) = x (nr_limit:nr, j);
-  endfor
-
-endfunction
-
-%!test
-%! a = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
-%! 
-%! l0 = [1, 0, 0; 4, 5, 0; 7, 8, 9; 10, 11, 12];
-%! l1 = [1, 2, 0; 4, 5, 6; 7, 8, 9; 10, 11, 12];
-%! l2 = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
-%! lm1 = [0, 0, 0; 4, 0, 0; 7, 8, 0; 10, 11, 12];
-%! lm2 = [0, 0, 0; 0, 0, 0; 7, 0, 0; 10, 11, 0];
-%! lm3 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 10, 0, 0];
-%! lm4 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 0, 0, 0];
-%! 
-%! assert((tril (a, -4) == lm4 && tril (a, -3) == lm3
-%! && tril (a, -2) == lm2 && tril (a, -1) == lm1
-%! && tril (a) == l0 && tril (a, 1) == l1 && tril (a, 2) == l2));
-
-%!error tril ();
-
-%!error tril (1, 2, 3);
-
--- a/scripts/general/triu.m	Tue Oct 20 20:22:10 2009 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,71 +0,0 @@
-## Copyright (C) 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2005, 2006,
-##               2007, 2008 John W. Eaton
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {} triu (@var{a}, @var{k})
-## See tril.
-## @end deftypefn
-
-## Author: jwe
-
-function retval = triu (x, k)
-
-  if (nargin > 0)
-    if (isstruct (x))
-       error ("tril: structure arrays not supported");
-     endif 
-    [nr, nc] = size (x);
-  endif
-  if (nargin == 1)
-    k = 0;
-  elseif (nargin == 2)
-    if ((k > 0 && k > nc) || (k < 0 && k < -nr))
-      error ("triu: requested diagonal out of range");
-    endif
-  else
-    print_usage ();
-  endif
-
-  retval = resize (resize (x, 0), nr, nc);
-  for j = max (1, k+1) : nc
-    nr_limit = min (nr, j-k);
-    retval (1:nr_limit, j) = x (1:nr_limit, j);
-  endfor
-
-endfunction
-
-%!test
-%! a = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
-%! 
-%! u0 = [1, 2, 3; 0, 5, 6; 0, 0, 9; 0, 0, 0];
-%! u1 = [0, 2, 3; 0, 0, 6; 0, 0, 0; 0, 0, 0];
-%! u2 = [0, 0, 3; 0, 0, 0; 0, 0, 0; 0, 0, 0];
-%! u3 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 0, 0, 0];
-%! um1 = [1, 2, 3; 4, 5, 6; 0, 8, 9; 0, 0, 12];
-%! um2 = [1, 2, 3; 4, 5, 6; 7, 8, 9; 0, 11, 12];
-%! um3 = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
-%! 
-%! assert((triu (a, -3) == um3 && triu (a, -2) == um2
-%! && triu (a, -1) == um1 && triu (a) == u0 && triu (a, 1) == u1
-%! && triu (a, 2) == u2 && triu (a, 3) == u3));
-
-%!error triu ();
-
-%!error triu (1, 2, 3);
-
--- a/src/ChangeLog	Tue Oct 20 20:22:10 2009 +0200
+++ b/src/ChangeLog	Fri Oct 23 10:10:37 2009 +0200
@@ -1,3 +1,8 @@
+2009-10-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/tril.cc: New source.
+	* Makefile.in: Include it.
+
 2009-10-22  Jaroslav Hajek  <highegg@gmail.com>
 
 	* error.cc (verror (bool, std::ostream&, ..., bool)): Add optional
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/tril.cc	Fri Oct 23 10:10:37 2009 +0200
@@ -0,0 +1,424 @@
+/*
+
+Copyright (C) 2004, 2007 David Bateman
+Copyright (C) 2009 VZLU Prague
+
+This program is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3, or (at your option)
+any later version.
+
+This program is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <algorithm>
+#include "Array.h"
+#include "Sparse.h"
+#include "mx-base.h"
+
+#include "ov.h"
+#include "Cell.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "oct-obj.h"
+
+// The bulk of the work.
+template <class T>
+static Array<T> 
+do_tril (const Array<T>& a, octave_idx_type k, bool pack)
+{
+  octave_idx_type nr = a.rows (), nc = a.columns ();
+  const T *avec = a.fortran_vec ();
+
+  if (pack)
+    {
+      octave_idx_type j1 = std::min (std::max (0, k), nc); 
+      octave_idx_type j2 = std::min (std::max (0, nr + k), nc);
+      octave_idx_type n = j1 * nr + ((j2 - j1) * (nr-(j1-k) + nr-(j2-1-k))) / 2;
+      Array<T> r (n);
+      T *rvec = r.fortran_vec ();
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type ii = std::min (std::max (0, j - k), nr);
+          rvec = std::copy (avec + ii, avec + nr, rvec);
+          avec += nr;
+        }
+
+      return r;
+    }
+  else
+    {
+      Array<T> r (a.dims ());
+      T *rvec = r.fortran_vec ();
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type ii = std::min (std::max (0, j - k), nr);
+          std::fill (rvec, rvec + ii, T());
+          std::copy (avec + ii, avec + nr, rvec + ii);
+          avec += nr;
+          rvec += nr;
+        }
+
+      return r;
+    }
+}
+
+template <class T>
+static Array<T> 
+do_triu (const Array<T>& a, octave_idx_type k, bool pack)
+{
+  octave_idx_type nr = a.rows (), nc = a.columns ();
+  const T *avec = a.fortran_vec ();
+
+  if (pack)
+    {
+      octave_idx_type j1 = std::min (std::max (0, k), nc); 
+      octave_idx_type j2 = std::min (std::max (0, nr + k), nc);
+      octave_idx_type n = ((j2 - j1) * ((j1+1-k) + (j2-k))) / 2 + (nc - j2) * nr;
+      Array<T> r (n);
+      T *rvec = r.fortran_vec ();
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type ii = std::min (std::max (0, j + 1 - k), nr);
+          rvec = std::copy (avec, avec + ii, rvec);
+          avec += nr;
+        }
+
+      return r;
+    }
+  else
+    {
+      NoAlias<Array<T> > r (a.dims ());
+      T *rvec = r.fortran_vec ();
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type ii = std::min (std::max (0, j + 1 - k), nr);
+          std::copy (avec, avec + ii, rvec);
+          std::fill (rvec + ii, rvec + nr, T());
+          avec += nr;
+          rvec += nr;
+        }
+
+      return r;
+    }
+}
+
+// These two are by David Bateman.
+// FIXME: optimizations possible. "pack" support missing.
+
+template <class T>
+static Sparse<T> 
+do_tril (const Sparse<T>& a, octave_idx_type k, bool pack)
+{
+  if (pack) // FIXME
+    {
+      error ("tril: \"pack\" not implemented for sparse matrices");
+      return Sparse<T> ();
+    }
+
+  Sparse<T> m = a;
+  octave_idx_type nc = m.cols();
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
+      if (m.ridx(i) < j-k)
+        m.data(i) = 0.;
+
+  m.maybe_compress (true);
+  return m;
+}
+
+template <class T>
+static Sparse<T> 
+do_triu (const Sparse<T>& a, octave_idx_type k, bool pack)
+{
+  if (pack) // FIXME
+    {
+      error ("triu: \"pack\" not implemented for sparse matrices");
+      return Sparse<T> ();
+    }
+
+  Sparse<T> m = a;
+  octave_idx_type nc = m.cols();
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
+      if (m.ridx(i) > j-k)
+        m.data(i) = 0.;
+
+  m.maybe_compress (true);
+  return m;
+}
+
+// Convenience dispatchers.
+template <class T>
+static Array<T>
+do_trilu (const Array<T>& a, octave_idx_type k, bool lower, bool pack)
+{
+  return lower ? do_tril (a, k, pack) : do_triu (a, k, pack);
+}
+
+template <class T>
+static Sparse<T>
+do_trilu (const Sparse<T>& a, octave_idx_type k, bool lower, bool pack)
+{
+  return lower ? do_tril (a, k, pack) : do_triu (a, k, pack);
+}
+
+static octave_value
+do_trilu (const std::string& name,
+          const octave_value_list& args)
+{
+  bool lower = name == "tril";
+
+  octave_value retval;
+  int nargin = args.length ();
+  octave_idx_type k = 0;
+  bool pack = false;
+  if (nargin >= 2 && args(nargin-1).is_string ())
+    {
+      pack = args(nargin-1).string_value () == "pack";
+      nargin--;
+    }
+
+  if (nargin == 2)
+    {
+      k = args(1).int_value (true);
+      
+      if (error_state)
+	return retval;
+    }
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+  else
+    {
+      octave_value arg = args (0);
+
+      dim_vector dims = arg.dims ();
+      if (dims.length () != 2)
+        error ("%s: needs a 2D matrix", name.c_str ());
+      else if (k < -dims (0) || k > dims(1))
+        error ("%s: requested diagonal out of range", name.c_str ());
+      else
+        {
+          switch (arg.builtin_type ())
+            {
+            case btyp_double:
+              if (arg.is_sparse_type ())
+                retval = do_trilu (arg.sparse_matrix_value (), k, lower, pack);
+              else
+                retval = do_trilu (arg.array_value (), k, lower, pack);
+              break;
+            case btyp_complex:
+              if (arg.is_sparse_type ())
+                retval = do_trilu (arg.sparse_complex_matrix_value (), k, lower, pack);
+              else
+                retval = do_trilu (arg.complex_array_value (), k, lower, pack);
+              break;
+            case btyp_bool:
+              if (arg.is_sparse_type ())
+                retval = do_trilu (arg.sparse_bool_matrix_value (), k, lower, pack);
+              else
+                retval = do_trilu (arg.bool_array_value (), k, lower, pack);
+              break;
+#define ARRAYCASE(TYP) \
+            case btyp_ ## TYP: \
+              retval = do_trilu (arg.TYP ## _array_value (), k, lower, pack); \
+              break
+            ARRAYCASE (float);
+            ARRAYCASE (float_complex);
+            ARRAYCASE (int8);
+            ARRAYCASE (int16);
+            ARRAYCASE (int32);
+            ARRAYCASE (int64);
+            ARRAYCASE (uint8);
+            ARRAYCASE (uint16);
+            ARRAYCASE (uint32);
+            ARRAYCASE (uint64);
+            ARRAYCASE (char);
+#undef ARRAYCASE
+            default:
+              {
+                // Generic code that works on octave-values, that is slow
+                // but will also work on arbitrary user types
+
+                if (pack) // FIXME
+                  {
+                    error ("%s: \"pack\" not implemented for class %s", 
+                           name.c_str (), arg.class_name ().c_str ());
+                    return octave_value ();
+                  }
+
+                octave_value tmp = arg;
+                if (arg.numel () == 0)
+                  return arg;
+
+                octave_idx_type nr = dims(0), nc = dims (1);
+
+                // The sole purpose of the below is to force the correct
+                // matrix size. This would not be necessary if the
+                // octave_value resize function allowed a fill_value.
+                // It also allows odd attributes in some user types
+                // to be handled. With a fill_value ot should be replaced 
+                // with
+                //
+                // octave_value_list ov_idx;
+                // tmp = tmp.resize(dim_vector (0,0)).resize (dims, fill_value);
+
+                octave_value_list ov_idx;
+                std::list<octave_value_list> idx_tmp;
+                ov_idx(1) = static_cast<double> (nc+1);
+                ov_idx(0) = Range (1, nr);
+                idx_tmp.push_back (ov_idx);
+                ov_idx(1) = static_cast<double> (nc);
+                tmp = tmp.resize (dim_vector (0,0));
+                tmp = tmp.subsasgn("(",idx_tmp, arg.do_index_op (ov_idx));
+                tmp = tmp.resize(dims);
+
+                if (lower)
+                  {
+                    octave_idx_type st = nc < nr + k ? nc : nr + k;
+
+                    for (octave_idx_type j = 1; j <= st; j++)
+                      {
+                        octave_idx_type nr_limit = 1 > j - k ? 1 : j - k;
+                        ov_idx(1) = static_cast<double> (j);
+                        ov_idx(0) = Range (nr_limit, nr);
+                        std::list<octave_value_list> idx;
+                        idx.push_back (ov_idx);
+
+                        tmp = tmp.subsasgn ("(", idx, arg.do_index_op(ov_idx));
+
+                        if (error_state)
+                          return retval;
+                      }
+                  }
+                else
+                  {
+                    octave_idx_type st = k + 1 > 1 ? k + 1 : 1;
+
+                    for (octave_idx_type j = st; j <= nc; j++)
+                      {
+                        octave_idx_type nr_limit = nr < j - k ? nr : j - k;
+                        ov_idx(1) = static_cast<double> (j);
+                        ov_idx(0) = Range (1, nr_limit);
+                        std::list<octave_value_list> idx;
+                        idx.push_back (ov_idx);
+
+                        tmp = tmp.subsasgn ("(", idx, arg.do_index_op(ov_idx));
+
+                        if (error_state)
+                          return retval;
+                      }
+                  }
+
+                retval = tmp;
+              }
+            }
+        }
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (tril, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Function File} {} tril (@var{a}, @var{k}, @var{pack})\n\
+@deftypefnx {Function File} {} triu (@var{a}, @var{k}, @var{pack})\n\
+Return a new matrix formed by extracting extract the lower (@code{tril})\n\
+or upper (@code{triu}) triangular part of the matrix @var{a}, and\n\
+setting all other elements to zero.  The second argument is optional,\n\
+and specifies how many diagonals above or below the main diagonal should\n\
+also be set to zero.\n\
+\n\
+The default value of @var{k} is zero, so that @code{triu} and\n\
+@code{tril} normally include the main diagonal as part of the result\n\
+matrix.\n\
+\n\
+If the value of @var{k} is negative, additional elements above (for\n\
+@code{tril}) or below (for @code{triu}) the main diagonal are also\n\
+selected.\n\
+\n\
+The absolute value of @var{k} must not be greater than the number of\n\
+sub- or super-diagonals.\n\
+\n\
+For example,\n\
+\n\
+@example\n\
+@group\n\
+tril (ones (3), -1)\n\
+     @result{}  0  0  0\n\
+         1  0  0\n\
+         1  1  0\n\
+@end group\n\
+@end example\n\
+\n\
+@noindent\n\
+and\n\
+\n\
+@example\n\
+@group\n\
+tril (ones (3), 1)\n\
+     @result{}  1  1  0\n\
+         1  1  1\n\
+         1  1  1\n\
+@end group\n\
+@end example\n\
+\n\
+If the option \"pack\" is given as third argument, the extracted elements\n\
+are not inserted into a matrix, but rather stacked column-wise one above other.\n\
+@seealso{triu, diag}\n\
+@end deftypefn")
+{
+  return do_trilu ("tril", args);
+}
+
+DEFUN_DLD (triu, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Function File} {} triu (@var{a}, @var{k})\n\
+See tril.\n\
+@end deftypefn")
+{
+  return do_trilu ("triu", args);
+}
+
+/*
+
+%!test
+%! a = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
+%! 
+%! l0 = [1, 0, 0; 4, 5, 0; 7, 8, 9; 10, 11, 12];
+%! l1 = [1, 2, 0; 4, 5, 6; 7, 8, 9; 10, 11, 12];
+%! l2 = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
+%! lm1 = [0, 0, 0; 4, 0, 0; 7, 8, 0; 10, 11, 12];
+%! lm2 = [0, 0, 0; 0, 0, 0; 7, 0, 0; 10, 11, 0];
+%! lm3 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 10, 0, 0];
+%! lm4 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 0, 0, 0];
+%! 
+%! assert((tril (a, -4) == lm4 && tril (a, -3) == lm3
+%! && tril (a, -2) == lm2 && tril (a, -1) == lm1
+%! && tril (a) == l0 && tril (a, 1) == l1 && tril (a, 2) == l2));
+
+%!error tril ();
+
+*/
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/src/Makefile.in	Tue Oct 20 20:22:10 2009 +0200
+++ b/src/Makefile.in	Fri Oct 23 10:10:37 2009 +0200
@@ -73,7 +73,7 @@
 	lsode.cc lu.cc luinc.cc matrix_type.cc max.cc md5sum.cc \
 	pinv.cc qr.cc quad.cc qz.cc rand.cc rcond.cc regexp.cc \
 	schur.cc sparse.cc spparms.cc sqrtm.cc sub2ind.cc svd.cc syl.cc \
-	symrcm.cc symbfact.cc time.cc tsearch.cc typecast.cc \
+	symrcm.cc symbfact.cc time.cc tril.cc tsearch.cc typecast.cc \
 	urlwrite.cc __contourc__.cc __delaunayn__.cc \
 	__dsearchn__.cc __glpk__.cc __lin_interpn__.cc \
 	__magick_read__.cc __pchip_deriv__.cc __qp__.cc \