pure-lang-svn Mailing List for Pure (Page 4)
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
pure-lang-svn — SVN commits.
You can subscribe to this list here.
| 2008 |
Jan
|
Feb
|
Mar
|
Apr
(5) |
May
(141) |
Jun
(184) |
Jul
(97) |
Aug
(232) |
Sep
(196) |
Oct
|
Nov
|
Dec
|
|---|
|
From: <ag...@us...> - 2008-09-20 17:56:50
|
Revision: 810
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=810&view=rev
Author: agraef
Date: 2008-09-20 17:56:26 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/interpreter.cc
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-20 16:46:14 UTC (rev 809)
+++ pure/trunk/interpreter.cc 2008-09-20 17:56:26 UTC (rev 810)
@@ -2758,13 +2758,15 @@
comp_clause_list::iterator cs,
comp_clause_list::iterator end)
{
- if (cs == end)
- return expr::cons(x, expr::nil());
- else {
+ if (cs == end) {
+ exprll *xs = new exprll(1, exprl(1, x));
+ return expr(EXPR::MATRIX, xs);
+ } else {
comp_clause& c = *cs;
if (c.second.is_null()) {
expr p = c.first;
- return expr::cond(p, mkmatcomp_expr(x, n, ++cs, end), expr::nil());
+ return expr::cond(p, mkmatcomp_expr(x, n, ++cs, end),
+ expr(EXPR::MATRIX, new exprll));
} else {
expr pat = c.first, body = mkmatcomp_expr(x, n-1, ++cs, end),
arg = c.second;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 16:46:36
|
Revision: 809
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=809&view=rev
Author: agraef
Date: 2008-09-20 16:46:14 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Fix deprecated list comprehension syntax.
Modified Paths:
--------------
pure/trunk/examples/hello.pure
pure/trunk/examples/libor/myutils.pure
pure/trunk/examples/libor/queens.pure
Modified: pure/trunk/examples/hello.pure
===================================================================
--- pure/trunk/examples/hello.pure 2008-09-20 16:39:35 UTC (rev 808)
+++ pure/trunk/examples/hello.pure 2008-09-20 16:46:14 UTC (rev 809)
@@ -141,7 +141,7 @@
the same little 2.5-liner. ;-) */
qsort p [] = [];
-qsort p (x:xs) = qsort p [l; l = xs; l<x] + (x : qsort p [r; r = xs; r>=x])
+qsort p (x:xs) = qsort p [l | l = xs; l<x] + (x : qsort p [r | r = xs; r>=x])
with x<y = p x y; x>=y = not p x y end;
qsort (<) (1..20); // ascending sort, no-op in this case
@@ -153,15 +153,15 @@
extern int rand();
-qsort (<) [rand; i = 1..20]; // sort 20 random numbers in ascending order
-qsort (>) [rand; i = 1..20]; // sort 20 random numbers in descending order
+qsort (<) [rand|i = 1..20]; // sort 20 random numbers in ascending order
+qsort (>) [rand|i = 1..20]; // sort 20 random numbers in descending order
/* Erathosthenes' classical prime sieve. Another example showing the beauty
and usefulness of list comprehensions. */
primes n = sieve (2..n) with
sieve [] = [];
- sieve (p:qs) = p : sieve [q; q = qs; q mod p];
+ sieve (p:qs) = p : sieve [q | q = qs; q mod p];
end;
primes 100;
@@ -170,7 +170,7 @@
prime numbers. */
all_primes = sieve (2..inf) with
- sieve (p:qs) = p : sieve [q; q = qs; q mod p] &;
+ sieve (p:qs) = p : sieve [q | q = qs; q mod p] &;
end;
// Assign this to a variable, so we can take advantage of memoization.
@@ -193,7 +193,7 @@
queens n = search n 1 []
with
search n i p = [reverse p] if i>n;
- = cat [search n (i+1) ((i,j):p); j = 1..n; safe (i,j) p];
+ = cat [search n (i+1) ((i,j):p) | j = 1..n; safe (i,j) p];
safe (i,j) p = not any (check (i,j)) p;
check (i1,j1) (i2,j2)
= i1==i2 || j1==j2 || i1+j1==i2+j2 || i1-j1==i2-j2;
@@ -206,7 +206,7 @@
queens1 n = catch reverse (search n 1 [])
with
search n i p = throw p if i>n;
- = void [search n (i+1) ((i,j):p); j = 1..n; safe (i,j) p];
+ = void [search n (i+1) ((i,j):p) | j = 1..n; safe (i,j) p];
safe (i,j) p = not any (check (i,j)) p;
check (i1,j1) (i2,j2)
= i1==i2 || j1==j2 || i1+j1==i2+j2 || i1-j1==i2-j2;
@@ -286,4 +286,4 @@
// A random example:
-let xs = [rand; i = 1..20]; let T = bintree xs; xs; T; members T;
+let xs = [rand|i = 1..20]; let T = bintree xs; xs; T; members T;
Modified: pure/trunk/examples/libor/myutils.pure
===================================================================
--- pure/trunk/examples/libor/myutils.pure 2008-09-20 16:39:35 UTC (rev 808)
+++ pure/trunk/examples/libor/myutils.pure 2008-09-20 16:46:14 UTC (rev 809)
@@ -4,7 +4,7 @@
MathIter1 op i1 i2 f = foldl1 op (map f (i1..i2));
MathIter2 op i1 i2 j1 j2 f =
- foldl1 op (map (uncurry f) [x,y; x = i1..i2; y = j1..j2]);
+ foldl1 op (map (uncurry f) [x,y | x = i1..i2; y = j1..j2]);
//Examples on how to use the mathematical iterators
Sigma i1 i2 f = MathIter1 (+) i1 i2 f;
Modified: pure/trunk/examples/libor/queens.pure
===================================================================
--- pure/trunk/examples/libor/queens.pure 2008-09-20 16:39:35 UTC (rev 808)
+++ pure/trunk/examples/libor/queens.pure 2008-09-20 16:46:14 UTC (rev 809)
@@ -20,7 +20,7 @@
with
searchall n::int 0 p = p;
searchall n::int i::int p =
- tuple [searchall n (i-1) (j:p); j = 1..n; safe 1 j p]
+ tuple [searchall n (i-1) (j:p) | j = 1..n; safe 1 j p]
end;
// the solution is only the rows permutation, without the ordered columns (1..n)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 16:39:45
|
Revision: 808
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=808&view=rev
Author: agraef
Date: 2008-09-20 16:39:35 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Syntax changes (list and array comprehension syntax).
Modified Paths:
--------------
pure/trunk/interpreter.cc
pure/trunk/interpreter.hh
pure/trunk/parser.yy
pure/trunk/symtable.cc
pure/trunk/symtable.hh
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-20 16:37:50 UTC (rev 807)
+++ pure/trunk/interpreter.cc 2008-09-20 16:39:35 UTC (rev 808)
@@ -2754,6 +2754,40 @@
return new expr(y);
}
+expr interpreter::mkmatcomp_expr(expr x, size_t n,
+ comp_clause_list::iterator cs,
+ comp_clause_list::iterator end)
+{
+ if (cs == end)
+ return expr::cons(x, expr::nil());
+ else {
+ comp_clause& c = *cs;
+ if (c.second.is_null()) {
+ expr p = c.first;
+ return expr::cond(p, mkmatcomp_expr(x, n, ++cs, end), expr::nil());
+ } else {
+ expr pat = c.first, body = mkmatcomp_expr(x, n-1, ++cs, end),
+ arg = c.second;
+ closure(pat, body);
+ expr f = (n&1)?symtab.colcatmap_sym().x:symtab.rowcatmap_sym().x;
+ return expr(f, expr::lambda(pat, body), arg);
+ }
+ }
+}
+
+expr *interpreter::mkmatcomp_expr(expr *x, comp_clause_list *cs)
+{
+ size_t n = 0;
+ for (comp_clause_list::iterator it = cs->begin(), end = cs->end();
+ it != end; it++) {
+ comp_clause& c = *it;
+ if (!c.second.is_null()) n++;
+ }
+ expr y = mkmatcomp_expr(*x, n, cs->begin(), cs->end());
+ delete x; delete cs;
+ return new expr(y);
+}
+
// Code generation.
#define Dbl(d) ConstantFP::get(Type::DoubleTy, d)
Modified: pure/trunk/interpreter.hh
===================================================================
--- pure/trunk/interpreter.hh 2008-09-20 16:37:50 UTC (rev 807)
+++ pure/trunk/interpreter.hh 2008-09-20 16:39:35 UTC (rev 808)
@@ -483,6 +483,9 @@
expr *mklistcomp_expr(expr *x, comp_clause_list *cs);
expr mklistcomp_expr(expr x, comp_clause_list::iterator cs,
comp_clause_list::iterator end);
+ expr *mkmatcomp_expr(expr *x, comp_clause_list *cs);
+ expr mkmatcomp_expr(expr x, size_t n, comp_clause_list::iterator cs,
+ comp_clause_list::iterator end);
// LLVM code generation and execution.
Modified: pure/trunk/parser.yy
===================================================================
--- pure/trunk/parser.yy 2008-09-20 16:37:50 UTC (rev 807)
+++ pure/trunk/parser.yy 2008-09-20 16:39:35 UTC (rev 808)
@@ -526,8 +526,14 @@
| DBL { $$ = new expr(EXPR::DBL, $1); }
| STR { $$ = new expr(EXPR::STR, $1); }
| '{' rows '}' { $$ = new expr(EXPR::MATRIX, $2); }
+| '{' expr '|' comp_clauses '}'
+ { $$ = interp.mkmatcomp_expr($2, $4); }
| '[' expr ']' { $$ = interp.mklist_expr($2); }
| '[' expr ';' comp_clauses ']'
+ { interp.warning(yyloc,
+ "warning: deprecated comprehension syntax");
+ $$ = interp.mklistcomp_expr($2, $4); }
+| '[' expr '|' comp_clauses ']'
{ $$ = interp.mklistcomp_expr($2, $4); }
| '(' expr ')' { $$ = $2;
if ($$->is_pair()) $$->flags() |= EXPR::PAREN; }
Modified: pure/trunk/symtable.cc
===================================================================
--- pure/trunk/symtable.cc 2008-09-20 16:37:50 UTC (rev 807)
+++ pure/trunk/symtable.cc 2008-09-20 16:39:35 UTC (rev 808)
@@ -34,11 +34,10 @@
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();
+ rowcatmap_sym();
+ colcatmap_sym();
failed_match_sym();
failed_cond_sym();
signal_sym();
Modified: pure/trunk/symtable.hh
===================================================================
--- pure/trunk/symtable.hh 2008-09-20 16:37:50 UTC (rev 807)
+++ pure/trunk/symtable.hh 2008-09-20 16:39:35 UTC (rev 808)
@@ -92,6 +92,8 @@
symbol& mod_sym();
symbol& catch_sym() { return sym("catch"); }
symbol& catmap_sym() { return sym("catmap"); }
+ symbol& rowcatmap_sym() { return sym("rowcatmap"); }
+ symbol& colcatmap_sym() { return sym("colcatmap"); }
symbol& failed_match_sym() { return sym("failed_match"); }
symbol& failed_cond_sym() { return sym("failed_cond"); }
symbol& signal_sym() { return sym("signal"); }
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 16:37:56
|
Revision: 807
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=807&view=rev
Author: agraef
Date: 2008-09-20 16:37:50 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Fix deprecated list comprehension syntax.
Modified Paths:
--------------
pure/trunk/lib/strings.pure
Modified: pure/trunk/lib/strings.pure
===================================================================
--- pure/trunk/lib/strings.pure 2008-09-20 16:36:47 UTC (rev 806)
+++ pure/trunk/lib/strings.pure 2008-09-20 16:37:50 UTC (rev 807)
@@ -137,7 +137,7 @@
join delim::string [] = "";
join delim::string (x::string:xs)
- = x + strcat [delim+x; x = xs]
+ = x + strcat [delim+x | x = xs]
if listp xs && all stringp xs;
/* Split a string into parts delimited by the given (nonempty) string. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 16:37:02
|
Revision: 806
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=806&view=rev
Author: agraef
Date: 2008-09-20 16:36:47 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Fix deprecated list comprehension syntax.
Modified Paths:
--------------
pure/trunk/test/test014.pure
pure/trunk/test/test015.pure
pure/trunk/test/test019.pure
pure/trunk/test/test020.pure
pure/trunk/test/test021.pure
Modified: pure/trunk/test/test014.pure
===================================================================
--- pure/trunk/test/test014.pure 2008-09-20 09:49:59 UTC (rev 805)
+++ pure/trunk/test/test014.pure 2008-09-20 16:36:47 UTC (rev 806)
@@ -3,7 +3,7 @@
// Random number generator.
-drop 97 [random; i=1..100] when () = srandom 0 end;
+drop 97 [random | i=1..100] when () = srandom 0 end;
// Some rational arithmetic tests, pilfered from Rob Hubbard's Q+Q manual.
Modified: pure/trunk/test/test015.pure
===================================================================
--- pure/trunk/test/test015.pure 2008-09-20 09:49:59 UTC (rev 805)
+++ pure/trunk/test/test015.pure 2008-09-20 16:36:47 UTC (rev 806)
@@ -123,18 +123,18 @@
********************************************/
// Create data structures
-let a = dict [(i => (double i)); i = 1..10];
-let b = dict [(i => (double i)); i = 11..20];
+let a = dict [i => (double i) | i = 1..10];
+let b = dict [i => (double i) | i = 11..20];
-let c = hdict $ zipwith (=>) [(i, double i, str i); i = 1..10] ( 1..10);
-let d = hdict $ zipwith (=>) [(i, double i, str i); i = 11..20] (11..20);
+let c = hdict $ zipwith (=>) [i, double i, str i | i = 1..10] ( 1..10);
+let d = hdict $ zipwith (=>) [i, double i, str i | i = 11..20] (11..20);
let e = dict $ zipwith (=>) (map str (1..10)) (map str (1..10));
a; b; c; d; e;
mkdict 1000 (1..10);
-mkhdict 1000 [(i, double i, str i); i = 1..10];
+mkhdict 1000 [i, double i, str i | i = 1..10];
// Type tests
@@ -167,7 +167,7 @@
// Slicing
a!!(5..15);
-c!![(i, double i, str i); i = 5..15];
+c!![i, double i, str i | i = 5..15];
// Convert data set to list
@@ -185,7 +185,7 @@
foldl delete a (1..10);
delete a 5000;
-foldl delete c [(i, double i, str i); i = 1..10];
+foldl delete c [i, double i, str i | i = 1..10];
delete c (5000, 5000.0, "5000");
// Relations
@@ -259,8 +259,8 @@
let a = array ( 1..10);
let b = array (11..20);
-let c = array2 [[str (i + j); i = j..10]; j = 1..10];
-let d = array2 [[str (i + j); i = j..10]; j = 11..20];
+let c = array2 [[str (i + j) | i = j..10] | j = 1..10];
+let d = array2 [[str (i + j) | i = j..10] | j = 11..20];
a; b; c; d;
@@ -293,7 +293,7 @@
// Slicing
a!!(5..15);
-c!![(i, j); i = 3..15; j = 3..15];
+c!![(i, j) | i = 3..15; j = 3..15];
// Convert data set to list
Modified: pure/trunk/test/test019.pure
===================================================================
--- pure/trunk/test/test019.pure 2008-09-20 09:49:59 UTC (rev 805)
+++ pure/trunk/test/test019.pure 2008-09-20 16:36:47 UTC (rev 806)
@@ -1,3 +1,3 @@
// const name capture regression, reported by Eddie Rucker
const i = -1;
-[(i,j); i=1..2; j=3..4];
+[(i,j) | i=1..2; j=3..4];
Modified: pure/trunk/test/test020.pure
===================================================================
--- pure/trunk/test/test020.pure 2008-09-20 09:49:59 UTC (rev 805)
+++ pure/trunk/test/test020.pure 2008-09-20 16:36:47 UTC (rev 806)
@@ -21,7 +21,7 @@
// binary operations
let f2 = [(+), (-), (*), (/), (^), atan2, pow];
-let x2 = [i,j; i=x; j=x];
+let x2 = [i,j | i=x; j=x];
// test instrumentation (AG)
@@ -58,7 +58,7 @@
print x = str x otherwise;
tests =
-puts "*** UNARY ***" $$ void [test (f,x); f=f; x=x] $$
-puts "*** BINARY ***" $$ void [test (f,x); f=f2; x=x2];
+puts "*** UNARY ***" $$ void [test (f,x) | f=f; x=x] $$
+puts "*** BINARY ***" $$ void [test (f,x) | f=f2; x=x2];
tests;
Modified: pure/trunk/test/test021.pure
===================================================================
--- pure/trunk/test/test021.pure 2008-09-20 09:49:59 UTC (rev 805)
+++ pure/trunk/test/test021.pure 2008-09-20 16:36:47 UTC (rev 806)
@@ -36,11 +36,11 @@
puts "*** EXACT/INEXACT ***" $$
// These should all return exact results, except +/- with polar operands, as
// well as / and ^ which always return inexact results.
-do test [op,2,a;op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]] $$
-do test [op,a,2;op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]] $$
+do test [op,2,a|op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]] $$
+do test [op,a,2|op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]] $$
puts "*** SYMBOLIC ***" $$
// If everything is all right here, these should all print __failed__.
-do test [op,x,a;op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]] $$
-do test [op,a,x;op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]];
+do test [op,x,a|op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]] $$
+do test [op,a,x|op=[(+),(-),(*),(%),(/),(^)];a=[2+:3,2<:3,2%3]];
tests;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 09:50:09
|
Revision: 805
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=805&view=rev
Author: agraef
Date: 2008-09-20 09:49:59 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Cosmetic changes to arithmetic sequence syntax, to make it easier to write matrix slices.
Modified Paths:
--------------
pure/trunk/lib/prelude.pure
pure/trunk/test/prelude.log
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-20 09:24:14 UTC (rev 804)
+++ pure/trunk/lib/prelude.pure 2008-09-20 09:49:59 UTC (rev 805)
@@ -45,8 +45,8 @@
infixl 0 $$ ; // sequence operator
infixr 0 $ ; // right-associative application
-infix 1 .. ; // arithmetic sequences
infixr 1 , ; // pair (tuple)
+infix 2 .. ; // arithmetic sequences
infix 2 => ; // mapsto constructor
infixr 2 || ; // logical or (short-circuit)
infixr 3 && ; // logical and (short-circuit)
@@ -254,7 +254,7 @@
/* Arithmetic sequences. */
-n1,n2..m = if m===s*inf then iterate (\x->x+k) n1
+n1:n2..m = if m===s*inf then iterate (\x->x+k) n1
else while (\i->s*i<=s*m) (\x->x+k) n1
when k = n2-n1; s = if k>0 then 1 else -1 end if n1!=n2;
n..m = if m===inf then iterate (\x->x+1) n
Modified: pure/trunk/test/prelude.log
===================================================================
--- pure/trunk/test/prelude.log 2008-09-20 09:24:14 UTC (rev 804)
+++ pure/trunk/test/prelude.log 2008-09-20 09:49:59 UTC (rev 805)
@@ -319,7 +319,7 @@
<var> state 2
state 2: #0
} end;
-n1/*0:0101*/,n2/*0:011*/..m/*0:1*/ = if m/*2:1*/===s/*0:*/*inf then iterate (\x/*0:*/ -> x/*0:*/+k/*2:*/ {
+n1/*0:0101*/:n2/*0:011*/..m/*0:1*/ = if m/*2:1*/===s/*0:*/*inf then iterate (\x/*0:*/ -> x/*0:*/+k/*2:*/ {
rule #0: x = x+k
state 0: #0
<var> state 1
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 09:24:41
|
Revision: 804
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=804&view=rev
Author: agraef
Date: 2008-09-20 09:24:14 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Comment change.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-20 09:21:34 UTC (rev 803)
+++ pure/trunk/lib/primitives.pure 2008-09-20 09:24:14 UTC (rev 804)
@@ -470,8 +470,7 @@
extern expr* matrix_rows(expr *x) = rowcat;
extern expr* matrix_columns(expr *x) = colcat;
-/* Combinations of rowcat/colcat and map, to be used in matrix
- comprehensions. */
+/* Combinations of rowcat/colcat and map. */
rowcatmap f [] = {};
rowcatmap f xs@(_:_) = rowcat (map f xs);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 09:21:41
|
Revision: 803
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=803&view=rev
Author: agraef
Date: 2008-09-20 09:21:34 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-20 08:57:43 UTC (rev 802)
+++ pure/trunk/runtime.cc 2008-09-20 09:21:34 UTC (rev 803)
@@ -973,6 +973,10 @@
data[i*tda+k] = (double)mat1->data[j*mat1->tda+k];
break;
}
+ case EXPR::CMATRIX:
+ case EXPR::MATRIX:
+ // empty matrix, skip
+ break;
default:
assert(0 && "bad matrix element");
break;
@@ -1021,6 +1025,10 @@
i += mat1->size2;
break;
}
+ case EXPR::CMATRIX:
+ case EXPR::MATRIX:
+ // empty matrix, skip
+ break;
default:
assert(0 && "bad matrix element");
break;
@@ -1096,6 +1104,9 @@
ncols*2*sizeof(double));
break;
}
+ case EXPR::MATRIX:
+ // empty matrix, skip
+ break;
default:
assert(0 && "bad matrix element");
break;
@@ -1174,6 +1185,9 @@
i += mat1->size2;
break;
}
+ case EXPR::MATRIX:
+ // empty matrix, skip
+ break;
default:
assert(0 && "bad matrix element");
break;
@@ -1219,6 +1233,10 @@
memcpy(data+i*tda, mat1->data+j*mat1->tda, ncols*sizeof(int));
break;
}
+ case EXPR::CMATRIX:
+ case EXPR::MATRIX:
+ // empty matrix, skip
+ break;
default:
assert(0 && "bad matrix element");
break;
@@ -1267,6 +1285,10 @@
i += mat1->size2;
break;
}
+ case EXPR::CMATRIX:
+ case EXPR::MATRIX:
+ // empty matrix, skip
+ break;
default:
assert(0 && "bad matrix element");
break;
@@ -1431,13 +1453,10 @@
switch (x->tag) {
case EXPR::MATRIX: {
gsl_matrix_symbolic *mp = (gsl_matrix_symbolic*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
return 0;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::MATRIX);
have_matrix = true;
}
@@ -1466,13 +1485,10 @@
}
case EXPR::DMATRIX: {
gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
return 0;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::DMATRIX);
have_matrix = true;
}
@@ -1480,13 +1496,10 @@
}
case EXPR::CMATRIX: {
gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
return 0;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::CMATRIX);
have_matrix = true;
}
@@ -1494,13 +1507,10 @@
}
case EXPR::IMATRIX: {
gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
return 0;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::IMATRIX);
have_matrix = true;
}
@@ -1556,13 +1566,10 @@
switch (x->tag) {
case EXPR::MATRIX: {
gsl_matrix_symbolic *mp = (gsl_matrix_symbolic*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
return 0;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::MATRIX);
have_matrix = true;
}
@@ -1591,13 +1598,10 @@
}
case EXPR::DMATRIX: {
gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
return 0;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::DMATRIX);
have_matrix = true;
}
@@ -1605,13 +1609,10 @@
}
case EXPR::CMATRIX: {
gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
return 0;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::CMATRIX);
have_matrix = true;
}
@@ -1619,13 +1620,10 @@
}
case EXPR::IMATRIX: {
gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
return 0;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- return 0;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::IMATRIX);
have_matrix = true;
}
@@ -2530,13 +2528,10 @@
switch (x->tag) {
case EXPR::MATRIX: {
gsl_matrix_symbolic *mp = (gsl_matrix_symbolic*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
goto err;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::MATRIX);
have_matrix = true;
}
@@ -2565,13 +2560,10 @@
}
case EXPR::DMATRIX: {
gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
goto err;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::DMATRIX);
have_matrix = true;
}
@@ -2579,13 +2571,10 @@
}
case EXPR::CMATRIX: {
gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
goto err;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::CMATRIX);
have_matrix = true;
}
@@ -2593,13 +2582,10 @@
}
case EXPR::IMATRIX: {
gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size2 != (size_t)k)
goto err;
nrows += mp->size1; k = mp->size2;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::IMATRIX);
have_matrix = true;
}
@@ -2662,13 +2648,10 @@
switch (x->tag) {
case EXPR::MATRIX: {
gsl_matrix_symbolic *mp = (gsl_matrix_symbolic*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
goto err;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::MATRIX);
have_matrix = true;
}
@@ -2697,13 +2680,10 @@
}
case EXPR::DMATRIX: {
gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
goto err;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::DMATRIX);
have_matrix = true;
}
@@ -2711,13 +2691,10 @@
}
case EXPR::CMATRIX: {
gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
goto err;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::CMATRIX);
have_matrix = true;
}
@@ -2725,13 +2702,10 @@
}
case EXPR::IMATRIX: {
gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
- if (mp) {
+ if (mp->size1 > 0 && mp->size2 > 0) {
if (k >= 0 && mp->size1 != (size_t)k)
goto err;
ncols += mp->size2; k = mp->size1;
- } else if (k>0)
- goto err;
- if (mp->size1 > 0 && mp->size2 > 0) {
set_target_type(target, EXPR::IMATRIX);
have_matrix = true;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 08:58:06
|
Revision: 802
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=802&view=rev
Author: agraef
Date: 2008-09-20 08:57:43 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/printer.cc
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-20 08:20:26 UTC (rev 801)
+++ pure/trunk/lib/primitives.pure 2008-09-20 08:57:43 UTC (rev 802)
@@ -460,7 +460,7 @@
/* Extract a submatrix of a given size at a given offset. */
submat x::matrix (i::int,j::int) (n::int,m::int)
- = matrix_slice x i j (i+n) (j+m);
+ = matrix_slice x i j (i+n-1) (j+m-1);
/* Construct matrices from lists of rows and columns. These take either
scalars or submatrices as inputs; corresponding dimensions must match.
Modified: pure/trunk/printer.cc
===================================================================
--- pure/trunk/printer.cc 2008-09-20 08:20:26 UTC (rev 801)
+++ pure/trunk/printer.cc 2008-09-20 08:57:43 UTC (rev 802)
@@ -764,33 +764,37 @@
}
case EXPR::PTR:
return os << "#<pointer " << x->data.p << ">";
-#ifdef HAVE_GSL
/* 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.p) {
- prec_t p = sym_nprec(interpreter::g_interp->symtab.pair_sym().f) + 1;
gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)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 << pure_paren(p, m->data[i * m->tda + j]);
+ if (m->size1>0 && m->size2>0) {
+ prec_t p = sym_nprec(interpreter::g_interp->symtab.pair_sym().f) + 1;
+ 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 << pure_paren(p, m->data[i * m->tda + j]);
+ }
}
}
}
return os << "}";
+#ifdef HAVE_GSL
case EXPR::DMATRIX:
os << "{";
if (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 << ",";
- print_double(os, m->data[i * m->tda + j]);
+ if (m->size1>0 && m->size2>0) {
+ 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 << ",";
+ print_double(os, m->data[i * m->tda + j]);
+ }
}
}
}
@@ -799,11 +803,15 @@
os << "{";
if (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];
+ if (m->size1>0 && m->size2>0) {
+ if (m->size1>0 && m->size2>0) {
+ 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];
+ }
+ }
}
}
}
@@ -812,39 +820,39 @@
os << "{";
if (x->data.mat.p) {
gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
- /* GSL represents complex matrices using pairs of double values, while
- Pure provides its own complex type in math.pure. If math.pure has
- been loaded, then the '+:' operator is defined and we use this
- representation. Otherwise, we print complex values as pairs of real
- and imaginary part. */
- symbol *rect = interpreter::g_interp->symtab.complex_rect_sym();
- if (rect)
- 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 << ",";
- print_double(os, m->data[2*(i * m->tda + j)]);
- os << rect->s;
- print_double(os, m->data[2*(i * m->tda + j) + 1]);
+ if (m->size1>0 && m->size2>0) {
+ /* GSL represents complex matrices using pairs of double values, while
+ Pure provides its own complex type in math.pure. If math.pure has
+ been loaded, then the '+:' operator is defined and we use this
+ representation. Otherwise, we print complex values as pairs of real
+ and imaginary part. */
+ symbol *rect = interpreter::g_interp->symtab.complex_rect_sym();
+ if (rect)
+ 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 << ",";
+ print_double(os, m->data[2*(i * m->tda + j)]);
+ os << rect->s;
+ print_double(os, m->data[2*(i * m->tda + j) + 1]);
+ }
}
- }
- else
- 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 << "(";
- print_double(os, m->data[2*(i * m->tda + j)]);
- os << ",";
- print_double(os, m->data[2*(i * m->tda + j) + 1]);
- os << ")";
+ else
+ 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 << "(";
+ print_double(os, m->data[2*(i * m->tda + j)]);
+ os << ",";
+ print_double(os, m->data[2*(i * m->tda + j) + 1]);
+ os << ")";
+ }
}
- }
+ }
}
return os << "}";
#else
- case EXPR::MATRIX:
- return os << "#<matrix " << x->data.mat.p << ">";
case EXPR::DMATRIX:
return os << "#<dmatrix " << x->data.mat.p << ">";
case EXPR::IMATRIX:
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-20 08:20:26 UTC (rev 801)
+++ pure/trunk/runtime.cc 2008-09-20 08:57:43 UTC (rev 802)
@@ -4053,7 +4053,7 @@
}
extern "C"
-pure_expr *matrix_elem_at(pure_expr *x, uint32_t i)
+pure_expr *matrix_elem_at(pure_expr *x, int32_t i)
{
switch (x->tag) {
case EXPR::MATRIX: {
@@ -4080,7 +4080,7 @@
}
extern "C"
-pure_expr *matrix_elem_at2(pure_expr *x, uint32_t i, uint32_t j)
+pure_expr *matrix_elem_at2(pure_expr *x, int32_t i, int32_t j)
{
switch (x->tag) {
case EXPR::MATRIX: {
@@ -4111,14 +4111,18 @@
}
extern "C"
-pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
- uint32_t i2, uint32_t j2)
+pure_expr *matrix_slice(pure_expr *x, int32_t i1, int32_t j1,
+ int32_t i2, int32_t j2)
{
void *p = 0;
+ if (i1<0) i1 = 0; if (j1<0) j1 = 0;
switch (x->tag) {
case EXPR::MATRIX: {
gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
- size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (i2 >= (int)m->size1) i2 = m->size1-1;
+ if (j2 >= (int)m->size2) j2 = m->size2-1;
+ size_t n1 = (i1<(int)m->size1 && i2>=i1)?(i2+1-i1):0,
+ n2 = (j1<(int)m->size2 && j2>=j1)?(j2+1-j1):0;
if (n1 == 0 || n2 == 0) // empty matrix
return pure_symbolic_matrix(create_symbolic_matrix(n1, n2));
gsl_matrix_symbolic_view v =
@@ -4134,7 +4138,10 @@
#ifdef HAVE_GSL
case EXPR::DMATRIX: {
gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
- size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (i2 >= (int)m->size1) i2 = m->size1-1;
+ if (j2 >= (int)m->size2) j2 = m->size2-1;
+ size_t n1 = (i1<(int)m->size1 && i2>=i1)?(i2+1-i1):0,
+ n2 = (j1<(int)m->size2 && j2>=j1)?(j2+1-j1):0;
if (n1 == 0 || n2 == 0) // empty matrix
return pure_double_matrix(create_double_matrix(n1, n2));
gsl_matrix_view v = gsl_matrix_submatrix(m, i1, j1, n1, n2);
@@ -4147,7 +4154,10 @@
}
case EXPR::CMATRIX: {
gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
- size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (i2 >= (int)m->size1) i2 = m->size1-1;
+ if (j2 >= (int)m->size2) j2 = m->size2-1;
+ size_t n1 = (i1<(int)m->size1 && i2>=i1)?(i2+1-i1):0,
+ n2 = (j1<(int)m->size2 && j2>=j1)?(j2+1-j1):0;
if (n1 == 0 || n2 == 0) // empty matrix
return pure_complex_matrix(create_complex_matrix(n1, n2));
gsl_matrix_complex_view v =
@@ -4162,7 +4172,10 @@
}
case EXPR::IMATRIX: {
gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
- size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (i2 >= (int)m->size1) i2 = m->size1-1;
+ if (j2 >= (int)m->size2) j2 = m->size2-1;
+ size_t n1 = (i1<(int)m->size1 && i2>=i1)?(i2+1-i1):0,
+ n2 = (j1<(int)m->size2 && j2>=j1)?(j2+1-j1):0;
if (n1 == 0 || n2 == 0) // empty matrix
return pure_int_matrix(create_int_matrix(n1, n2));
gsl_matrix_int_view v = gsl_matrix_int_submatrix(m, i1, j1, n1, n2);
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-20 08:20:26 UTC (rev 801)
+++ pure/trunk/runtime.h 2008-09-20 08:57:43 UTC (rev 802)
@@ -659,17 +659,17 @@
aren't range-checked, if this is needed you have to do it beforehand,
checking against matrix_size or matrix_dim above. */
-pure_expr *matrix_elem_at(pure_expr *x, uint32_t i);
-pure_expr *matrix_elem_at2(pure_expr *x, uint32_t i, uint32_t j);
+pure_expr *matrix_elem_at(pure_expr *x, int32_t i);
+pure_expr *matrix_elem_at2(pure_expr *x, int32_t i, int32_t j);
/* The following operation retrieves a slice a.k.a. submatrix of a matrix and
returns it as a new matrix object. The result matrix shares the underlying
storage with the input matrix (i.e., matrix elements are *not* copied) and
- so this is a comparatively cheap operation. Indices are zero-based and must
- be checked by the caller if necessary. */
+ so this is a comparatively cheap operation. Indices are zero-based and are
+ clamped to the available index range automatically. */
-pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
- uint32_t i2, uint32_t j2);
+pure_expr *matrix_slice(pure_expr *x, int32_t i1, int32_t j1,
+ int32_t i2, int32_t j2);
/* Matrix construction. These work like the pure_matrix_rows/
pure_matrix_columns functions in the public API, but take their input from
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 08:20:55
|
Revision: 801
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=801&view=rev
Author: agraef
Date: 2008-09-20 08:20:26 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Add submatrix operation.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-20 07:56:32 UTC (rev 800)
+++ pure/trunk/lib/primitives.pure 2008-09-20 08:20:26 UTC (rev 801)
@@ -457,6 +457,11 @@
cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
+/* Extract a submatrix of a given size at a given offset. */
+
+submat x::matrix (i::int,j::int) (n::int,m::int)
+ = matrix_slice x i j (i+n) (j+m);
+
/* Construct matrices from lists of rows and columns. These take either
scalars or submatrices as inputs; corresponding dimensions must match.
rowcat combines submatrices vertically, like {x;y}; colcat combines them
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 07:56:36
|
Revision: 800
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=800&view=rev
Author: agraef
Date: 2008-09-20 07:56:32 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Implement array slicing operations.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-20 07:03:48 UTC (rev 799)
+++ pure/trunk/lib/primitives.pure 2008-09-20 07:56:32 UTC (rev 800)
@@ -430,6 +430,16 @@
when n::int,m::int = dim x end);
= throw out_of_bounds otherwise;
+/* Slices. */
+
+x::matrix!!(ns,ms) = colcatmap (mth (rowcatmap (nth x) ns)) ms with
+ nth x n = catch (cst {}) (row x n);
+ mth x m = catch (cst {}) (col x m);
+ end;
+x::matrix!!ns = colcatmap (nth x) ns with
+ nth x n = catch (cst {}) {x!n};
+ end;
+
/* Extract rows and columns from a matrix. */
private matrix_slice;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 07:03:55
|
Revision: 799
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=799&view=rev
Author: agraef
Date: 2008-09-20 07:03:48 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-20 06:59:05 UTC (rev 798)
+++ pure/trunk/runtime.cc 2008-09-20 07:03:48 UTC (rev 799)
@@ -4183,7 +4183,7 @@
y->tag = x->tag;
y->data.mat.p = p;
y->data.mat.refc = x->data.mat.refc;
- *y->data.mat.refc++;
+ (*y->data.mat.refc)++;
MEMDEBUG_NEW(y)
return y;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 06:59:11
|
Revision: 798
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=798&view=rev
Author: agraef
Date: 2008-09-20 06:59:05 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Add some more basic matrix construction and deconstruction operations.
Modified Paths:
--------------
pure/trunk/lib/prelude.pure
pure/trunk/lib/primitives.pure
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-20 05:42:23 UTC (rev 797)
+++ pure/trunk/lib/prelude.pure 2008-09-20 06:59:05 UTC (rev 798)
@@ -45,6 +45,7 @@
infixl 0 $$ ; // sequence operator
infixr 0 $ ; // right-associative application
+infix 1 .. ; // arithmetic sequences
infixr 1 , ; // pair (tuple)
infix 2 => ; // mapsto constructor
infixr 2 || ; // logical or (short-circuit)
@@ -253,8 +254,6 @@
/* Arithmetic sequences. */
-infix 1 .. ;
-
n1,n2..m = if m===s*inf then iterate (\x->x+k) n1
else while (\i->s*i<=s*m) (\x->x+k) n1
when k = n2-n1; s = if k>0 then 1 else -1 end if n1!=n2;
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-20 05:42:23 UTC (rev 797)
+++ pure/trunk/lib/primitives.pure 2008-09-20 06:59:05 UTC (rev 798)
@@ -409,7 +409,7 @@
x::pointer==y::pointer = bigint x == bigint y;
x::pointer!=y::pointer = bigint x != bigint y;
-/* Basic matrix operations. */
+/* Basic matrix operations: size, dimensions and indexing. */
private matrix_size matrix_dim;
extern int matrix_size(expr *x), expr* matrix_dim(expr *x);
@@ -430,14 +430,49 @@
when n::int,m::int = dim x end);
= throw out_of_bounds otherwise;
-extern expr* matrix_rows(expr *x) = rowmatrix;
-extern expr* matrix_columns(expr *x) = colmatrix;
+/* Extract rows and columns from a matrix. */
+private matrix_slice;
+extern expr* matrix_slice(expr* x, int i1, int j1, int i2, int j2);
+
+row x::matrix i::int = if i>=0 && i<n then matrix_slice x i 0 i (m-1)
+ else throw out_of_bounds
+ when n::int,m::int = dim x end;
+
+col x::matrix j::int = if j>=0 && j<m then matrix_slice x 0 j (n-1) j
+ else throw out_of_bounds
+ when n::int,m::int = dim x end;
+
+rows x::matrix = map (row x) (0..n-1) when n::int,_ = dim x end;
+
+cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
+
+/* Construct matrices from lists of rows and columns. These take either
+ scalars or submatrices as inputs; corresponding dimensions must match.
+ rowcat combines submatrices vertically, like {x;y}; colcat combines them
+ horizontally, like {x,y}. */
+
+extern expr* matrix_rows(expr *x) = rowcat;
+extern expr* matrix_columns(expr *x) = colcat;
+
+/* Combinations of rowcat/colcat and map, to be used in matrix
+ comprehensions. */
+
+rowcatmap f [] = {};
+rowcatmap f xs@(_:_) = rowcat (map f xs);
+
+colcatmap f [] = {};
+colcatmap f xs@(_:_) = colcat (map f xs);
+
+/* Transpose a matrix. */
+
private matrix_transpose;
extern expr* matrix_transpose(expr *x);
x::matrix' = matrix_transpose x;
+/* Matrix conversions. */
+
private matrix_double matrix_complex matrix_int;
extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
expr* matrix_int(expr *x);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-20 05:42:41
|
Revision: 797
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=797&view=rev
Author: agraef
Date: 2008-09-20 05:42:23 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Allow trailing semicolon in row list of matrix values.
Modified Paths:
--------------
pure/trunk/parser.yy
Modified: pure/trunk/parser.yy
===================================================================
--- pure/trunk/parser.yy 2008-09-19 21:34:08 UTC (rev 796)
+++ pure/trunk/parser.yy 2008-09-20 05:42:23 UTC (rev 797)
@@ -555,6 +555,7 @@
rows
: row_list
+| row_list ';' { $$ = $1; }
| /* empty */ { $$ = new exprll; }
;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 22:26:27
|
Revision: 787
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=787&view=rev
Author: agraef
Date: 2008-09-19 05:26:38 +0000 (Fri, 19 Sep 2008)
Log Message:
-----------
Add matrix transposition and conversion operations.
Modified Paths:
--------------
pure/trunk/lib/prelude.pure
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/lib/prelude.pure 2008-09-19 05:26:38 UTC (rev 787)
@@ -58,6 +58,7 @@
infixl 6 + - or ; // addition, bitwise or
infixl 7 * / div mod and ; // multiplication, bitwise and
prefix 7 ~ ; // bitwise not
+postfix 7 ' ; // matrix transposition
infixr 8 ^ ; // exponentiation
prefix 8 # ; // size operator
infixl 9 ! !! ; // indexing, slicing
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/lib/primitives.pure 2008-09-19 05:26:38 UTC (rev 787)
@@ -428,6 +428,25 @@
when n::int,m::int = dim x end);
= throw out_of_bounds otherwise;
+private matrix_transpose;
+extern expr* matrix_transpose(expr *x);
+
+x::matrix' | x::cmatrix' | x::imatrix' = matrix_transpose x;
+
+private matrix_double matrix_complex matrix_int;
+extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
+ expr* matrix_int(expr *x);
+
+dmatrix x::matrix | dmatrix x::imatrix = matrix_double x;
+imatrix x::matrix | imatrix x::imatrix = matrix_int x;
+cmatrix x::matrix | cmatrix x::cmatrix | cmatrix x::imatrix = matrix_complex x;
+
+private matrix_re matrix_im;
+extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
+
+re x::matrix | re x::cmatrix | re x::imatrix = matrix_re x;
+im x::matrix | im x::cmatrix | im x::imatrix = matrix_im x;
+
/* IEEE floating point infinities and NaNs. Place these after the definitions
of the built-in operators so that the double arithmetic works. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/runtime.cc 2008-09-19 05:26:38 UTC (rev 787)
@@ -3789,6 +3789,201 @@
#endif
}
+extern "C"
+pure_expr *matrix_transpose(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(m, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[j*m2->tda+i] = m1->data[i*m1->tda+j];
+ return pure_double_matrix(m2);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_complex *m2 = create_complex_matrix(m, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m1->tda+j), l = 2*(j*m2->tda+i);
+ m2->data[l] = m1->data[k];
+ m2->data[l+1] = m1->data[k+1];
+ }
+ return pure_complex_matrix(m2);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(m, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[j*m2->tda+i] = m1->data[i*m1->tda+j];
+ return pure_int_matrix(m2);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_double(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX:
+ return x;
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[i*m2->tda+j] = (double)m1->data[i*m1->tda+j];
+ return pure_double_matrix(m2);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_complex(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_complex *m2 = create_complex_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m2->tda+j);
+ m2->data[k] = m1->data[i*m1->tda+j];
+ m2->data[k+1] = 0.0;
+ }
+ return pure_complex_matrix(m2);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_complex *m2 = create_complex_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m2->tda+j);
+ m2->data[k] = (double)m1->data[i*m1->tda+j];
+ m2->data[k+1] = 0.0;
+ }
+ return pure_complex_matrix(m2);
+ }
+ case EXPR::CMATRIX:
+ return x;
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_int(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[i*m2->tda+j] = (int)m1->data[i*m1->tda+j];
+ return pure_int_matrix(m2);
+ }
+ case EXPR::IMATRIX:
+ return x;
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_re(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m1->tda+j), l = i*m2->tda+j;
+ m2->data[l] = m1->data[k];
+ }
+ return pure_double_matrix(m2);
+ }
+ case EXPR::MATRIX:
+ case EXPR::IMATRIX:
+ return x;
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_im(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m1->tda+j), l = i*m2->tda+j;
+ m2->data[l] = m1->data[k+1];
+ }
+ return pure_double_matrix(m2);
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ memset(m2->data, 0, n*m*sizeof(double));
+ return pure_double_matrix(m2);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(n, m);
+ memset(m2->data, 0, n*m*sizeof(int));
+ return pure_int_matrix(m2);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
static uint32_t mpz_hash(const mpz_t z)
{
uint32_t h = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/runtime.h 2008-09-19 05:26:38 UTC (rev 787)
@@ -626,6 +626,27 @@
pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
uint32_t i2, uint32_t j2);
+/* Transpose a matrix. The resulting matrix has the rows of the original
+ matrix as its columns, and vice versa. */
+
+pure_expr *matrix_transpose(pure_expr *x);
+
+/* Convert an existing matrix to a double, complex or int matrix,
+ respectively. Any kind of matrix can be converted to a complex matrix, but
+ the input must be a double or integer matrix for the other conversions (see
+ matrix_re and matrix_im below to handle the complex->double case). */
+
+pure_expr *matrix_double(pure_expr *x);
+pure_expr *matrix_complex(pure_expr *x);
+pure_expr *matrix_int(pure_expr *x);
+
+/* Extract the real and imaginary parts of a matrix. If the input is a complex
+ matrix, the result is a double matrix. Otherwise the type of the result is
+ the same as that of the input matrix. */
+
+pure_expr *matrix_re(pure_expr *x);
+pure_expr *matrix_im(pure_expr *x);
+
/* Compute a 32 bit hash code of a Pure expression. This makes it possible to
use arbitary Pure values as keys in a hash table. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 13:35:50
|
Revision: 786
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=786&view=rev
Author: agraef
Date: 2008-09-18 20:35:54 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Add some basic matrix operations (type checking predicates for vectors).
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-18 20:19:31 UTC (rev 785)
+++ pure/trunk/lib/primitives.pure 2008-09-18 20:35:54 UTC (rev 786)
@@ -46,6 +46,22 @@
cmatrixp x = case x of _::cmatrix = 1; _ = 0 end;
imatrixp x = case x of _::imatrix = 1; _ = 0 end;
+/* Pure represents row and column vectors as matrices with 1 row or column,
+ respectively. The following predicates allow you to check for these special
+ kinds of matrices. */
+
+vectorp x = matrixp x && (n==1 || m==1 when n::int,m::int = dim x end);
+rowvectorp x = matrixp x && dim x!0==1;
+colvectorp x = matrixp x && dim x!1==1;
+
+cvectorp x = cmatrixp x && (n==1 || m==1 when n::int,m::int = dim x end);
+rowcvectorp x = cmatrixp x && dim x!0==1;
+colcvectorp x = cmatrixp x && dim x!1==1;
+
+ivectorp x = imatrixp x && (n==1 || m==1 when n::int,m::int = dim x end);
+rowivectorp x = imatrixp x && dim x!0==1;
+colivectorp x = imatrixp x && dim x!1==1;
+
/* Predicates to check for function objects, global (unbound) variables,
function applications, proper lists, list nodes and tuples. */
@@ -409,7 +425,7 @@
x::imatrix!(i::int,j::int)
= matrix_elem x i j
if (i>=0 && i<n && j>=0 && j<m
- when n,m = dim x end);
+ when n::int,m::int = dim x end);
= throw out_of_bounds otherwise;
/* IEEE floating point infinities and NaNs. Place these after the definitions
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 13:19:21
|
Revision: 785
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=785&view=rev
Author: agraef
Date: 2008-09-18 20:19:31 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Add some basic matrix operations (matrix size and dimensions, indexing, slicing).
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-18 17:35:01 UTC (rev 784)
+++ pure/trunk/ChangeLog 2008-09-18 20:19:31 UTC (rev 785)
@@ -1,17 +1,36 @@
2008-09-18 Albert Graef <Dr....@t-...>
+ * lib/primitives.pure: Add definitions of basic matrix operations.
+ Currently only size, dimensions and indexing are supported.
+
+ Synopsis:
+ - #x total number of elements
+ - dim x number of rows and columns (as a pair)
+ - x!i ith matrix element in row-major order
+ - x!(i,j) jth element in ith row of the matrix
+
* expr.cc, interpreter.cc, runtime.cc, printer.cc: Add basic GSL
matrix support. GSL double, complex and integer matrices can be
created with the new {x,y;u,v} syntax, which works more or less
like Octave/MATLAB matrices, but using curly braces instead of
- brackets.
+ brackets. Moreover, various basic operations to handle this kind
+ of objects (conversions, determining sizes, indexing, slicing)
+ have been added to the runtime, including some public API
+ operations to create and inspect matrix objects.
+ Note that the {...} syntax can be used only on the right-hand side
+ of a definition, matrix patterns are not supported right now. As a
+ remedy, there are three new type tags, matrix, cmatrix and
+ imatrix, which can be used in patterns to match double, complex
+ and integer matrices, respectively.
+
GSL matrices are always homogeneous, i.e., they only contain
values from one numeric type. Integer matrices can be created from
any combination of Pure machine int and bigint values (the latter
are converted to machine ints automatically). Matrices with at
least one double or complex element become double and complex
- matrices, respectively.
+ matrices, respectively, casting the other matrix elements as
+ needed.
Complex matrices can be created from either pairs of double or
integer values, such as {(1,2),(3,4)}, or from Pure complex
@@ -20,18 +39,10 @@
math.pure has been loaded in which case complex values are printed
in rectangular format x+:y. Also note that, for performance
reasons, the expression printer doesn't use the __show__ function
- to print matrix elements, but of course it is possible to override
- the default print representation of matrix values as a whole.
+ to print individual matrix elements, but of course it is possible
+ to override the default print representation of matrix values as a
+ whole.
- Note that the {...} syntax can be used only on the right-hand side
- of a definition, matrix patterns are not supported right now. As a
- remedy, there are three new type tags, matrix, cmatrix and
- imatrix, which can be used in patterns to match double, complex
- and integer matrices, respectively.
-
- Operations to handle Pure matrix expressions have been added to
- the public runtime API.
-
2008-09-15 Albert Graef <Dr....@t-...>
* configure.ac: Bump version number.
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-18 17:35:01 UTC (rev 784)
+++ pure/trunk/lib/primitives.pure 2008-09-18 20:19:31 UTC (rev 785)
@@ -42,6 +42,10 @@
stringp x = case x of _::string = 1; _ = 0 end;
pointerp x = case x of _::pointer = 1; _ = 0 end;
+matrixp x = case x of _::matrix = 1; _ = 0 end;
+cmatrixp x = case x of _::cmatrix = 1; _ = 0 end;
+imatrixp x = case x of _::imatrix = 1; _ = 0 end;
+
/* Predicates to check for function objects, global (unbound) variables,
function applications, proper lists, list nodes and tuples. */
@@ -383,6 +387,31 @@
x::pointer==y::pointer = bigint x == bigint y;
x::pointer!=y::pointer = bigint x != bigint y;
+/* Basic matrix operations. */
+
+private matrix_size matrix_dim;
+extern int matrix_size(expr *x), expr* matrix_dim(expr *x);
+
+#x::matrix | #x::cmatrix | #x::imatrix = matrix_size x;
+dim x::matrix | dim x::cmatrix | dim x::imatrix = matrix_dim x;
+
+private matrix_elem_at;
+extern expr* matrix_elem_at(expr* x, int i);
+
+x::matrix!i::int | x::cmatrix!i::int | x::imatrix!i::int
+ = matrix_elem_at x i if i>=0 && i<#x;
+ = throw out_of_bounds otherwise;
+
+private matrix_elem;
+extern expr* matrix_elem(expr* x, int i, int j);
+
+x::matrix!(i::int,j::int) | x::cmatrix!(i::int,j::int) |
+x::imatrix!(i::int,j::int)
+ = matrix_elem x i j
+ if (i>=0 && i<n && j>=0 && j<m
+ when n,m = dim x end);
+ = throw out_of_bounds otherwise;
+
/* IEEE floating point infinities and NaNs. Place these after the definitions
of the built-in operators so that the double arithmetic works. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-18 17:35:01 UTC (rev 784)
+++ pure/trunk/runtime.cc 2008-09-18 20:19:31 UTC (rev 785)
@@ -3611,6 +3611,184 @@
return interp.errmsg.c_str();
}
+extern "C"
+uint32_t matrix_size(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return m->size1*m->size2;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return m->size1*m->size2;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return m->size1*m->size2;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_dim(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return pure_tuplel(2, pure_int(m->size1), pure_int(m->size2));
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return pure_tuplel(2, pure_int(m->size1), pure_int(m->size2));
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return pure_tuplel(2, pure_int(m->size1), pure_int(m->size2));
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+static inline pure_expr *make_complex(double a, double b)
+{
+ interpreter& interp = *interpreter::g_interp;
+ symbol *rect = interp.symtab.complex_rect_sym();
+ if (rect)
+ return pure_appl(pure_symbol(rect->f), 2, pure_double(a), pure_double(b));
+ else
+ return pure_tuplel(2, pure_double(a), pure_double(b));
+}
+
+extern "C"
+pure_expr *matrix_elem_at(pure_expr *x, uint32_t i)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return pure_double(m->data[i]);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return make_complex(m->data[2*i], m->data[2*i+1]);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return pure_int(m->data[i]);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_elem(pure_expr *x, uint32_t i, uint32_t j)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ size_t k = i*m->tda+j;
+ return pure_double(m->data[k]);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ size_t k = 2*(i*m->tda+j);
+ return make_complex(m->data[k], m->data[k+1]);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ size_t k = i*m->tda+j;
+ return pure_int(m->data[k]);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
+ uint32_t i2, uint32_t j2)
+{
+#ifdef HAVE_GSL
+ void *p = 0;
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ gsl_matrix_view v = gsl_matrix_submatrix(m, i1, j1, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m1 = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ assert(m1 && v.matrix.data);
+ *m1 = v.matrix;
+ p = m1;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_submatrix(m, i1, j1, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m1 =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ assert(m1 && v.matrix.data);
+ *m1 = v.matrix;
+ p = m1;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ gsl_matrix_int_view v = gsl_matrix_int_submatrix(m, i1, j1, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m1 = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ assert(m1 && v.matrix.data);
+ *m1 = v.matrix;
+ p = m1;
+ break;
+ }
+ default:
+ return 0;
+ }
+ // create a new expression for the slice, update the reference counter for
+ // the underlying GSL matrix
+ pure_expr *y = new_expr();
+ y->tag = EXPR::MATRIX;
+ y->data.mat.p = p;
+ y->data.mat.refc = x->data.mat.refc;
+ *y->data.mat.refc++;
+ MEMDEBUG_NEW(y)
+ return y;
+#else
+ return 0;
+#endif
+}
+
static uint32_t mpz_hash(const mpz_t z)
{
uint32_t h = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-18 17:35:01 UTC (rev 784)
+++ pure/trunk/runtime.h 2008-09-18 20:19:31 UTC (rev 785)
@@ -602,6 +602,30 @@
const char *lasterr();
+/* Basic matrix operations. These work with all supported GSL matrix types.
+ matrix_size determines the number of elements in a matrix, matrix_dim the
+ number of rows and columns, which are returned as a pair (n,m). */
+
+uint32_t matrix_size(pure_expr *x);
+pure_expr *matrix_dim(pure_expr *x);
+
+/* Matrix elements can be retrieved either by a single index (using row-major
+ order), or by row and column index. All indices are zero-based. Indices
+ aren't range-checked, if this is needed you have to do it beforehand using
+ matrix_size or matrix_dim above. */
+
+pure_expr *matrix_elem_at(pure_expr *x, uint32_t i);
+pure_expr *matrix_elem(pure_expr *x, uint32_t i, uint32_t j);
+
+/* The following operation retrieves a slice a.k.a. submatrix of a matrix and
+ returns it as a matrix object. The new matrix object shares the underlying
+ storage with the original matrix (i.e., matrix elements are *not* copied)
+ and so this is a comparatively cheap operation. Indices are zero-based and
+ not checked. */
+
+pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
+ uint32_t i2, uint32_t j2);
+
/* Compute a 32 bit hash code of a Pure expression. This makes it possible to
use arbitary Pure values as keys in a hash table. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 11:44:19
|
Revision: 782
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=782&view=rev
Author: agraef
Date: 2008-09-18 11:44:30 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Implement the remaining matrix functions in the public runtime API. This completes basic GSL matrix support, but note that not all functionality has been fully tested yet, and operations to inspect matrix values are still missing.
Modified Paths:
--------------
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-18 11:41:11 UTC (rev 781)
+++ pure/trunk/runtime.cc 2008-09-18 11:44:30 UTC (rev 782)
@@ -1128,19 +1128,97 @@
pure_expr *pure_matrix_rowsl(uint32_t n, ...)
{
#ifdef HAVE_GSL
- // XXXTODO
- return 0;
+ va_list ap;
+ pure_expr **xs = (pure_expr**)alloca(n*sizeof(pure_expr*));
+ va_start(ap, n);
+ for (size_t i = 0; i < n; i++)
+ xs[i] = va_arg(ap, pure_expr*);
+ va_end(ap);
+ return pure_matrix_rowsv(n, xs);
#else
return 0;
#endif
}
extern "C"
-pure_expr *pure_matrix_rowsv(uint32_t n, pure_expr **elems)
+pure_expr *pure_matrix_rowsv(uint32_t n, pure_expr **xs)
{
#ifdef HAVE_GSL
- // XXXTODO
- return 0;
+ int k = -1;
+ size_t nrows = 0, ncols = 0;
+ int32_t target = EXPR::IMATRIX;
+ bool have_matrix = false;
+ for (size_t i = 0; i < n; i++) {
+ pure_expr *x = xs[i];
+ switch (x->tag) {
+ case EXPR::DBL:
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ case EXPR::INT:
+ case EXPR::BIGINT:
+ if (k >= 0 && k != 1) return 0;
+ nrows++; k = 1;
+ break;
+ case EXPR::APP: {
+ double a, b;
+ if (!get_complex(x, a, b)) return 0;
+ if (k >= 0 && k != 1) return 0;
+ target = EXPR::CMATRIX;
+ nrows++; k = 1;
+ break;
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ return 0;
+ nrows += mp->size1; k = mp->size2;
+ } else if (k>0)
+ return 0;
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ have_matrix = true;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ return 0;
+ nrows += mp->size1; k = mp->size2;
+ } else if (k>0)
+ return 0;
+ target = EXPR::CMATRIX;
+ have_matrix = true;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ return 0;
+ nrows += mp->size1; k = mp->size2;
+ } else if (k>0)
+ return 0;
+ have_matrix = true;
+ break;
+ }
+ default:
+ return 0;
+ }
+ }
+ if (n == 1 && have_matrix) return xs[0];
+ if (k < 0) k = 0;
+ ncols = k;
+ if (nrows == 0 && ncols == 0) target = EXPR::MATRIX;
+ switch (target) {
+ case EXPR::MATRIX:
+ return double_matrix_rows(nrows, ncols, n, xs);
+ case EXPR::CMATRIX:
+ return complex_matrix_rows(nrows, ncols, n, xs);
+ case EXPR::IMATRIX:
+ return int_matrix_rows(nrows, ncols, n, xs);
+ default:
+ return 0;
+ }
#else
return 0;
#endif
@@ -1150,19 +1228,97 @@
pure_expr *pure_matrix_columnsl(uint32_t n, ...)
{
#ifdef HAVE_GSL
- // XXXTODO
- return 0;
+ va_list ap;
+ pure_expr **xs = (pure_expr**)alloca(n*sizeof(pure_expr*));
+ va_start(ap, n);
+ for (size_t i = 0; i < n; i++)
+ xs[i] = va_arg(ap, pure_expr*);
+ va_end(ap);
+ return pure_matrix_columnsv(n, xs);
#else
return 0;
#endif
}
extern "C"
-pure_expr *pure_matrix_columnsv(uint32_t n, pure_expr **elems)
+pure_expr *pure_matrix_columnsv(uint32_t n, pure_expr **xs)
{
#ifdef HAVE_GSL
- // XXXTODO
- return 0;
+ int k = -1;
+ size_t nrows = 0, ncols = 0;
+ int32_t target = EXPR::IMATRIX;
+ bool have_matrix = false;
+ for (size_t i = 0; i < n; i++) {
+ pure_expr *x = xs[i];
+ switch (x->tag) {
+ case EXPR::DBL:
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ case EXPR::INT:
+ case EXPR::BIGINT:
+ if (k >= 0 && k != 1) return 0;
+ ncols++; k = 1;
+ break;
+ case EXPR::APP: {
+ double a, b;
+ if (!get_complex(x, a, b)) return 0;
+ if (k >= 0 && k != 1) return 0;
+ target = EXPR::CMATRIX;
+ ncols++; k = 1;
+ break;
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ return 0;
+ ncols += mp->size2; k = mp->size1;
+ } else if (k>0)
+ return 0;
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ have_matrix = true;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ return 0;
+ ncols += mp->size2; k = mp->size1;
+ } else if (k>0)
+ return 0;
+ target = EXPR::CMATRIX;
+ have_matrix = true;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ return 0;
+ ncols += mp->size2; k = mp->size1;
+ } else if (k>0)
+ return 0;
+ have_matrix = true;
+ break;
+ }
+ default:
+ return 0;
+ }
+ }
+ if (n == 1 && have_matrix) return xs[0];
+ if (k < 0) k = 0;
+ nrows = k;
+ if (nrows == 0 && ncols == 0) target = EXPR::MATRIX;
+ switch (target) {
+ case EXPR::MATRIX:
+ return double_matrix_columns(nrows, ncols, n, xs);
+ case EXPR::CMATRIX:
+ return complex_matrix_columns(nrows, ncols, n, xs);
+ case EXPR::IMATRIX:
+ return int_matrix_columns(nrows, ncols, n, xs);
+ default:
+ return 0;
+ }
#else
return 0;
#endif
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-18 11:41:11 UTC (rev 781)
+++ pure/trunk/runtime.h 2008-09-18 11:44:30 UTC (rev 782)
@@ -149,10 +149,10 @@
The pure_matrix_rows functions arrange the elements vertically, while the
pure_matrix_columns functions arrange them horizontally, given that the
other dimensions match. A null expression is returned in case of an error
- (no matrix support, dimension mismatch, or invalid element type). Otherwise
- a new matrix expression is returned. Temporary element expressions are
- taken over by the callee and will be garbage-collected, but the elems
- vectors are owned by the caller and won't be freed. */
+ (no matrix support, dimension mismatch, or invalid element type), leaving
+ the input expressions untouched. Otherwise a new matrix expression is
+ returned and temporary element expressions are garbage-collected. In any
+ case, the elems vectors are owned by the caller and won't be freed. */
pure_expr *pure_matrix_rowsl(uint32_t n, ...);
pure_expr *pure_matrix_rowsv(uint32_t n, pure_expr **elems);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 10:56:32
|
Revision: 780
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=780&view=rev
Author: agraef
Date: 2008-09-18 10:56:42 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Finish matrix construction operations (complex matrices).
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-18 09:57:08 UTC (rev 779)
+++ pure/trunk/runtime.cc 2008-09-18 10:56:42 UTC (rev 780)
@@ -28,6 +28,7 @@
#include <unistd.h>
#include <limits.h>
#include <locale.h>
+#include <math.h>
#include <iostream>
#include <sstream>
@@ -827,7 +828,209 @@
return pure_double_matrix(mat);
}
+static inline bool get_complex(pure_expr *x, double& a, double& b)
+{
+ if (x->tag != EXPR::APP) return false;
+ pure_expr *u = x->data.x[0], *v = x->data.x[1];
+ if (u->tag == EXPR::APP) {
+ interpreter& interp = *interpreter::g_interp;
+ pure_expr *f = u->data.x[0];
+ symbol *rect = interp.symtab.complex_rect_sym(),
+ *polar = interp.symtab.complex_polar_sym();
+ if ((!rect || f->tag != rect->f) &&
+ (!polar || f->tag != polar->f) &&
+ f->tag != interp.symtab.pair_sym().f)
+ return false;
+ u = u->data.x[1];
+ switch (u->tag) {
+ case EXPR::INT:
+ a = (double)u->data.i;
+ break;
+ case EXPR::BIGINT:
+ a = mpz_get_d(x->data.z);
+ break;
+ case EXPR::DBL:
+ a = u->data.d;
+ break;
+ default:
+ return false;
+ }
+ switch (v->tag) {
+ case EXPR::INT:
+ b = (double)v->data.i;
+ break;
+ case EXPR::BIGINT:
+ b = mpz_get_d(x->data.z);
+ break;
+ case EXPR::DBL:
+ b = v->data.d;
+ break;
+ default:
+ return false;
+ }
+ if (polar && f->tag == polar->f) {
+ double r = a, t = b;
+ a = r*cos(t); b = r*sin(t);
+ }
+ return true;
+ } else
+ return false;
+}
+
static pure_expr*
+complex_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
+{
+ gsl_matrix_complex *mat = create_complex_matrix(nrows, ncols);
+ if (!mat) return 0;
+ double *data = mat->data;
+ size_t tda = mat->tda;
+ for (size_t count = 0, i = 0; count < n; count++) {
+ pure_expr *x = xs[count];
+ pure_new_internal(x);
+ switch (x->tag) {
+ case EXPR::INT:
+ data[2*i*tda] = (double)x->data.i;
+ data[2*i*tda+1] = 0.0;
+ i++;
+ break;
+ case EXPR::BIGINT:
+ data[2*i*tda] = mpz_get_d(x->data.z);
+ data[2*i*tda+1] = 0.0;
+ i++;
+ break;
+ case EXPR::DBL:
+ data[2*i*tda] = x->data.d;
+ data[2*i*tda+1] = 0.0;
+ i++;
+ break;
+ case EXPR::APP: {
+ double a, b;
+ if (get_complex(x, a, b)) {
+ data[2*i*tda] = a;
+ data[2*i*tda+1] = b;
+ i++;
+ } else {
+ assert(0 && "bad matrix element");
+ }
+ break;
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ for (size_t k = 0; k < mat1->size2; k++) {
+ data[2*(i*tda+k)] = mat1->data[j*mat1->tda+k];
+ data[2*(i*tda+k)+1] = 0.0;
+ }
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ for (size_t k = 0; k < mat1->size2; k++) {
+ data[2*(i*tda+k)] = (double)mat1->data[j*mat1->tda+k];
+ data[2*(i*tda+k)+1] = 0.0;
+ }
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mat1 = (gsl_matrix_complex*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ memcpy(data+2*i*tda, mat1->data+2*j*mat1->tda,
+ ncols*2*sizeof(double));
+ break;
+ }
+ default:
+ assert(0 && "bad matrix element");
+ break;
+ }
+ }
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ return pure_complex_matrix(mat);
+}
+
+static pure_expr*
+complex_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
+{
+ gsl_matrix_complex *mat = create_complex_matrix(nrows, ncols);
+ if (!mat) return 0;
+ double *data = mat->data;
+ size_t tda = mat->tda;
+ for (size_t count = 0, i = 0; count < n; count++) {
+ pure_expr *x = xs[count];
+ pure_new_internal(x);
+ switch (x->tag) {
+ case EXPR::INT:
+ data[2*i] = (double)x->data.i;
+ data[2*i+1] = 0.0;
+ i++;
+ break;
+ case EXPR::BIGINT:
+ data[2*i] = mpz_get_d(x->data.z);
+ data[2*i+1] = 0.0;
+ i++;
+ break;
+ case EXPR::DBL:
+ data[2*i] = x->data.d;
+ data[2*i+1] = 0.0;
+ i++;
+ break;
+ case EXPR::APP: {
+ double a, b;
+ if (get_complex(x, a, b)) {
+ data[2*i] = a;
+ data[2*i+1] = b;
+ i++;
+ } else {
+ assert(0 && "bad matrix element");
+ }
+ break;
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; j++)
+ for (size_t k = 0; k < mat1->size2; k++) {
+ data[2*(j*tda+k+i)] = mat1->data[j*mat1->tda+k];
+ data[2*(j*tda+k+i)+1] = 0.0;
+ }
+ i += mat1->size2;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; j++)
+ for (size_t k = 0; k < mat1->size2; k++) {
+ data[2*(j*tda+k+i)] = (double)mat1->data[j*mat1->tda+k];
+ data[2*(j*tda+k+i)+1] = 0.0;
+ }
+ i += mat1->size2;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mat1 = (gsl_matrix_complex*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; j++)
+ memcpy(data+2*(j*tda+i), mat1->data+2*j*mat1->tda,
+ mat1->size2*2*sizeof(double));
+ i += mat1->size2;
+ break;
+ }
+ default:
+ assert(0 && "bad matrix element");
+ break;
+ }
+ }
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ return pure_complex_matrix(mat);
+}
+
+static pure_expr*
int_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
gsl_matrix_int *mat = create_int_matrix(nrows, ncols);
@@ -1845,25 +2048,8 @@
nrows++; k = 1;
break;
case EXPR::APP: {
- pure_expr *u = x->data.x[0], *v = x->data.x[1];
- if (u->tag == EXPR::APP) {
- interpreter& interp = *interpreter::g_interp;
- pure_expr *f = u->data.x[0];
- symbol *rect = interp.symtab.complex_rect_sym(),
- *polar = interp.symtab.complex_polar_sym();
- if ((!rect || f->tag != rect->f) &&
- (!polar || f->tag != polar->f) &&
- f->tag != interp.symtab.pair_sym().f)
- goto err;
- u = u->data.x[1];
- if (u->tag != EXPR::INT && u->tag != EXPR::BIGINT &&
- u->tag != EXPR::DBL)
- goto err;
- if (v->tag != EXPR::INT && v->tag != EXPR::BIGINT &&
- v->tag != EXPR::DBL)
- goto err;
- } else
- goto err;
+ double a, b;
+ if (!get_complex(x, a, b)) goto err;
if (k >= 0 && k != 1) goto err;
target = EXPR::CMATRIX;
nrows++; k = 1;
@@ -1917,11 +2103,9 @@
case EXPR::MATRIX:
ret = double_matrix_rows(nrows, ncols, n, xs);
break;
-#if 0 // XXXTODO
case EXPR::CMATRIX:
ret = complex_matrix_rows(nrows, ncols, n, xs);
break;
-#endif
case EXPR::IMATRIX:
ret = int_matrix_rows(nrows, ncols, n, xs);
break;
@@ -1973,25 +2157,8 @@
ncols++; k = 1;
break;
case EXPR::APP: {
- pure_expr *u = x->data.x[0], *v = x->data.x[1];
- if (u->tag == EXPR::APP) {
- interpreter& interp = *interpreter::g_interp;
- pure_expr *f = u->data.x[0];
- symbol *rect = interp.symtab.complex_rect_sym(),
- *polar = interp.symtab.complex_polar_sym();
- if ((!rect || f->tag != rect->f) &&
- (!polar || f->tag != polar->f) &&
- f->tag != interp.symtab.pair_sym().f)
- goto err;
- u = u->data.x[1];
- if (u->tag != EXPR::INT && u->tag != EXPR::BIGINT &&
- u->tag != EXPR::DBL)
- goto err;
- if (v->tag != EXPR::INT && v->tag != EXPR::BIGINT &&
- v->tag != EXPR::DBL)
- goto err;
- } else
- goto err;
+ double a, b;
+ if (!get_complex(x, a, b)) goto err;
if (k >= 0 && k != 1) goto err;
target = EXPR::CMATRIX;
ncols++; k = 1;
@@ -2045,11 +2212,9 @@
case EXPR::MATRIX:
ret = double_matrix_columns(nrows, ncols, n, xs);
break;
-#if 0 // XXXTODO
case EXPR::CMATRIX:
ret = complex_matrix_columns(nrows, ncols, n, xs);
break;
-#endif
case EXPR::IMATRIX:
ret = int_matrix_columns(nrows, ncols, n, xs);
break;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 10:34:50
|
Revision: 784
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=784&view=rev
Author: agraef
Date: 2008-09-18 17:35:01 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-18 12:13:09 UTC (rev 783)
+++ pure/trunk/runtime.cc 2008-09-18 17:35:01 UTC (rev 784)
@@ -638,6 +638,7 @@
#ifdef HAVE_GSL
gsl_matrix *m = (gsl_matrix*)p;
if (!m || !m->owner) return 0;
+ m->owner = 0;
pure_expr *x = new_expr();
x->tag = EXPR::MATRIX;
x->data.mat.p = p;
@@ -656,6 +657,7 @@
#ifdef HAVE_GSL
gsl_matrix_complex *m = (gsl_matrix_complex*)p;
if (!m || !m->owner) return 0;
+ m->owner = 0;
pure_expr *x = new_expr();
x->tag = EXPR::CMATRIX;
x->data.mat.p = p;
@@ -674,6 +676,7 @@
#ifdef HAVE_GSL
gsl_matrix_int *m = (gsl_matrix_int*)p;
if (!m || !m->owner) return 0;
+ m->owner = 0;
pure_expr *x = new_expr();
x->tag = EXPR::IMATRIX;
x->data.mat.p = p;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-18 09:56:57
|
Revision: 779
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=779&view=rev
Author: agraef
Date: 2008-09-18 09:57:08 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Add support for empty matrices.
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-17 21:17:01 UTC (rev 778)
+++ pure/trunk/runtime.cc 2008-09-18 09:57:08 UTC (rev 779)
@@ -585,6 +585,52 @@
return x;
}
+#ifdef HAVE_GSL
+/* GSL doesn't really support empty matrices, so we fake them by allocating 1
+ dummy row or column if the corresponding dimension is actually zero. */
+static inline gsl_matrix*
+create_double_matrix(size_t nrows, size_t ncols)
+{
+ if (nrows == 0 || ncols == 0 ) {
+ size_t nrows1 = (nrows>0)?nrows:1;
+ size_t ncols1 = (ncols>0)?ncols:1;
+ gsl_matrix *m = gsl_matrix_calloc(nrows1, ncols1);
+ if (!m) return 0;
+ m->size1 = nrows; m->size2 = ncols;
+ return m;
+ } else
+ return gsl_matrix_alloc(nrows, ncols);
+}
+
+static inline gsl_matrix_complex*
+create_complex_matrix(size_t nrows, size_t ncols)
+{
+ if (nrows == 0 || ncols == 0 ) {
+ size_t nrows1 = (nrows>0)?nrows:1;
+ size_t ncols1 = (ncols>0)?ncols:1;
+ gsl_matrix_complex *m = gsl_matrix_complex_calloc(nrows1, ncols1);
+ if (!m) return 0;
+ m->size1 = nrows; m->size2 = ncols;
+ return m;
+ } else
+ return gsl_matrix_complex_alloc(nrows, ncols);
+}
+
+static inline gsl_matrix_int*
+create_int_matrix(size_t nrows, size_t ncols)
+{
+ if (nrows == 0 || ncols == 0 ) {
+ size_t nrows1 = (nrows>0)?nrows:1;
+ size_t ncols1 = (ncols>0)?ncols:1;
+ gsl_matrix_int *m = gsl_matrix_int_calloc(nrows1, ncols1);
+ if (!m) return 0;
+ m->size1 = nrows; m->size2 = ncols;
+ return m;
+ } else
+ return gsl_matrix_int_alloc(nrows, ncols);
+}
+#endif
+
extern "C"
pure_expr *pure_double_matrix(void *p)
{
@@ -645,8 +691,10 @@
#ifdef HAVE_GSL
gsl_matrix *m1 = (gsl_matrix*)p;
if (!m1) return 0;
- gsl_matrix *m2 = gsl_matrix_alloc(m1->size1, m1->size2);
- gsl_matrix_memcpy(m2, m1);
+ gsl_matrix *m2 = create_double_matrix(m1->size1, m1->size2);
+ if (!m2) return 0;
+ if (m1->size1 > 0 && m1->size2 > 0)
+ gsl_matrix_memcpy(m2, m1);
return pure_double_matrix(m2);
#else
return 0;
@@ -659,8 +707,10 @@
#ifdef HAVE_GSL
gsl_matrix_complex *m1 = (gsl_matrix_complex*)p;
if (!m1) return 0;
- gsl_matrix_complex *m2 = gsl_matrix_complex_alloc(m1->size1, m1->size2);
- gsl_matrix_complex_memcpy(m2, m1);
+ gsl_matrix_complex *m2 = create_complex_matrix(m1->size1, m1->size2);
+ if (!m2) return 0;
+ if (m1->size1 > 0 && m1->size2 > 0)
+ gsl_matrix_complex_memcpy(m2, m1);
return pure_complex_matrix(m2);
#else
return 0;
@@ -673,8 +723,10 @@
#ifdef HAVE_GSL
gsl_matrix_int *m1 = (gsl_matrix_int*)p;
if (!m1) return 0;
- gsl_matrix_int *m2 = gsl_matrix_int_alloc(m1->size1, m1->size2);
- gsl_matrix_int_memcpy(m2, m1);
+ gsl_matrix_int *m2 = create_int_matrix(m1->size1, m1->size2);
+ if (!m2) return 0;
+ if (m1->size1 > 0 && m1->size2 > 0)
+ gsl_matrix_int_memcpy(m2, m1);
return pure_int_matrix(m2);
#else
return 0;
@@ -685,15 +737,8 @@
static pure_expr*
double_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
- gsl_matrix *mat = gsl_matrix_alloc(nrows, ncols);
- if (!mat) {
- // XXXTODO: empty matrices
- for (size_t i = 0; i < n; i++)
- pure_new_internal(xs[i]);
- for (size_t i = 0; i < n; i++)
- pure_free_internal(xs[i]);
- return 0;
- }
+ gsl_matrix *mat = create_double_matrix(nrows, ncols);
+ if (!mat) return 0;
double *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
@@ -737,15 +782,8 @@
static pure_expr*
double_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
- gsl_matrix *mat = gsl_matrix_alloc(nrows, ncols);
- if (!mat) {
- // XXXTODO: empty matrices
- for (size_t i = 0; i < n; i++)
- pure_new_internal(xs[i]);
- for (size_t i = 0; i < n; i++)
- pure_free_internal(xs[i]);
- return 0;
- }
+ gsl_matrix *mat = create_double_matrix(nrows, ncols);
+ if (!mat) return 0;
double *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
@@ -792,15 +830,8 @@
static pure_expr*
int_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
- gsl_matrix_int *mat = gsl_matrix_int_alloc(nrows, ncols);
- if (!mat) {
- // XXXTODO: empty matrices
- for (size_t i = 0; i < n; i++)
- pure_new_internal(xs[i]);
- for (size_t i = 0; i < n; i++)
- pure_free_internal(xs[i]);
- return 0;
- }
+ gsl_matrix_int *mat = create_int_matrix(nrows, ncols);
+ if (!mat) return 0;
int *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
@@ -844,15 +875,8 @@
static pure_expr*
int_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
- gsl_matrix_int *mat = gsl_matrix_int_alloc(nrows, ncols);
- if (!mat) {
- // XXXTODO: empty matrices
- for (size_t i = 0; i < n; i++)
- pure_new_internal(xs[i]);
- for (size_t i = 0; i < n; i++)
- pure_free_internal(xs[i]);
- return 0;
- }
+ gsl_matrix_int *mat = create_int_matrix(nrows, ncols);
+ if (!mat) return 0;
int *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-17 21:16:50
|
Revision: 778
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=778&view=rev
Author: agraef
Date: 2008-09-17 21:17:01 +0000 (Wed, 17 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/runtime.cc
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-17 15:05:12 UTC (rev 777)
+++ pure/trunk/runtime.cc 2008-09-17 21:17:01 UTC (rev 778)
@@ -817,7 +817,7 @@
data[i++*tda] = (int)x->data.d;
break;
case EXPR::MATRIX: {
- gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
if (mat1)
for (size_t j = 0; j < mat1->size1; i++, j++)
for (size_t k = 0; k < mat1->size2; k++)
@@ -825,7 +825,7 @@
break;
}
case EXPR::IMATRIX: {
- gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
if (mat1)
for (size_t j = 0; j < mat1->size1; i++, j++)
memcpy(data+i*tda, mat1->data+j*mat1->tda, ncols*sizeof(int));
@@ -869,7 +869,7 @@
data[i++] = (int)x->data.d;
break;
case EXPR::MATRIX: {
- gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
if (mat1)
for (size_t j = 0; j < mat1->size1; j++)
for (size_t k = 0; k < mat1->size2; k++)
@@ -878,7 +878,7 @@
break;
}
case EXPR::IMATRIX: {
- gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
if (mat1)
for (size_t j = 0; j < mat1->size1; j++)
memcpy(data+j*tda+i, mat1->data+j*mat1->tda,
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-17 15:26:10
|
Revision: 777
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=777&view=rev
Author: agraef
Date: 2008-09-17 15:05:12 +0000 (Wed, 17 Sep 2008)
Log Message:
-----------
Bugfixes.
Modified Paths:
--------------
pure/trunk/interpreter.cc
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-17 08:32:48 UTC (rev 776)
+++ pure/trunk/interpreter.cc 2008-09-17 15:05:12 UTC (rev 777)
@@ -5033,10 +5033,7 @@
us[i] =
act_env().CreateCall(module->getFunction("pure_matrix_columns"), vs);
}
- if (n == 1)
- return us[1];
- else
- return act_env().CreateCall(module->getFunction("pure_matrix_rows"), us);
+ return act_env().CreateCall(module->getFunction("pure_matrix_rows"), us);
}
// application:
case EXPR::APP:
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-17 08:32:48 UTC (rev 776)
+++ pure/trunk/runtime.cc 2008-09-17 15:05:12 UTC (rev 777)
@@ -686,11 +686,19 @@
double_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
gsl_matrix *mat = gsl_matrix_alloc(nrows, ncols);
- if (!mat) return 0; // XXXTODO: empty matrices
+ if (!mat) {
+ // XXXTODO: empty matrices
+ for (size_t i = 0; i < n; i++)
+ pure_new_internal(xs[i]);
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ return 0;
+ }
double *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
pure_expr *x = xs[count];
+ pure_new_internal(x);
switch (x->tag) {
case EXPR::INT:
data[i++*tda] = (double)x->data.i;
@@ -720,8 +728,9 @@
assert(0 && "bad matrix element");
break;
}
- pure_freenew(x);
}
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
return pure_double_matrix(mat);
}
@@ -729,11 +738,19 @@
double_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
gsl_matrix *mat = gsl_matrix_alloc(nrows, ncols);
- if (!mat) return 0; // XXXTODO: empty matrices
+ if (!mat) {
+ // XXXTODO: empty matrices
+ for (size_t i = 0; i < n; i++)
+ pure_new_internal(xs[i]);
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ return 0;
+ }
double *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
pure_expr *x = xs[count];
+ pure_new_internal(x);
switch (x->tag) {
case EXPR::INT:
data[i++] = (double)x->data.i;
@@ -747,8 +764,9 @@
case EXPR::MATRIX: {
gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
if (mat1)
- for (size_t j = 0; j < mat1->size1; i++, j++)
- memcpy(data+j*tda+i, mat1->data+j*mat1->tda, ncols*sizeof(double));
+ for (size_t j = 0; j < mat1->size1; j++)
+ memcpy(data+j*tda+i, mat1->data+j*mat1->tda,
+ mat1->size2*sizeof(double));
i += mat1->size2;
break;
}
@@ -765,8 +783,9 @@
assert(0 && "bad matrix element");
break;
}
- pure_freenew(x);
}
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
return pure_double_matrix(mat);
}
@@ -774,11 +793,19 @@
int_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
gsl_matrix_int *mat = gsl_matrix_int_alloc(nrows, ncols);
- if (!mat) return 0; // XXXTODO: empty matrices
+ if (!mat) {
+ // XXXTODO: empty matrices
+ for (size_t i = 0; i < n; i++)
+ pure_new_internal(xs[i]);
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ return 0;
+ }
int *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
pure_expr *x = xs[count];
+ pure_new_internal(x);
switch (x->tag) {
case EXPR::INT:
data[i++*tda] = x->data.i;
@@ -808,8 +835,9 @@
assert(0 && "bad matrix element");
break;
}
- pure_freenew(x);
}
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
return pure_int_matrix(mat);
}
@@ -817,11 +845,19 @@
int_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
{
gsl_matrix_int *mat = gsl_matrix_int_alloc(nrows, ncols);
- if (!mat) return 0; // XXXTODO: empty matrices
+ if (!mat) {
+ // XXXTODO: empty matrices
+ for (size_t i = 0; i < n; i++)
+ pure_new_internal(xs[i]);
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ return 0;
+ }
int *data = mat->data;
size_t tda = mat->tda;
for (size_t count = 0, i = 0; count < n; count++) {
pure_expr *x = xs[count];
+ pure_new_internal(x);
switch (x->tag) {
case EXPR::INT:
data[i++] = x->data.i;
@@ -844,8 +880,9 @@
case EXPR::IMATRIX: {
gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
if (mat1)
- for (size_t j = 0; j < mat1->size1; i++, j++)
- memcpy(data+j*tda+i, mat1->data+j*mat1->tda, ncols*sizeof(int));
+ for (size_t j = 0; j < mat1->size1; j++)
+ memcpy(data+j*tda+i, mat1->data+j*mat1->tda,
+ mat1->size2*sizeof(int));
i += mat1->size2;
break;
}
@@ -853,8 +890,9 @@
assert(0 && "bad matrix element");
break;
}
- pure_freenew(x);
}
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
return pure_int_matrix(mat);
}
#endif
@@ -1762,21 +1800,25 @@
{
#ifdef HAVE_GSL
va_list ap;
- va_start(ap, n);
pure_expr **xs = (pure_expr**)alloca(n*sizeof(pure_expr*));
- size_t m = 0; int k = -1;
+ int k = -1;
size_t nrows = 0, ncols = 0;
int32_t target = EXPR::IMATRIX;
+ bool have_matrix = false;
+ pure_expr *x = 0, *ret = 0;
+ va_start(ap, n);
+ for (size_t i = 0; i < n; i++)
+ xs[i] = va_arg(ap, pure_expr*);
+ va_end(ap);
for (size_t i = 0; i < n; i++) {
- pure_expr *x = va_arg(ap, pure_expr*);
+ x = xs[i];
switch (x->tag) {
case EXPR::DBL:
if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
case EXPR::INT:
case EXPR::BIGINT:
- if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
- xs[m++] = x; nrows++;
- k = 1;
+ if (k >= 0 && k != 1) goto err;
+ nrows++; k = 1;
break;
case EXPR::APP: {
pure_expr *u = x->data.x[0], *v = x->data.x[1];
@@ -1788,87 +1830,93 @@
if ((!rect || f->tag != rect->f) &&
(!polar || f->tag != polar->f) &&
f->tag != interp.symtab.pair_sym().f)
- pure_throw(bad_matrix_exception(x));
+ goto err;
u = u->data.x[1];
if (u->tag != EXPR::INT && u->tag != EXPR::BIGINT &&
u->tag != EXPR::DBL)
- pure_throw(bad_matrix_exception(x));
+ goto err;
if (v->tag != EXPR::INT && v->tag != EXPR::BIGINT &&
v->tag != EXPR::DBL)
- pure_throw(bad_matrix_exception(x));
+ goto err;
} else
- pure_throw(bad_matrix_exception(x));
- if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
+ goto err;
+ if (k >= 0 && k != 1) goto err;
target = EXPR::CMATRIX;
- xs[m++] = x; nrows++;
- k = 1;
+ nrows++; k = 1;
break;
}
case EXPR::MATRIX: {
gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
if (mp) {
if (k >= 0 && mp->size2 != (size_t)k)
- pure_throw(bad_matrix_exception(x));
- if (mp->size1 > 0) xs[m++] = x;
- nrows += mp->size1;
- k = mp->size2;
+ goto err;
+ nrows += mp->size1; k = mp->size2;
} else if (k>0)
- pure_throw(bad_matrix_exception(x));
+ goto err;
if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ have_matrix = true;
break;
}
case EXPR::CMATRIX: {
gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
if (mp) {
if (k >= 0 && mp->size2 != (size_t)k)
- pure_throw(bad_matrix_exception(x));
- if (mp->size1 > 0) xs[m++] = x;
- nrows += mp->size1;
- k = mp->size2;
+ goto err;
+ nrows += mp->size1; k = mp->size2;
} else if (k>0)
- pure_throw(bad_matrix_exception(x));
+ goto err;
target = EXPR::CMATRIX;
+ have_matrix = true;
break;
}
case EXPR::IMATRIX: {
gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
if (mp) {
if (k >= 0 && mp->size2 != (size_t)k)
- pure_throw(bad_matrix_exception(x));
- if (mp->size1 > 0) xs[m++] = x;
- nrows += mp->size1;
- k = mp->size2;
+ goto err;
+ nrows += mp->size1; k = mp->size2;
} else if (k>0)
- pure_throw(bad_matrix_exception(x));
+ goto err;
+ have_matrix = true;
break;
}
default:
- pure_throw(bad_matrix_exception(x));
+ goto err;
break;
}
}
- va_end(ap);
+ if (n == 1 && have_matrix) return xs[0];
if (k < 0) k = 0;
ncols = k;
if (nrows == 0 && ncols == 0) target = EXPR::MATRIX;
- pure_expr *ret = 0;
switch (target) {
case EXPR::MATRIX:
- ret = double_matrix_rows(nrows, ncols, m, xs);
+ ret = double_matrix_rows(nrows, ncols, n, xs);
break;
#if 0 // XXXTODO
case EXPR::CMATRIX:
- ret = complex_matrix_rows(nrows, ncols, m, xs);
+ ret = complex_matrix_rows(nrows, ncols, n, xs);
break;
#endif
case EXPR::IMATRIX:
- ret = int_matrix_rows(nrows, ncols, m, xs);
+ ret = int_matrix_rows(nrows, ncols, n, xs);
break;
default:
break;
}
if (!ret) pure_throw(not_implemented_exception());
return ret;
+ err:
+ /* This is called without a shadow stack frame, so we do our own cleanup
+ here to avoid having temporaries hanging around indefinitely. */
+ if (x) x->refc++;
+ for (size_t i = 0; i < n; i++)
+ pure_new_internal(xs[i]);
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ pure_unref_internal(x);
+ pure_throw(bad_matrix_exception(x));
+ return 0;
#else
pure_throw(not_implemented_exception());
return 0;
@@ -1880,21 +1928,25 @@
{
#ifdef HAVE_GSL
va_list ap;
- va_start(ap, n);
pure_expr **xs = (pure_expr**)alloca(n*sizeof(pure_expr*));
- size_t m = 0; int k = -1;
+ int k = -1;
size_t nrows = 0, ncols = 0;
int32_t target = EXPR::IMATRIX;
+ bool have_matrix = false;
+ pure_expr *x = 0, *ret = 0;
+ va_start(ap, n);
+ for (size_t i = 0; i < n; i++)
+ xs[i] = va_arg(ap, pure_expr*);
+ va_end(ap);
for (size_t i = 0; i < n; i++) {
- pure_expr *x = va_arg(ap, pure_expr*);
+ x = xs[i];
switch (x->tag) {
case EXPR::DBL:
if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
case EXPR::INT:
case EXPR::BIGINT:
- if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
- xs[m++] = x; ncols++;
- k = 1;
+ if (k >= 0 && k != 1) goto err;
+ ncols++; k = 1;
break;
case EXPR::APP: {
pure_expr *u = x->data.x[0], *v = x->data.x[1];
@@ -1906,87 +1958,93 @@
if ((!rect || f->tag != rect->f) &&
(!polar || f->tag != polar->f) &&
f->tag != interp.symtab.pair_sym().f)
- pure_throw(bad_matrix_exception(x));
+ goto err;
u = u->data.x[1];
if (u->tag != EXPR::INT && u->tag != EXPR::BIGINT &&
u->tag != EXPR::DBL)
- pure_throw(bad_matrix_exception(x));
+ goto err;
if (v->tag != EXPR::INT && v->tag != EXPR::BIGINT &&
v->tag != EXPR::DBL)
- pure_throw(bad_matrix_exception(x));
+ goto err;
} else
- pure_throw(bad_matrix_exception(x));
- if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
+ goto err;
+ if (k >= 0 && k != 1) goto err;
target = EXPR::CMATRIX;
- xs[m++] = x; ncols++;
- k = 1;
+ ncols++; k = 1;
break;
}
case EXPR::MATRIX: {
gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
if (mp) {
if (k >= 0 && mp->size1 != (size_t)k)
- pure_throw(bad_matrix_exception(x));
- if (mp->size2 > 0) xs[m++] = x;
- ncols += mp->size2;
- k = mp->size1;
+ goto err;
+ ncols += mp->size2; k = mp->size1;
} else if (k>0)
- pure_throw(bad_matrix_exception(x));
+ goto err;
if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ have_matrix = true;
break;
}
case EXPR::CMATRIX: {
gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
if (mp) {
if (k >= 0 && mp->size1 != (size_t)k)
- pure_throw(bad_matrix_exception(x));
- if (mp->size2 > 0) xs[m++] = x;
- ncols += mp->size2;
- k = mp->size1;
+ goto err;
+ ncols += mp->size2; k = mp->size1;
} else if (k>0)
- pure_throw(bad_matrix_exception(x));
+ goto err;
+ target = EXPR::CMATRIX;
+ have_matrix = true;
break;
}
case EXPR::IMATRIX: {
gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
if (mp) {
- if (k >= 0 && mp->size2 != (size_t)k)
- pure_throw(bad_matrix_exception(x));
- if (mp->size1 > 0) xs[m++] = x;
- ncols += mp->size2;
- k = mp->size2;
+ if (k >= 0 && mp->size1 != (size_t)k)
+ goto err;
+ ncols += mp->size2; k = mp->size1;
} else if (k>0)
- pure_throw(bad_matrix_exception(x));
- target = EXPR::CMATRIX;
+ goto err;
+ have_matrix = true;
break;
}
default:
- pure_throw(bad_matrix_exception(x));
+ goto err;
break;
}
}
- va_end(ap);
+ if (n == 1 && have_matrix) return xs[0];
if (k < 0) k = 0;
nrows = k;
if (nrows == 0 && ncols == 0) target = EXPR::MATRIX;
- pure_expr *ret = 0;
switch (target) {
case EXPR::MATRIX:
- ret = double_matrix_columns(nrows, ncols, m, xs);
+ ret = double_matrix_columns(nrows, ncols, n, xs);
break;
#if 0 // XXXTODO
case EXPR::CMATRIX:
- ret = complex_matrix_columns(nrows, ncols, m, xs);
+ ret = complex_matrix_columns(nrows, ncols, n, xs);
break;
#endif
case EXPR::IMATRIX:
- ret = int_matrix_columns(nrows, ncols, m, xs);
+ ret = int_matrix_columns(nrows, ncols, n, xs);
break;
default:
break;
}
if (!ret) pure_throw(not_implemented_exception());
return ret;
+ err:
+ /* This is called without a shadow stack frame, so we do our own cleanup
+ here to avoid having temporaries hanging around indefinitely. */
+ if (x) x->refc++;
+ for (size_t i = 0; i < n; i++)
+ pure_new_internal(xs[i]);
+ for (size_t i = 0; i < n; i++)
+ pure_free_internal(xs[i]);
+ pure_unref_internal(x);
+ pure_throw(bad_matrix_exception(x));
+ return 0;
#else
pure_throw(not_implemented_exception());
return 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-17 08:32:48 UTC (rev 776)
+++ pure/trunk/runtime.h 2008-09-17 15:05:12 UTC (rev 777)
@@ -148,9 +148,11 @@
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). */
+ other dimensions match. A null expression is returned in case of an error
+ (no matrix support, dimension mismatch, or invalid element type). Otherwise
+ a new matrix expression is returned. Temporary element expressions are
+ taken over by the callee and will be garbage-collected, but the elems
+ vectors are owned by the caller and won't be freed. */
pure_expr *pure_matrix_rowsl(uint32_t n, ...);
pure_expr *pure_matrix_rowsv(uint32_t n, pure_expr **elems);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-17 08:32:38
|
Revision: 776
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=776&view=rev
Author: agraef
Date: 2008-09-17 08:32:48 +0000 (Wed, 17 Sep 2008)
Log Message:
-----------
Various fixes, partial implementation of matrix constructors.
Modified Paths:
--------------
pure/trunk/interpreter.cc
pure/trunk/printer.cc
pure/trunk/runtime.cc
pure/trunk/symtable.cc
pure/trunk/symtable.hh
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-16 18:46:30 UTC (rev 775)
+++ pure/trunk/interpreter.cc 2008-09-17 08:32:48 UTC (rev 776)
@@ -913,13 +913,15 @@
if (x->data.mat.p) {
gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
exprll *xs = new exprll;
+ symbol *rect = symtab.complex_rect_sym();
+ expr f = rect?rect->x:symtab.pair_sym().x;
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));
+ ys.push_back(expr(f, u, v));
}
}
return expr(EXPR::MATRIX, m);
@@ -5031,7 +5033,10 @@
us[i] =
act_env().CreateCall(module->getFunction("pure_matrix_columns"), vs);
}
- return act_env().CreateCall(module->getFunction("pure_matrix_rows"), us);
+ if (n == 1)
+ return us[1];
+ else
+ return act_env().CreateCall(module->getFunction("pure_matrix_rows"), us);
}
// application:
case EXPR::APP:
Modified: pure/trunk/printer.cc
===================================================================
--- pure/trunk/printer.cc 2008-09-16 18:46:30 UTC (rev 775)
+++ pure/trunk/printer.cc 2008-09-17 08:32:48 UTC (rev 776)
@@ -713,6 +713,26 @@
return false;
}
+static inline ostream& print_double(ostream& os, double d)
+{
+ char buf[64];
+ if (is_inf(d))
+ if (d > 0)
+ strcpy(buf, "inf");
+ else
+ strcpy(buf, "-inf");
+ else if (is_nan(d))
+ strcpy(buf, "nan");
+ else
+ my_formatd(buf, "%0.15g", d);
+ // make sure that the output conforms to Pure syntax
+ os << buf;
+ if (strchr("0123456789", buf[buf[0]=='-'?1:0]) &&
+ !strchr(buf, '.') && !strchr(buf, 'e') && !strchr(buf, 'E'))
+ os << ".0";
+ return os;
+}
+
ostream& operator << (ostream& os, const pure_expr *x)
{
char test;
@@ -720,7 +740,6 @@
interpreter::stackdir*(&test - interpreter::baseptr) >=
interpreter::stackmax)
throw err("stack overflow in printer");
- char buf[64];
assert(x);
if (pstr(os, (pure_expr*)x)) return os;
//os << "{" << x->refc << "}";
@@ -732,24 +751,8 @@
os << s << "L"; free(s);
return os;
}
- case EXPR::DBL: {
- double d = x->data.d;
- if (is_inf(d))
- if (d > 0)
- strcpy(buf, "inf");
- else
- strcpy(buf, "-inf");
- else if (is_nan(d))
- strcpy(buf, "nan");
- else
- my_formatd(buf, "%0.15g", d);
- // make sure that the output conforms to Pure syntax
- os << buf;
- if (strchr("0123456789", buf[buf[0]=='-'?1:0]) &&
- !strchr(buf, '.') && !strchr(buf, 'e') && !strchr(buf, 'E'))
- os << ".0";
- return os;
- }
+ case EXPR::DBL:
+ return print_double(os, x->data.d);
case EXPR::STR: {
char *s = printstr(x->data.s);
os << '"' << s << '"';
@@ -758,6 +761,7 @@
}
case EXPR::PTR:
return os << "#<pointer " << x->data.p << ">";
+#ifdef HAVE_GSL
/* 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. */
@@ -769,7 +773,7 @@
if (i > 0) os << ";";
for (size_t j = 0; j < m->size2; j++) {
if (j > 0) os << ",";
- os << m->data[i * m->tda + j];
+ print_double(os, m->data[i * m->tda + j]);
}
}
}
@@ -791,21 +795,44 @@
os << "{";
if (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];
+ /* GSL represents complex matrices using pairs of double values, while
+ Pure provides its own complex type in math.pure. If math.pure has
+ been loaded, then the '+:' operator is defined and we use this
+ representation. Otherwise, we print complex values as pairs of real
+ and imaginary part. */
+ symbol *rect = interpreter::g_interp->symtab.complex_rect_sym();
+ if (rect)
+ 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 << ",";
+ print_double(os, m->data[2*(i * m->tda + j)]);
+ os << rect->s;
+ print_double(os, m->data[2*(i * m->tda + j) + 1]);
+ }
}
- }
+ else
+ 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 << "(";
+ print_double(os, m->data[2*(i * m->tda + j)]);
+ os << ",";
+ print_double(os, m->data[2*(i * m->tda + j) + 1]);
+ os << ")";
+ }
+ }
}
return os << "}";
+#else
+ case EXPR::MATRIX:
+ return os << "#<matrix " << x->data.mat.p << ">";
+ case EXPR::IMATRIX:
+ return os << "#<imatrix " << x->data.mat.p << ">";
+ case EXPR::CMATRIX:
+ return os << "#<cmatrix " << x->data.mat.p << ">";
+#endif
case EXPR::APP: {
list<const pure_expr*> xs;
prec_t p;
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-16 18:46:30 UTC (rev 775)
+++ pure/trunk/runtime.cc 2008-09-17 08:32:48 UTC (rev 776)
@@ -301,7 +301,7 @@
break;
}
#endif
- if (owner) free(x->data.mat.refc);
+ if (owner) delete x->data.mat.refc;
}
#if 1
@@ -589,7 +589,15 @@
pure_expr *pure_double_matrix(void *p)
{
#ifdef HAVE_GSL
- return 0; // XXXTODO
+ gsl_matrix *m = (gsl_matrix*)p;
+ if (!m || !m->owner) return 0;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::MATRIX;
+ x->data.mat.p = p;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
#else
return 0;
#endif
@@ -599,7 +607,15 @@
pure_expr *pure_complex_matrix(void *p)
{
#ifdef HAVE_GSL
- return 0; // XXXTODO
+ gsl_matrix_complex *m = (gsl_matrix_complex*)p;
+ if (!m || !m->owner) return 0;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = p;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
#else
return 0;
#endif
@@ -609,7 +625,15 @@
pure_expr *pure_int_matrix(void *p)
{
#ifdef HAVE_GSL
- return 0; // XXXTODO
+ gsl_matrix_int *m = (gsl_matrix_int*)p;
+ if (!m || !m->owner) return 0;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = p;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
#else
return 0;
#endif
@@ -619,7 +643,11 @@
pure_expr *pure_double_matrix_dup(const void *p)
{
#ifdef HAVE_GSL
- return 0; // XXXTODO
+ gsl_matrix *m1 = (gsl_matrix*)p;
+ if (!m1) return 0;
+ gsl_matrix *m2 = gsl_matrix_alloc(m1->size1, m1->size2);
+ gsl_matrix_memcpy(m2, m1);
+ return pure_double_matrix(m2);
#else
return 0;
#endif
@@ -629,7 +657,11 @@
pure_expr *pure_complex_matrix_dup(const void *p)
{
#ifdef HAVE_GSL
- return 0; // XXXTODO
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)p;
+ if (!m1) return 0;
+ gsl_matrix_complex *m2 = gsl_matrix_complex_alloc(m1->size1, m1->size2);
+ gsl_matrix_complex_memcpy(m2, m1);
+ return pure_complex_matrix(m2);
#else
return 0;
#endif
@@ -639,12 +671,194 @@
pure_expr *pure_int_matrix_dup(const void *p)
{
#ifdef HAVE_GSL
- return 0; // XXXTODO
+ gsl_matrix_int *m1 = (gsl_matrix_int*)p;
+ if (!m1) return 0;
+ gsl_matrix_int *m2 = gsl_matrix_int_alloc(m1->size1, m1->size2);
+ gsl_matrix_int_memcpy(m2, m1);
+ return pure_int_matrix(m2);
#else
return 0;
#endif
}
+#ifdef HAVE_GSL
+static pure_expr*
+double_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
+{
+ gsl_matrix *mat = gsl_matrix_alloc(nrows, ncols);
+ if (!mat) return 0; // XXXTODO: empty matrices
+ double *data = mat->data;
+ size_t tda = mat->tda;
+ for (size_t count = 0, i = 0; count < n; count++) {
+ pure_expr *x = xs[count];
+ switch (x->tag) {
+ case EXPR::INT:
+ data[i++*tda] = (double)x->data.i;
+ break;
+ case EXPR::BIGINT:
+ data[i++*tda] = mpz_get_d(x->data.z);
+ break;
+ case EXPR::DBL:
+ data[i++*tda] = x->data.d;
+ break;
+ case EXPR::MATRIX: {
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ memcpy(data+i*tda, mat1->data+j*mat1->tda, ncols*sizeof(double));
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ for (size_t k = 0; k < mat1->size2; k++)
+ data[i*tda+k] = (double)mat1->data[j*mat1->tda+k];
+ break;
+ }
+ default:
+ assert(0 && "bad matrix element");
+ break;
+ }
+ pure_freenew(x);
+ }
+ return pure_double_matrix(mat);
+}
+
+static pure_expr*
+double_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
+{
+ gsl_matrix *mat = gsl_matrix_alloc(nrows, ncols);
+ if (!mat) return 0; // XXXTODO: empty matrices
+ double *data = mat->data;
+ size_t tda = mat->tda;
+ for (size_t count = 0, i = 0; count < n; count++) {
+ pure_expr *x = xs[count];
+ switch (x->tag) {
+ case EXPR::INT:
+ data[i++] = (double)x->data.i;
+ break;
+ case EXPR::BIGINT:
+ data[i++] = mpz_get_d(x->data.z);
+ break;
+ case EXPR::DBL:
+ data[i++] = x->data.d;
+ break;
+ case EXPR::MATRIX: {
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ memcpy(data+j*tda+i, mat1->data+j*mat1->tda, ncols*sizeof(double));
+ i += mat1->size2;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; j++)
+ for (size_t k = 0; k < mat1->size2; k++)
+ data[j*tda+k+i] = (double)mat1->data[j*mat1->tda+k];
+ i += mat1->size2;
+ break;
+ }
+ default:
+ assert(0 && "bad matrix element");
+ break;
+ }
+ pure_freenew(x);
+ }
+ return pure_double_matrix(mat);
+}
+
+static pure_expr*
+int_matrix_rows(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
+{
+ gsl_matrix_int *mat = gsl_matrix_int_alloc(nrows, ncols);
+ if (!mat) return 0; // XXXTODO: empty matrices
+ int *data = mat->data;
+ size_t tda = mat->tda;
+ for (size_t count = 0, i = 0; count < n; count++) {
+ pure_expr *x = xs[count];
+ switch (x->tag) {
+ case EXPR::INT:
+ data[i++*tda] = x->data.i;
+ break;
+ case EXPR::BIGINT:
+ data[i++*tda] = pure_get_int(x);
+ break;
+ case EXPR::DBL:
+ data[i++*tda] = (int)x->data.d;
+ break;
+ case EXPR::MATRIX: {
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ for (size_t k = 0; k < mat1->size2; k++)
+ data[i*tda+k] = (int)mat1->data[j*mat1->tda+k];
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ memcpy(data+i*tda, mat1->data+j*mat1->tda, ncols*sizeof(int));
+ break;
+ }
+ default:
+ assert(0 && "bad matrix element");
+ break;
+ }
+ pure_freenew(x);
+ }
+ return pure_int_matrix(mat);
+}
+
+static pure_expr*
+int_matrix_columns(size_t nrows, size_t ncols, size_t n, pure_expr **xs)
+{
+ gsl_matrix_int *mat = gsl_matrix_int_alloc(nrows, ncols);
+ if (!mat) return 0; // XXXTODO: empty matrices
+ int *data = mat->data;
+ size_t tda = mat->tda;
+ for (size_t count = 0, i = 0; count < n; count++) {
+ pure_expr *x = xs[count];
+ switch (x->tag) {
+ case EXPR::INT:
+ data[i++] = x->data.i;
+ break;
+ case EXPR::BIGINT:
+ data[i++] = pure_get_int(x);
+ break;
+ case EXPR::DBL:
+ data[i++] = (int)x->data.d;
+ break;
+ case EXPR::MATRIX: {
+ gsl_matrix_int *mat1 = (gsl_matrix_int*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; j++)
+ for (size_t k = 0; k < mat1->size2; k++)
+ data[j*tda+k+i] = (int)mat1->data[j*mat1->tda+k];
+ i += mat1->size2;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix *mat1 = (gsl_matrix*)x->data.mat.p;
+ if (mat1)
+ for (size_t j = 0; j < mat1->size1; i++, j++)
+ memcpy(data+j*tda+i, mat1->data+j*mat1->tda, ncols*sizeof(int));
+ i += mat1->size2;
+ break;
+ }
+ default:
+ assert(0 && "bad matrix element");
+ break;
+ }
+ pure_freenew(x);
+ }
+ return pure_int_matrix(mat);
+}
+#endif
+
extern "C"
pure_expr *pure_matrix_rowsl(uint32_t n, ...)
{
@@ -1547,9 +1761,114 @@
pure_expr *pure_matrix_rows(uint32_t n, ...)
{
#ifdef HAVE_GSL
- // XXXTODO
- pure_throw(not_implemented_exception());
- return 0;
+ va_list ap;
+ va_start(ap, n);
+ pure_expr **xs = (pure_expr**)alloca(n*sizeof(pure_expr*));
+ size_t m = 0; int k = -1;
+ size_t nrows = 0, ncols = 0;
+ int32_t target = EXPR::IMATRIX;
+ for (size_t i = 0; i < n; i++) {
+ pure_expr *x = va_arg(ap, pure_expr*);
+ switch (x->tag) {
+ case EXPR::DBL:
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ case EXPR::INT:
+ case EXPR::BIGINT:
+ if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
+ xs[m++] = x; nrows++;
+ k = 1;
+ break;
+ case EXPR::APP: {
+ pure_expr *u = x->data.x[0], *v = x->data.x[1];
+ if (u->tag == EXPR::APP) {
+ interpreter& interp = *interpreter::g_interp;
+ pure_expr *f = u->data.x[0];
+ symbol *rect = interp.symtab.complex_rect_sym(),
+ *polar = interp.symtab.complex_polar_sym();
+ if ((!rect || f->tag != rect->f) &&
+ (!polar || f->tag != polar->f) &&
+ f->tag != interp.symtab.pair_sym().f)
+ pure_throw(bad_matrix_exception(x));
+ u = u->data.x[1];
+ if (u->tag != EXPR::INT && u->tag != EXPR::BIGINT &&
+ u->tag != EXPR::DBL)
+ pure_throw(bad_matrix_exception(x));
+ if (v->tag != EXPR::INT && v->tag != EXPR::BIGINT &&
+ v->tag != EXPR::DBL)
+ pure_throw(bad_matrix_exception(x));
+ } else
+ pure_throw(bad_matrix_exception(x));
+ if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
+ target = EXPR::CMATRIX;
+ xs[m++] = x; nrows++;
+ k = 1;
+ break;
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ pure_throw(bad_matrix_exception(x));
+ if (mp->size1 > 0) xs[m++] = x;
+ nrows += mp->size1;
+ k = mp->size2;
+ } else if (k>0)
+ pure_throw(bad_matrix_exception(x));
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ pure_throw(bad_matrix_exception(x));
+ if (mp->size1 > 0) xs[m++] = x;
+ nrows += mp->size1;
+ k = mp->size2;
+ } else if (k>0)
+ pure_throw(bad_matrix_exception(x));
+ target = EXPR::CMATRIX;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ pure_throw(bad_matrix_exception(x));
+ if (mp->size1 > 0) xs[m++] = x;
+ nrows += mp->size1;
+ k = mp->size2;
+ } else if (k>0)
+ pure_throw(bad_matrix_exception(x));
+ break;
+ }
+ default:
+ pure_throw(bad_matrix_exception(x));
+ break;
+ }
+ }
+ va_end(ap);
+ if (k < 0) k = 0;
+ ncols = k;
+ if (nrows == 0 && ncols == 0) target = EXPR::MATRIX;
+ pure_expr *ret = 0;
+ switch (target) {
+ case EXPR::MATRIX:
+ ret = double_matrix_rows(nrows, ncols, m, xs);
+ break;
+#if 0 // XXXTODO
+ case EXPR::CMATRIX:
+ ret = complex_matrix_rows(nrows, ncols, m, xs);
+ break;
+#endif
+ case EXPR::IMATRIX:
+ ret = int_matrix_rows(nrows, ncols, m, xs);
+ break;
+ default:
+ break;
+ }
+ if (!ret) pure_throw(not_implemented_exception());
+ return ret;
#else
pure_throw(not_implemented_exception());
return 0;
@@ -1560,9 +1879,114 @@
pure_expr *pure_matrix_columns(uint32_t n, ...)
{
#ifdef HAVE_GSL
- // XXXTODO
- pure_throw(not_implemented_exception());
- return 0;
+ va_list ap;
+ va_start(ap, n);
+ pure_expr **xs = (pure_expr**)alloca(n*sizeof(pure_expr*));
+ size_t m = 0; int k = -1;
+ size_t nrows = 0, ncols = 0;
+ int32_t target = EXPR::IMATRIX;
+ for (size_t i = 0; i < n; i++) {
+ pure_expr *x = va_arg(ap, pure_expr*);
+ switch (x->tag) {
+ case EXPR::DBL:
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ case EXPR::INT:
+ case EXPR::BIGINT:
+ if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
+ xs[m++] = x; ncols++;
+ k = 1;
+ break;
+ case EXPR::APP: {
+ pure_expr *u = x->data.x[0], *v = x->data.x[1];
+ if (u->tag == EXPR::APP) {
+ interpreter& interp = *interpreter::g_interp;
+ pure_expr *f = u->data.x[0];
+ symbol *rect = interp.symtab.complex_rect_sym(),
+ *polar = interp.symtab.complex_polar_sym();
+ if ((!rect || f->tag != rect->f) &&
+ (!polar || f->tag != polar->f) &&
+ f->tag != interp.symtab.pair_sym().f)
+ pure_throw(bad_matrix_exception(x));
+ u = u->data.x[1];
+ if (u->tag != EXPR::INT && u->tag != EXPR::BIGINT &&
+ u->tag != EXPR::DBL)
+ pure_throw(bad_matrix_exception(x));
+ if (v->tag != EXPR::INT && v->tag != EXPR::BIGINT &&
+ v->tag != EXPR::DBL)
+ pure_throw(bad_matrix_exception(x));
+ } else
+ pure_throw(bad_matrix_exception(x));
+ if (k >= 0 && k != 1) pure_throw(bad_matrix_exception(x));
+ target = EXPR::CMATRIX;
+ xs[m++] = x; ncols++;
+ k = 1;
+ break;
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ pure_throw(bad_matrix_exception(x));
+ if (mp->size2 > 0) xs[m++] = x;
+ ncols += mp->size2;
+ k = mp->size1;
+ } else if (k>0)
+ pure_throw(bad_matrix_exception(x));
+ if (target == EXPR::IMATRIX) target = EXPR::MATRIX;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ pure_throw(bad_matrix_exception(x));
+ if (mp->size2 > 0) xs[m++] = x;
+ ncols += mp->size2;
+ k = mp->size1;
+ } else if (k>0)
+ pure_throw(bad_matrix_exception(x));
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
+ if (mp) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ pure_throw(bad_matrix_exception(x));
+ if (mp->size1 > 0) xs[m++] = x;
+ ncols += mp->size2;
+ k = mp->size2;
+ } else if (k>0)
+ pure_throw(bad_matrix_exception(x));
+ target = EXPR::CMATRIX;
+ break;
+ }
+ default:
+ pure_throw(bad_matrix_exception(x));
+ break;
+ }
+ }
+ va_end(ap);
+ if (k < 0) k = 0;
+ nrows = k;
+ if (nrows == 0 && ncols == 0) target = EXPR::MATRIX;
+ pure_expr *ret = 0;
+ switch (target) {
+ case EXPR::MATRIX:
+ ret = double_matrix_columns(nrows, ncols, m, xs);
+ break;
+#if 0 // XXXTODO
+ case EXPR::CMATRIX:
+ ret = complex_matrix_columns(nrows, ncols, m, xs);
+ break;
+#endif
+ case EXPR::IMATRIX:
+ ret = int_matrix_columns(nrows, ncols, m, xs);
+ break;
+ default:
+ break;
+ }
+ if (!ret) pure_throw(not_implemented_exception());
+ return ret;
#else
pure_throw(not_implemented_exception());
return 0;
Modified: pure/trunk/symtable.cc
===================================================================
--- pure/trunk/symtable.cc 2008-09-16 18:46:30 UTC (rev 775)
+++ pure/trunk/symtable.cc 2008-09-17 08:32:48 UTC (rev 776)
@@ -368,24 +368,6 @@
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 18:46:30 UTC (rev 775)
+++ pure/trunk/symtable.hh 2008-09-17 08:32:48 UTC (rev 776)
@@ -90,8 +90,6 @@
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"); }
@@ -101,6 +99,9 @@
symbol& not_implemented_sym() { return sym("not_implemented"); }
symbol& bad_matrix_sym() { return sym("bad_matrix_value"); }
symbol& amp_sym();
+ // these may be undefined
+ symbol* complex_rect_sym() { return lookup("+:"); }
+ symbol* complex_polar_sym() { return lookup("<:"); }
};
#endif // ! SYMTABLE_HH
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-16 18:46:20
|
Revision: 775
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=775&view=rev
Author: agraef
Date: 2008-09-16 18:46:30 +0000 (Tue, 16 Sep 2008)
Log Message:
-----------
Add support for cmatrix and imatrix tags.
Modified Paths:
--------------
pure/trunk/lexer.ll
pure/trunk/printer.cc
Modified: pure/trunk/lexer.ll
===================================================================
--- pure/trunk/lexer.ll 2008-09-16 18:36:59 UTC (rev 774)
+++ pure/trunk/lexer.ll 2008-09-16 18:46:30 UTC (rev 775)
@@ -255,6 +255,8 @@
strtag ::{blank}*string
ptrtag ::{blank}*pointer
mattag ::{blank}*matrix
+cmattag ::{blank}*cmatrix
+imattag ::{blank}*imatrix
%x comment xdecl xdecl_comment xusing xusing_comment
@@ -1187,6 +1189,8 @@
{strtag}/[^a-zA-Z_0-9] yylval->ival = EXPR::STR; return token::TAG;
{ptrtag}/[^a-zA-Z_0-9] yylval->ival = EXPR::PTR; return token::TAG;
{mattag}/[^a-zA-Z_0-9] yylval->ival = EXPR::MATRIX; return token::TAG;
+{cmattag}/[^a-zA-Z_0-9] yylval->ival = EXPR::CMATRIX; return token::TAG;
+{imattag}/[^a-zA-Z_0-9] yylval->ival = EXPR::IMATRIX; return token::TAG;
extern BEGIN(xdecl); return token::EXTERN;
infix yylval->fix = infix; return token::FIX;
infixl yylval->fix = infixl; return token::FIX;
Modified: pure/trunk/printer.cc
===================================================================
--- pure/trunk/printer.cc 2008-09-16 18:36:59 UTC (rev 774)
+++ pure/trunk/printer.cc 2008-09-16 18:46:30 UTC (rev 775)
@@ -169,6 +169,10 @@
return os << "::string";
case EXPR::MATRIX:
return os << "::matrix";
+ case EXPR::CMATRIX:
+ return os << "::cmatrix";
+ case EXPR::IMATRIX:
+ return os << "::imatrix";
default:
return os;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
