[pure-lang-svn] SF.net SVN: pure-lang:[769] pure/trunk
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
[pure-lang-svn] SF.net SVN: pure-lang:[769] pure/trunk
|
From: <ag...@us...> - 2008-09-16 12:07:39
|
Revision: 769
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=769&view=rev
Author: agraef
Date: 2008-09-16 12:07:48 +0000 (Tue, 16 Sep 2008)
Log Message:
-----------
Add basic GSL matrix infrastructure (not quite finished yet).
Modified Paths:
--------------
pure/trunk/expr.cc
pure/trunk/expr.hh
pure/trunk/interpreter.cc
pure/trunk/interpreter.hh
pure/trunk/lib/prelude.pure
pure/trunk/printer.cc
pure/trunk/runtime.cc
pure/trunk/runtime.h
pure/trunk/symtable.cc
pure/trunk/symtable.hh
Modified: pure/trunk/expr.cc
===================================================================
--- pure/trunk/expr.cc 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/expr.cc 2008-09-16 12:07:48 UTC (rev 769)
@@ -29,6 +29,9 @@
if (data.x[1]) data.x[1]->del();
if (data.x[2]) data.x[2]->del();
break;
+ case MATRIX:
+ if (data.xs) delete data.xs;
+ break;
case CASE:
if (data.c.x) data.c.x->del();
if (data.c.r) delete data.c.r;
Modified: pure/trunk/expr.hh
===================================================================
--- pure/trunk/expr.hh 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/expr.hh 2008-09-16 12:07:48 UTC (rev 769)
@@ -72,6 +72,21 @@
size_t len() const { return size; }
};
+/* Smart expression pointers (see below for the full definition). These are
+ used recursively as components in matrix representations and rule lists in
+ the expression data structure. */
+
+class expr;
+
+/* Expression lists and lists of those. These are used to represent
+ collections of expressions and generic matrix data in a structured way
+ which facilitates code generation. In the case of exprll, each list member
+ represents a matrix "row" which is in turn described by a list of
+ "columns". */
+
+typedef list<expr> exprl;
+typedef list<exprl> exprll;
+
/* Rule lists are used to encode collections of equations and other rule sets
in 'case' expressions and the like. See the definition of the rule struct
at the end of this header. */
@@ -113,6 +128,16 @@
CASE = -10, // case expression
WHEN = -11, // when expression
WITH = -12, // with expression
+ // GSL matrix types:
+ MATRIX = -32, // generic GSL matrix, double matrix in runtime exprs
+ CMATRIX = -31, // complex matrix in runtime exprs
+ IMATRIX = -30, // integer matrix in runtime exprs
+ /* Other values in the range -17..-29 reserved for later use in the
+ runtime expression data structure. Note that all GSL-related tags,
+ taken as an unsigned binary quantity, are of the form 0xffffffe0+t,
+ where the least significant nibble t=0x0..0xf denotes corresponding
+ subtypes in runtime matrix data. For compile time expressions only the
+ EXPR::MATRIX tag (t=0) is used. */
};
// special flag values used during compilation:
@@ -137,6 +162,7 @@
uint8_t idx; // de Bruin index
} v;
EXPR *x[3]; // APP, LAMBDA, COND
+ exprll *xs; // MATRIX
struct { // CASE, WHEN, WITH
EXPR *x; // expression
union {
@@ -204,6 +230,9 @@
refc(0), tag(_tag), m(0), flags(0), ttag(0), astag(0), aspath(0)
{ assert(_tag == WITH);
data.c.x = newref(_arg); data.c.e = _e; }
+ EXPR(int32_t _tag, exprll *_args) :
+ refc(0), tag(_tag), m(0), flags(0), ttag(0), astag(0), aspath(0)
+ { assert(_tag == MATRIX); data.xs = _args; }
EXPR(EXPR *_fun, EXPR *_arg) :
refc(0), tag(APP), m(0), flags(0), ttag(0), astag(0), aspath(0)
{ data.x[0] = newref(_fun); data.x[1] = newref(_arg); }
@@ -222,9 +251,6 @@
/* Smart expression pointers. These take care of reference counting
automagically. */
-class expr;
-typedef list<expr> exprl;
-
class expr {
EXPR* p;
// debug helper
@@ -278,6 +304,8 @@
p(new EXPR(tag, &*arg, rules)) { p->incref(); }
expr(int32_t tag, expr arg, env *e) :
p(new EXPR(tag, &*arg, e)) { p->incref(); }
+ expr(int32_t tag, exprll *args) :
+ p(new EXPR(tag, args)) { p->incref(); }
expr(expr fun, expr arg) :
p(new EXPR(&*fun, &*arg)) { p->incref(); }
expr(expr fun, expr arg1, expr arg2) :
@@ -337,6 +365,8 @@
p->tag == EXPR::WHEN ||
p->tag == EXPR::WITH);
return expr(p->data.c.x); }
+ exprll *xvals() const { assert(p->tag == EXPR::MATRIX);
+ return p->data.xs; }
rulel *rules() const { assert(p->tag == EXPR::CASE ||
p->tag == EXPR::WHEN);
return p->data.c.r; }
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/interpreter.cc 2008-09-16 12:07:48 UTC (rev 769)
@@ -16,6 +16,7 @@
#ifdef HAVE_GSL
#include <gsl/gsl_version.h>
+#include <gsl/gsl_matrix.h>
#endif
uint8_t interpreter::g_verbose = 0;
@@ -246,6 +247,11 @@
declare_extern((void*)pure_apply,
"pure_apply", "expr*", 2, "expr*", "expr*");
+ declare_extern((void*)pure_matrix_rows,
+ "pure_matrix_rows", "expr*", -1, "int");
+ declare_extern((void*)pure_matrix_columns,
+ "pure_matrix_columns", "expr*", -1, "int");
+
declare_extern((void*)pure_listl,
"pure_listl", "expr*", -1, "int");
declare_extern((void*)pure_tuplel,
@@ -862,6 +868,68 @@
// Only null pointer constants permitted right now.
throw err("pointer must be null in constant definition");
return expr(EXPR::PTR, x->data.p);
+ case EXPR::MATRIX: {
+#ifdef HAVE_GSL
+ if (x->data.mat && x->data.mat->p) {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat->p;
+ exprll *xs = new exprll;
+ for (size_t i = 0; i < m->size1; i++) {
+ xs->push_back(exprl());
+ exprl& ys = xs->back();
+ for (size_t j = 0; j < m->size2; j++) {
+ ys.push_back(expr(EXPR::DBL, m->data[i * m->tda + j]));
+ }
+ }
+ return expr(EXPR::MATRIX, m);
+ } else
+ return expr(EXPR::MATRIX, new exprll);
+#else
+ throw err("GSL matrices not supported in this implementation");
+ return expr(EXPR::MATRIX, new exprll);
+#endif
+ }
+ case EXPR::IMATRIX: {
+#ifdef HAVE_GSL
+ if (x->data.mat && x->data.mat->p) {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat->p;
+ exprll *xs = new exprll;
+ for (size_t i = 0; i < m->size1; i++) {
+ xs->push_back(exprl());
+ exprl& ys = xs->back();
+ for (size_t j = 0; j < m->size2; j++) {
+ ys.push_back(expr(EXPR::INT, m->data[i * m->tda + j]));
+ }
+ }
+ return expr(EXPR::MATRIX, m);
+ } else
+ return expr(EXPR::MATRIX, new exprll);
+#else
+ throw err("GSL matrices not supported in this implementation");
+ return expr(EXPR::MATRIX, new exprll);
+#endif
+ }
+ case EXPR::CMATRIX: {
+#ifdef HAVE_GSL
+ if (x->data.mat && x->data.mat->p) {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat->p;
+ exprll *xs = new exprll;
+ for (size_t i = 0; i < m->size1; i++) {
+ xs->push_back(exprl());
+ exprl& ys = xs->back();
+ for (size_t j = 0; j < m->size2; j++) {
+ expr u = expr(EXPR::DBL, m->data[2*(i * m->tda + j)]);
+ expr v = expr(EXPR::DBL, m->data[2*(i * m->tda + j) + 1]);
+ ys.push_back(expr(symtab.complex_rect_sym().x, u, v));
+ }
+ }
+ return expr(EXPR::MATRIX, m);
+ } else
+ return expr(EXPR::MATRIX, new exprll);
+#else
+ throw err("GSL matrices not supported in this implementation");
+ return expr(EXPR::MATRIX, new exprll);
+#endif
+ }
default:
assert(x->tag > 0);
if (x->data.clos && x->data.clos->local)
@@ -1126,6 +1194,14 @@
{
if (x.is_null()) return;
switch (x.tag()) {
+ case EXPR::MATRIX:
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++)
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ compile(*ys);
+ }
+ break;
case EXPR::APP:
compile(x.xval1());
compile(x.xval2());
@@ -1578,6 +1654,9 @@
break;
}
// these must not occur on the lhs:
+ case EXPR::MATRIX:
+ throw err("matrix expression not permitted in pattern");
+ break;
case EXPR::LAMBDA:
throw err("lambda expression not permitted in pattern");
break;
@@ -1717,6 +1796,20 @@
case EXPR::STR:
case EXPR::PTR:
return x;
+ // matrix:
+ case EXPR::MATRIX: {
+ exprll *us = new exprll;
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++) {
+ us->push_back(exprl());
+ exprl& vs = us->back();
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ vs.push_back(subst(vars, *ys, idx));
+ }
+ }
+ return expr(EXPR::MATRIX, us);
+ }
// application:
case EXPR::APP:
if (x.xval1().tag() == symtab.amp_sym().f) {
@@ -1827,6 +1920,20 @@
case EXPR::STR:
case EXPR::PTR:
return x;
+ // matrix:
+ case EXPR::MATRIX: {
+ exprll *us = new exprll;
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++) {
+ us->push_back(exprl());
+ exprl& vs = us->back();
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ vs.push_back(fsubst(funs, *ys, idx));
+ }
+ }
+ return expr(EXPR::MATRIX, us);
+ }
// application:
case EXPR::APP:
if (x.xval1().tag() == symtab.amp_sym().f) {
@@ -1929,6 +2036,20 @@
case EXPR::STR:
case EXPR::PTR:
return x;
+ // matrix:
+ case EXPR::MATRIX: {
+ exprll *us = new exprll;
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++) {
+ us->push_back(exprl());
+ exprl& vs = us->back();
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ vs.push_back(csubst(*ys));
+ }
+ }
+ return expr(EXPR::MATRIX, us);
+ }
// application:
case EXPR::APP:
if (x.xval1().tag() == symtab.amp_sym().f) {
@@ -2044,6 +2165,20 @@
case EXPR::STR:
case EXPR::PTR:
return x;
+ // matrix:
+ case EXPR::MATRIX: {
+ exprll *us = new exprll;
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++) {
+ us->push_back(exprl());
+ exprl& vs = us->back();
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ vs.push_back(macsubst(*ys));
+ }
+ }
+ return expr(EXPR::MATRIX, us);
+ }
// application:
case EXPR::APP: {
expr u = macsubst(x.xval1()),
@@ -2132,6 +2267,20 @@
case EXPR::STR:
case EXPR::PTR:
return x;
+ // matrix:
+ case EXPR::MATRIX: {
+ exprll *us = new exprll;
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++) {
+ us->push_back(exprl());
+ exprl& vs = us->back();
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ vs.push_back(varsubst(*ys, offs));
+ }
+ }
+ return expr(EXPR::MATRIX, us);
+ }
// application:
case EXPR::APP: {
expr u = varsubst(x.xval1(), offs),
@@ -2226,6 +2375,20 @@
return v;
} else
return y;
+ // matrix:
+ case EXPR::MATRIX: {
+ exprll *us = new exprll;
+ for (exprll::iterator xs = y.xvals()->begin(), end = y.xvals()->end();
+ xs != end; xs++) {
+ us->push_back(exprl());
+ exprl& vs = us->back();
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ vs.push_back(macred(x, *ys, idx));
+ }
+ }
+ return expr(EXPR::MATRIX, us);
+ }
// application:
case EXPR::APP:
if (y.xval1().tag() == symtab.amp_sym().f) {
@@ -3034,6 +3197,14 @@
}
}
break;
+ case EXPR::MATRIX:
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++)
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++) {
+ build_map(*ys);
+ }
+ break;
case EXPR::APP: {
expr f; uint32_t n = count_args(x, f);
interpreter& interp = *interpreter::g_interp;
@@ -3919,6 +4090,8 @@
return pure_string_dup(x.sval());
case EXPR::PTR:
return pure_pointer(x.pval());
+ case EXPR::MATRIX:
+ return const_matrix_value(x);
case EXPR::APP:
return const_app_value(x);
default: {
@@ -3959,6 +4132,36 @@
}
}
+pure_expr *interpreter::const_matrix_value(expr x)
+{
+ size_t n = x.xvals()->size(), m = 0, i = 0, j = 0;
+ pure_expr **us = new pure_expr*[n], **vs = 0, *ret;
+ assert(us);
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++, i++) {
+ m = xs->size(); j = 0; vs = new pure_expr*[m];
+ assert(vs);
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++, j++) {
+ vs[j] = const_value(*ys);
+ if (!vs[j]) goto err;
+ }
+ us[i] = pure_matrix_columnsv(m, vs);
+ if (!us[i]) goto err;
+ delete[] vs;
+ }
+ ret = pure_matrix_rowsv(n, us);
+ delete[] us;
+ return ret;
+ err:
+ // bail out
+ for (size_t k = 0; k < j; k++) pure_freenew(vs[k]);
+ if (vs) delete[] vs;
+ for (size_t k = 0; k < i; k++) pure_freenew(us[k]);
+ if (us) delete[] us;
+ return 0;
+}
+
pure_expr *interpreter::const_app_value(expr x)
{
if (x.tag() == EXPR::APP) {
@@ -4802,6 +5005,25 @@
// FIXME: Only null pointers are supported right now.
assert(x.pval() == 0);
return pbox(x.pval());
+ // matrix:
+ case EXPR::MATRIX: {
+ size_t n = x.xvals()->size(), i = 1;
+ vector<Value*> us(n+1);
+ us[0] = UInt(n);
+ for (exprll::iterator xs = x.xvals()->begin(), end = x.xvals()->end();
+ xs != end; xs++, i++) {
+ size_t m = xs->size(), j = 1;
+ vector<Value*> vs(m+1);
+ vs[0] = UInt(m);
+ for (exprl::iterator ys = xs->begin(), end = xs->end();
+ ys != end; ys++, j++) {
+ vs[j] = codegen(*ys);
+ }
+ us[i] =
+ act_env().CreateCall(module->getFunction("pure_matrix_columns"), vs);
+ }
+ return act_env().CreateCall(module->getFunction("pure_matrix_rows"), us);
+ }
// application:
case EXPR::APP:
if (x.ttag() != 0) {
@@ -5928,6 +6150,9 @@
case EXPR::PTR:
assert(0 && "not implemented");
break;
+ case EXPR::MATRIX:
+ assert(0 && "not implemented");
+ break;
case EXPR::APP: {
// first match the tag...
BasicBlock *ok1bb = BasicBlock::Create("arg1");
Modified: pure/trunk/interpreter.hh
===================================================================
--- pure/trunk/interpreter.hh 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/interpreter.hh 2008-09-16 12:07:48 UTC (rev 769)
@@ -520,6 +520,7 @@
Env& act_env() { assert(!envstk.empty()); return *envstk.front(); }
Builder& act_builder() { return act_env().builder; }
pure_expr *const_value(expr x);
+ pure_expr *const_matrix_value(expr x);
pure_expr *const_app_value(expr x);
expr pure_expr_to_expr(pure_expr *x);
pure_expr *doeval(expr x, pure_expr*& e);
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/lib/prelude.pure 2008-09-16 12:07:48 UTC (rev 769)
@@ -27,6 +27,8 @@
nullary failed_cond; // failed conditional (guard, if-then-else)
nullary failed_match; // failed pattern match (lambda, case, etc.)
nullary stack_fault; // not enough stack space (PURE_STACK limit)
+nullary not_implemented; // operation not implemented
+// bad_matrix_value x; // error in matrix construction
/* Other exceptions defined by the prelude. */
Modified: pure/trunk/printer.cc
===================================================================
--- pure/trunk/printer.cc 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/printer.cc 2008-09-16 12:07:48 UTC (rev 769)
@@ -6,6 +6,12 @@
#include <sstream>
+#include "config.h"
+
+#ifdef HAVE_GSL
+#include <gsl/gsl_matrix.h>
+#endif
+
static inline const string& pname(int32_t f)
{
assert(f > 0);
@@ -60,6 +66,7 @@
case EXPR::VAR:
case EXPR::STR:
case EXPR::PTR:
+ case EXPR::MATRIX:
return 100;
case EXPR::FVAR:
return sym_nprec(x.vtag());
@@ -160,6 +167,8 @@
return os << "::double";
case EXPR::STR:
return os << "::string";
+ case EXPR::MATRIX:
+ return os << "::matrix";
default:
return os;
}
@@ -246,6 +255,30 @@
}
case EXPR::PTR:
return os << "#<pointer " << x.pval() << ">";
+ case EXPR::MATRIX: {
+ os << "{";
+ for (exprll::const_iterator xs = x.xvals()->begin(),
+ end = x.xvals()->end(); xs != end; ) {
+ size_t n = xs->size();
+ if (n>1 || n==1 && xs->front().is_pair()) {
+ // matrix elements at a precedence not larger than ',' have to be
+ // parenthesized
+ prec_t p = sym_nprec(interpreter::g_interp->symtab.pair_sym().f) + 1;
+ for (exprl::const_iterator it = xs->begin(), end = xs->end();
+ it != end; ) {
+ os << paren(p, *it, pat);
+ if (++it != end) os << ",";
+ }
+ } else
+ for (exprl::const_iterator it = xs->begin(), end = xs->end();
+ it != end; ) {
+ printx(os, *it, pat);
+ if (++it != end) os << ",";
+ }
+ if (++xs != end) os << ";";
+ }
+ return os << "}";
+ }
case EXPR::APP: {
expr u, v, w, y;
exprl xs;
@@ -565,6 +598,9 @@
switch (x->tag) {
case EXPR::STR:
case EXPR::PTR:
+ case EXPR::MATRIX:
+ case EXPR::CMATRIX:
+ case EXPR::IMATRIX:
return 100;
case EXPR::INT:
if (x->data.i < 0)
@@ -718,6 +754,54 @@
}
case EXPR::PTR:
return os << "#<pointer " << x->data.p << ">";
+ /* NOTE: For performance reasons, we don't do any custom representations for
+ matrix elements. As a workaround, you can define __show__ on matrices as
+ a whole. */
+ case EXPR::MATRIX:
+ os << "{";
+ if (x->data.mat && x->data.mat->p) {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat->p;
+ for (size_t i = 0; i < m->size1; i++) {
+ if (i > 0) os << ";";
+ for (size_t j = 0; j < m->size2; j++) {
+ if (j > 0) os << ",";
+ os << m->data[i * m->tda + j];
+ }
+ }
+ }
+ return os << "}";
+ case EXPR::IMATRIX:
+ os << "{";
+ if (x->data.mat && x->data.mat->p) {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat->p;
+ for (size_t i = 0; i < m->size1; i++) {
+ if (i > 0) os << ";";
+ for (size_t j = 0; j < m->size2; j++) {
+ if (j > 0) os << ",";
+ os << m->data[i * m->tda + j];
+ }
+ }
+ }
+ return os << "}";
+ case EXPR::CMATRIX:
+ os << "{";
+ if (x->data.mat && x->data.mat->p) {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat->p;
+ for (size_t i = 0; i < m->size1; i++) {
+ if (i > 0) os << ";";
+ for (size_t j = 0; j < m->size2; j++) {
+ if (j > 0) os << ",";
+ /* GSL represents complex matrices using pairs of double values.
+ FIXME: We take a shortcut here by just printing complex numbers
+ in rectangular format using the +: operator defined in math.pure.
+ This has to be adapted when the representation in math.pure
+ changes. */
+ os << m->data[2*(i * m->tda + j)] << "+:"
+ << m->data[2*(i * m->tda + j) + 1];
+ }
+ }
+ }
+ return os << "}";
case EXPR::APP: {
list<const pure_expr*> xs;
prec_t p;
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/runtime.cc 2008-09-16 12:07:48 UTC (rev 769)
@@ -36,6 +36,7 @@
#ifdef HAVE_GSL
#include <gsl/gsl_errno.h>
+#include <gsl/gsl_matrix.h>
#endif
// Hooks to report stack overflows and other kinds of hard errors.
@@ -127,6 +128,12 @@
return pure_const(interpreter::g_interp->symtab.segfault_sym().f);
}
+static inline pure_expr* not_implemented_exception()
+{
+ if (!interpreter::g_interp) return 0;
+ return pure_const(interpreter::g_interp->symtab.not_implemented_sym().f);
+}
+
static inline pure_expr *get_sentry(pure_expr *x)
{
if (x==0)
@@ -263,6 +270,40 @@
return ret;
}
+static void pure_free_matrix(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ assert(x->data.mat && "pure_free_matrix: null data");
+ assert(x->data.mat->refc > 0 && "pure_free_matrix: unreferenced data");
+ assert(x->data.mat->p && "pure_free_matrix: corrupt data");
+ if (--x->data.mat->refc == 0) {
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat->p;
+ m->owner = 1;
+ gsl_matrix_free(m);
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat->p;
+ m->owner = 1;
+ gsl_matrix_complex_free(m);
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat->p;
+ m->owner = 1;
+ gsl_matrix_int_free(m);
+ break;
+ }
+ default:
+ assert(0 && "pure_free_matrix: corrupt data");
+ break;
+ }
+ }
+#endif
+}
+
#if 1
/* This is implemented (mostly) non-recursively to prevent stack overflows,
@@ -286,6 +327,7 @@
goto loop;
case EXPR::INT:
case EXPR::DBL:
+ case EXPR::PTR:
// nothing to do
break;
case EXPR::BIGINT:
@@ -294,8 +336,10 @@
case EXPR::STR:
free(x->data.s);
break;
- case EXPR::PTR:
- // noop right now, should provide installable hook in the future
+ case EXPR::MATRIX:
+ case EXPR::CMATRIX:
+ case EXPR::IMATRIX:
+ if (x->data.mat) pure_free_matrix(x);
break;
default:
assert(x->tag >= 0);
@@ -542,6 +586,110 @@
}
extern "C"
+pure_expr *pure_double_matrix(void *p)
+{
+#ifdef HAVE_GSL
+ return 0; // XXXTODO
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_complex_matrix(void *p)
+{
+#ifdef HAVE_GSL
+ return 0; // XXXTODO
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_int_matrix(void *p)
+{
+#ifdef HAVE_GSL
+ return 0; // XXXTODO
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_double_matrix_dup(const void *p)
+{
+#ifdef HAVE_GSL
+ return 0; // XXXTODO
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_complex_matrix_dup(const void *p)
+{
+#ifdef HAVE_GSL
+ return 0; // XXXTODO
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_int_matrix_dup(const void *p)
+{
+#ifdef HAVE_GSL
+ return 0; // XXXTODO
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_matrix_rowsl(uint32_t n, ...)
+{
+#ifdef HAVE_GSL
+ // XXXTODO
+ return 0;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_matrix_rowsv(uint32_t n, pure_expr **elems)
+{
+#ifdef HAVE_GSL
+ // XXXTODO
+ return 0;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_matrix_columnsl(uint32_t n, ...)
+{
+#ifdef HAVE_GSL
+ // XXXTODO
+ return 0;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_matrix_columnsv(uint32_t n, pure_expr **elems)
+{
+#ifdef HAVE_GSL
+ // XXXTODO
+ return 0;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
pure_expr *pure_app(pure_expr *fun, pure_expr *arg)
{
return pure_apply2(fun, arg);
@@ -726,6 +874,36 @@
}
extern "C"
+bool pure_is_double_matrix(const pure_expr *x, const void **p)
+{
+ if (x->tag == EXPR::MATRIX && x->data.mat) {
+ *p = x->data.mat->p;
+ return true;
+ } else
+ return false;
+}
+
+extern "C"
+bool pure_is_complex_matrix(const pure_expr *x, const void **p)
+{
+ if (x->tag == EXPR::CMATRIX && x->data.mat) {
+ *p = x->data.mat->p;
+ return true;
+ } else
+ return false;
+}
+
+extern "C"
+bool pure_is_int_matrix(const pure_expr *x, const void **p)
+{
+ if (x->tag == EXPR::IMATRIX && x->data.mat) {
+ *p = x->data.mat->p;
+ return true;
+ } else
+ return false;
+}
+
+extern "C"
bool pure_is_app(const pure_expr *x, pure_expr **fun, pure_expr **arg)
{
assert(x);
@@ -1355,7 +1533,43 @@
return &x->data.z;
}
+static inline pure_expr* bad_matrix_exception(pure_expr *x)
+{
+ if (!interpreter::g_interp) return 0;
+ pure_expr *f = pure_const(interpreter::g_interp->symtab.bad_matrix_sym().f);
+ if (x)
+ return pure_apply2(f, x);
+ else
+ return f;
+}
+
extern "C"
+pure_expr *pure_matrix_rows(uint32_t n, ...)
+{
+#ifdef HAVE_GSL
+ // XXXTODO
+ pure_throw(not_implemented_exception());
+ return 0;
+#else
+ pure_throw(not_implemented_exception());
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *pure_matrix_columns(uint32_t n, ...)
+{
+#ifdef HAVE_GSL
+ // XXXTODO
+ pure_throw(not_implemented_exception());
+ return 0;
+#else
+ pure_throw(not_implemented_exception());
+ return 0;
+#endif
+}
+
+extern "C"
pure_expr *pure_call(pure_expr *x)
{
char test;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/runtime.h 2008-09-16 12:07:48 UTC (rev 769)
@@ -30,6 +30,16 @@
bool thunked; // thunked closure? (kept unevaluated)
} pure_closure;
+/* Matrix data. The GSL matrix data is represented as a void* whose actual
+ type depends on the expression tag. Different expressions may share the
+ same underlying memory block, so we do our own reference counting to manage
+ these. */
+
+typedef struct {
+ uint32_t refc; // reference counter
+ void *p; // pointer to GSL matrix struct
+} pure_matrix;
+
/* The runtime expression data structure. Keep this lean and mean. */
typedef struct _pure_expr {
@@ -44,6 +54,7 @@
double d; // double (EXPR::DBL)
char *s; // C string (EXPR::STR)
void *p; // generic pointer (EXPR::PTR)
+ pure_matrix *mat; // matrix data (EXPR::MATRIX et al)
pure_closure *clos; // closure (0 if none)
} data;
/* Internal fields (DO NOT TOUCH). The JIT doesn't care about these. */
@@ -115,6 +126,37 @@
pure_expr *pure_string(char *s);
pure_expr *pure_cstring(char *s);
+/* Matrix constructors. The given pointer must point to a valid GSL matrix
+ struct of the corresponding GSL matrix type (gsl_matrix, gsl_matrix_complex,
+ gsl_matrix_int). (These are just given as void* here to avoid depending on
+ the GSL headers which might not be available for some implementations.) In
+ the case of the _matrix routines, the matrix must be allocated dynamically
+ and Pure takes ownership of the matrix. The matrix_dup routines first take
+ a copy of the matrix, so the ownership of the original matrix remains with
+ the caller. The result is a Pure expression representing the matrix object,
+ or null if GSL matrix support is not available or some other error
+ occurs. */
+
+pure_expr *pure_double_matrix(void *p);
+pure_expr *pure_complex_matrix(void *p);
+pure_expr *pure_int_matrix(void *p);
+pure_expr *pure_double_matrix_dup(const void *p);
+pure_expr *pure_complex_matrix_dup(const void *p);
+pure_expr *pure_int_matrix_dup(const void *p);
+
+/* Convenience functions to construct a Pure matrix from a vector or a varargs
+ list of element expressions, which can be component matrices or scalars.
+ The pure_matrix_rows functions arrange the elements vertically, while the
+ pure_matrix_columns functions arrange them horizontally, given that the
+ other dimensions match. The elems vectors are owned by the caller and won't
+ be freed. A null expression is returned in case of an error (no matrix
+ support, dimension mismatch, or invalid element type). */
+
+pure_expr *pure_matrix_rowsl(uint32_t n, ...);
+pure_expr *pure_matrix_rowsv(uint32_t n, pure_expr **elems);
+pure_expr *pure_matrix_columnsl(uint32_t n, ...);
+pure_expr *pure_matrix_columnsv(uint32_t n, pure_expr **elems);
+
/* Function applications. pure_app applies the given function to the given
argument. The result is evaluated if possible (i.e., if it is a saturated
function call). Otherwise, the result is a literal application and
@@ -172,6 +214,14 @@
bool pure_is_string_dup(const pure_expr *x, char **s);
bool pure_is_cstring_dup(const pure_expr *x, char **s);
+/* Matrix deconstructors. The returned GSL matrix pointer (represented as a
+ const void*) points to memory owned by Pure which should be considered
+ read-only and must not be freed. */
+
+bool pure_is_double_matrix(const pure_expr *x, const void **p);
+bool pure_is_complex_matrix(const pure_expr *x, const void **p);
+bool pure_is_int_matrix(const pure_expr *x, const void **p);
+
/* Deconstruct literal applications. */
bool pure_is_app(const pure_expr *x, pure_expr **fun, pure_expr **arg);
@@ -345,6 +395,14 @@
int64_t pure_get_long(pure_expr *x);
int32_t pure_get_int(pure_expr *x);
+/* Additional matrix constructors. These work like pure_matrix_rowsl and
+ pure_matrix_columnsl in the public API, but are intended to be called
+ directly from generated code and raise the appropriate Pure exceptions in
+ case of an error condition. */
+
+pure_expr *pure_matrix_rows(uint32_t n, ...);
+pure_expr *pure_matrix_columns(uint32_t n, ...);
+
/* Execute a closure. If the given expression x (or x y in the case of
pure_apply) is a parameterless closure (or a saturated application of a
closure), call it with the appropriate parameters and environment, if
Modified: pure/trunk/symtable.cc
===================================================================
--- pure/trunk/symtable.cc 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/symtable.cc 2008-09-16 12:07:48 UTC (rev 769)
@@ -34,12 +34,17 @@
fdiv_sym();
div_sym();
mod_sym();
+ // complex_rect_sym() and complex_polar_sym() are not initialized here, as
+ // they're supposed to come from math.pure which is not included in the
+ // prelude
catch_sym();
catmap_sym();
failed_match_sym();
failed_cond_sym();
signal_sym();
segfault_sym();
+ not_implemented_sym();
+ bad_matrix_sym();
amp_sym();
}
@@ -363,6 +368,24 @@
return sym("mod", 7, infixl);
}
+symbol& symtable::complex_rect_sym()
+{
+ symbol *_sym = lookup("+:");
+ if (_sym)
+ return *_sym;
+ else
+ return sym("+", 5, infix);
+}
+
+symbol& symtable::complex_polar_sym()
+{
+ symbol *_sym = lookup("<:");
+ if (_sym)
+ return *_sym;
+ else
+ return sym("+", 5, infix);
+}
+
symbol& symtable::amp_sym()
{
symbol *_sym = lookup("&");
Modified: pure/trunk/symtable.hh
===================================================================
--- pure/trunk/symtable.hh 2008-09-16 07:35:33 UTC (rev 768)
+++ pure/trunk/symtable.hh 2008-09-16 12:07:48 UTC (rev 769)
@@ -90,12 +90,16 @@
symbol& fdiv_sym();
symbol& div_sym();
symbol& mod_sym();
+ symbol& complex_rect_sym();
+ symbol& complex_polar_sym();
symbol& catch_sym() { return sym("catch"); }
symbol& catmap_sym() { return sym("catmap"); }
symbol& failed_match_sym() { return sym("failed_match"); }
symbol& failed_cond_sym() { return sym("failed_cond"); }
symbol& signal_sym() { return sym("signal"); }
symbol& segfault_sym() { return sym("stack_fault"); }
+ symbol& not_implemented_sym() { return sym("not_implemented"); }
+ symbol& bad_matrix_sym() { return sym("bad_matrix_value"); }
symbol& amp_sym();
};
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
