Thread: [pure-lang-svn] SF.net SVN: pure-lang:[812] pure/trunk (Page 6)
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
pure-lang-svn
|
From: <ag...@us...> - 2008-09-20 20:33:23
|
Revision: 812
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=812&view=rev
Author: agraef
Date: 2008-09-20 20:33:16 +0000 (Sat, 20 Sep 2008)
Log Message:
-----------
Update ChangeLog and NEWS.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/NEWS
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-20 18:50:23 UTC (rev 811)
+++ pure/trunk/ChangeLog 2008-09-20 20:33:16 UTC (rev 812)
@@ -1,69 +1,15 @@
-2008-09-19 Albert Graef <Dr....@t-...>
-
- * runtime.cc, interpreter.cc et al: Support for symbolic matrices.
+2008-09-20 Albert Graef <Dr....@t-...>
- These work like the numeric GSL matrices, but can contain an
- arbitrary mixture of Pure values, and also work if GSL support
- isn't available. Numeric matrices will now only be created for
- matrices consisting of machine int, double and complex values, if
- all values are of the same type (and GSL support is available). In
- all other cases a symbolic matrix is created (this also includes
- the case of bigints, which are now considered as symbolic values
- in matrix creation).
+ * Implemented basic GSL matrix support, including support for
+ symbolic matrices (which is independent from GSL, so these will
+ also work when building the interpreter without GSL) and matrix
+ comprehensions. This required many additions and changes to the
+ parser, interpreter, compiler, runtime and the prelude; details
+ can be found in the svn log (see r759 and r769ff.).
- There's now only a single 'matrix' type tag which matches all
- kinds of symbolic and numeric matrices. The matrixp family of
- predicates in primitives.pure can be used to distinguish the
- different kinds of matrices. Predicates to check for row and
- column vectors (which are simply matrices with one row or column)
- are also provided, as are some routines for transposing and
- converting matrix values, see primitives.pure for details.
+ Preliminary documentation is in the NEWS file for now (the manual
+ still needs to be updated).
-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. 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, 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
- values, such as {1+:2,3<:4} (the latter requires math.pure to be
- loaded). The expression printer uses the former notation, unless
- 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 individual matrix elements, but of course it is possible
- to override the default print representation of matrix values as a
- whole.
-
2008-09-15 Albert Graef <Dr....@t-...>
* configure.ac: Bump version number.
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-20 18:50:23 UTC (rev 811)
+++ pure/trunk/NEWS 2008-09-20 20:33:16 UTC (rev 812)
@@ -1,8 +1,106 @@
** Pure 0.7 (in progress)
-GSL matrix support. (Work in progress.)
+Basic GSL (GNU Scientific Library) matrix support has been implemented,
+including matrix comprehensions. Here's some preliminary documentation (at the
+time of this writing, the manual still needs to be updated). For more
+information on GSL please refer to http://www.gnu.org/software/gsl.
+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. Also, indices are zero-based (rather than
+one-based) to be consistent with Pure's list and tuple structures.
+
+Here are some simple examples of matrices:
+
+- {1,2,3} a row vector consisting of machine ints
+- {1.0;2.0;3.0} a column vector of double values (';' separates rows)
+- {1,2;3,4} a 2x2 int matrix
+- {1L,y+1;foo,bar} a "symbolic" matrix (see explanation below)
+
+Note that the {...} syntax can be used only on the right-hand side of a
+definition to construct new matrices. Matrix patterns are not supported.
+However, the new 'matrix' type tag can be used in patterns to match matrix
+values, e.g.:
+
+foo a::matrix = {a!(i,j)+1 | i=0..n-1; j=0..m-1} when n,m = dim a end;
+
+(Additional predicates to check for the different subtypes of matrices are
+available in the prelude.)
+
+As mentioned above, there's also a special built-in type of "symbolic"
+matrices. These work like the numeric matrices, but can contain an arbitrary
+mixture of Pure values, and also work if GSL support isn't available.
+
+GSL matrices are always homogeneous, i.e., they only contain values from one
+numeric type, which in the Pure implementation can be machine ints, double and
+complex double values. (The latter are represented using Pure's infix notation
+for complex values defined in math.pure; this requires math.pure to be
+loaded. The same notation is also used to print complex matrices. Note that,
+for performance reasons, the expression printer doesn't use the __show__
+function to print individual matrix elements. It is possible to override the
+default print representation of matrix values as a whole, though.)
+
+If a matrix contains values of different types, or Pure values which cannot be
+stored in a GSL matrix, then a symbolic matrix is created instead (this also
+includes the case of bigints, which are considered as symbolic values as far
+as matrix construction is concerned). If the interpreter was built without GSL
+support then symbolic matrices are the only kind of matrices supported by the
+interpreter.
+
+Pure also provides so-called matrix comprehensions as a convenient means to
+create matrices from a template expression (which can denote either a scalar
+or a submatrix), drawing values from lists and (optionally) filtering the
+elements with predicates. These work pretty much like list comprehensions, but
+return matrices instead of lists. Generator clauses in matrix comprehensions
+alternate between row and column generation so that most common mathematical
+abbreviations carry over quite easily. E.g., here's a simple example which
+illustrates how you can define a function which returns a square identity
+matrix of a given dimension:
+
+> eye n = {i==j|i=1..n;j=1..n};
+> eye 3;
+{1,0,0;0,1,0;0,0,1}
+
+The prelude provides some additional basic matrix operations like determining
+the size and dimensions of a matrix, indexing and slicing, transposition,
+type-checking and various conversions between the different kinds of
+matrices. To these ends, various basic operations to deal with matrix objects
+have been added to the runtime, including some public API operations to create
+and inspect Pure matrix values from external C modules.
+
+Here is a brief synopsis of some important operations which are already
+implemented in the prelude (you can find the definitions of these and the
+other matrix operations in primitives.pure):
+
+- #x total number of elements
+- dim x number of rows and columns (as a pair)
+- x' transposed matrix
+- x!i ith matrix element in row-major order (zero-based)
+- x!(i,j) jth element in ith row of the matrix (zero-based)
+- x!!is, x!!(is,js) slicing (is and js are lists of machine ints)
+- x==y, x!=y matrix comparisons
+
+Other user-visible changes prompted by the introduction of the matrix syntax
+are listed below:
+
+- Changes in the comprehension syntax: '|' is now used to separate the
+template expression and the generator and filter clauses. This change was
+necessary to accommodate the matrix syntax which uses ';' to separate
+different rows in a matrix. For consistency, this applies to both list and
+matrix comprehensions. The old [x;...] syntax is still supported in list
+comprehensions, but is flagged with a "deprecated" warning by the compiler.
+
+- Arithmetic sequences with a stepsize different from 1 are now written
+x:y..z instead of x,y..z. This makes it possible to give the .. operator a
+lower precedence than the ',' operator, which makes writing matrix slices like
+x!!(i..j,k..l) much more convenient.
+
+Stuff that's still missing:
+
+- Marshalling of GSL matrices in the C interface, so that it becomes possible
+to access all the advanced functionality available in GSL.
+
** Pure 0.6 2008-09-12
New stuff in this release (please see the ChangeLog and the manual for
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-21 05:24:41
|
Revision: 814
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=814&view=rev
Author: agraef
Date: 2008-09-21 05:24:34 +0000 (Sun, 21 Sep 2008)
Log Message:
-----------
Marshalling of matrix values in the C interface.
Modified Paths:
--------------
pure/trunk/interpreter.cc
pure/trunk/interpreter.hh
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-20 20:53:03 UTC (rev 813)
+++ pure/trunk/interpreter.cc 2008-09-21 05:24:34 UTC (rev 814)
@@ -126,6 +126,98 @@
// Char pointer type.
CharPtrTy = PointerType::get(Type::Int8Ty, 0);
+ // int and double pointers.
+ IntPtrTy = PointerType::get(Type::Int32Ty, 0);
+ DoublePtrTy = PointerType::get(Type::DoubleTy, 0);
+
+ // Complex numbers (complex double).
+ {
+ std::vector<const Type*> elts;
+ elts.push_back(Type::DoubleTy);
+ elts.push_back(Type::DoubleTy);
+ ComplexTy = StructType::get(elts);
+ ComplexPtrTy = PointerType::get(ComplexTy, 0);
+ }
+
+ // GSL matrix types. These are used to marshall GSL matrices in the C
+ // interface.
+ {
+ std::vector<const Type*> elts;
+ if (sizeof(size_t) == 4) {
+ elts.push_back(Type::Int32Ty); // size1
+ elts.push_back(Type::Int32Ty); // size2
+ elts.push_back(Type::Int32Ty); // tda
+ } else {
+ assert(sizeof(size_t) == 8);
+ elts.push_back(Type::Int64Ty); // size1
+ elts.push_back(Type::Int64Ty); // size2
+ elts.push_back(Type::Int64Ty); // tda
+ }
+ elts.push_back(VoidPtrTy); // data
+ elts.push_back(VoidPtrTy); // block
+ elts.push_back(Type::Int32Ty); // owner
+ GSLMatrixTy = StructType::get(elts);
+ module->addTypeName("struct.gsl_matrix", GSLMatrixTy);
+ GSLMatrixPtrTy = PointerType::get(GSLMatrixTy, 0);
+ }
+ {
+ std::vector<const Type*> elts;
+ if (sizeof(size_t) == 4) {
+ elts.push_back(Type::Int32Ty); // size1
+ elts.push_back(Type::Int32Ty); // size2
+ elts.push_back(Type::Int32Ty); // tda
+ } else {
+ assert(sizeof(size_t) == 8);
+ elts.push_back(Type::Int64Ty); // size1
+ elts.push_back(Type::Int64Ty); // size2
+ elts.push_back(Type::Int64Ty); // tda
+ }
+ elts.push_back(DoublePtrTy); // data
+ elts.push_back(VoidPtrTy); // block
+ elts.push_back(Type::Int32Ty); // owner
+ GSLDoubleMatrixTy = StructType::get(elts);
+ module->addTypeName("struct.gsl_matrix_double", GSLDoubleMatrixTy);
+ GSLDoubleMatrixPtrTy = PointerType::get(GSLDoubleMatrixTy, 0);
+ }
+ {
+ std::vector<const Type*> elts;
+ if (sizeof(size_t) == 4) {
+ elts.push_back(Type::Int32Ty); // size1
+ elts.push_back(Type::Int32Ty); // size2
+ elts.push_back(Type::Int32Ty); // tda
+ } else {
+ assert(sizeof(size_t) == 8);
+ elts.push_back(Type::Int64Ty); // size1
+ elts.push_back(Type::Int64Ty); // size2
+ elts.push_back(Type::Int64Ty); // tda
+ }
+ elts.push_back(ComplexPtrTy); // data
+ elts.push_back(VoidPtrTy); // block
+ elts.push_back(Type::Int32Ty); // owner
+ GSLComplexMatrixTy = StructType::get(elts);
+ module->addTypeName("struct.gsl_matrix_complex", GSLComplexMatrixTy);
+ GSLComplexMatrixPtrTy = PointerType::get(GSLComplexMatrixTy, 0);
+ }
+ {
+ std::vector<const Type*> elts;
+ if (sizeof(size_t) == 4) {
+ elts.push_back(Type::Int32Ty); // size1
+ elts.push_back(Type::Int32Ty); // size2
+ elts.push_back(Type::Int32Ty); // tda
+ } else {
+ assert(sizeof(size_t) == 8);
+ elts.push_back(Type::Int64Ty); // size1
+ elts.push_back(Type::Int64Ty); // size2
+ elts.push_back(Type::Int64Ty); // tda
+ }
+ elts.push_back(IntPtrTy); // data
+ elts.push_back(VoidPtrTy); // block
+ elts.push_back(Type::Int32Ty); // owner
+ GSLIntMatrixTy = StructType::get(elts);
+ module->addTypeName("struct.gsl_matrix_int", GSLIntMatrixTy);
+ GSLIntMatrixPtrTy = PointerType::get(GSLIntMatrixTy, 0);
+ }
+
// Create the expr struct type.
/* NOTE: This is in fact just a part of the prologue of the expression data
@@ -248,9 +340,17 @@
"pure_apply", "expr*", 2, "expr*", "expr*");
declare_extern((void*)pure_matrix_rows,
- "pure_matrix_rows", "expr*", -1, "int");
+ "pure_matrix_rows", "expr*", -1, "int");
declare_extern((void*)pure_matrix_columns,
"pure_matrix_columns", "expr*", -1, "int");
+ declare_extern((void*)pure_symbolic_matrix,
+ "pure_symbolic_matrix", "expr*", 1, "void*");
+ declare_extern((void*)pure_double_matrix,
+ "pure_double_matrix", "expr*", 1, "void*");
+ declare_extern((void*)pure_complex_matrix,
+ "pure_complex_matrix", "expr*", 1, "void*");
+ declare_extern((void*)pure_int_matrix,
+ "pure_int_matrix", "expr*", 1, "void*");
declare_extern((void*)pure_listl,
"pure_listl", "expr*", -1, "int");
@@ -273,6 +373,8 @@
"pure_get_long", "long", 1, "expr*");
declare_extern((void*)pure_get_int,
"pure_get_int", "int", 1, "expr*");
+ declare_extern((void*)pure_get_matrix,
+ "pure_get_matrix", "void*", 1, "expr*");
declare_extern((void*)pure_catch,
"pure_catch", "expr*", 2, "expr*", "expr*");
@@ -3562,6 +3664,10 @@
return Type::FloatTy;
else if (name == "double")
return Type::DoubleTy;
+#if 0 // no marshalling available yet, does LLVM support these?
+ else if (name == "complex")
+ return ComplexTy;
+#endif
else if (name == "char*")
return CharPtrTy;
else if (name == "short*")
@@ -3572,10 +3678,22 @@
return PointerType::get(Type::Int64Ty, 0);
else if (name == "double*")
return PointerType::get(Type::DoubleTy, 0);
+#if 0
+ else if (name == "complex*")
+ return ComplexPtrTy;
+#endif
else if (name == "expr*")
return ExprPtrTy;
else if (name == "expr**")
return ExprPtrPtrTy;
+ else if (name == "matrix*")
+ return GSLMatrixPtrTy;
+ else if (name == "dmatrix*")
+ return GSLDoubleMatrixPtrTy;
+ else if (name == "cmatrix*")
+ return GSLComplexMatrixPtrTy;
+ else if (name == "imatrix*")
+ return GSLIntMatrixPtrTy;
else if (name == "void*")
return VoidPtrTy;
else if (name.size() > 0 && name[name.size()-1] == '*')
@@ -3603,6 +3721,8 @@
return "float";
else if (type == Type::DoubleTy)
return "double";
+ else if (type == ComplexTy)
+ return "complex";
else if (type == CharPtrTy)
return "char*";
else if (type == PointerType::get(Type::Int16Ty, 0))
@@ -3613,10 +3733,20 @@
return "long*";
else if (type == PointerType::get(Type::DoubleTy, 0))
return "double*";
+ else if (type == ComplexPtrTy)
+ return "complex*";
else if (type == ExprPtrTy)
return "expr*";
else if (type == ExprPtrPtrTy)
return "expr**";
+ else if (type == GSLMatrixPtrTy)
+ return "matrix*";
+ else if (type == GSLDoubleMatrixPtrTy)
+ return "dmatrix*";
+ else if (type == GSLComplexMatrixPtrTy)
+ return "cmatrix*";
+ else if (type == GSLIntMatrixPtrTy)
+ return "imatrix*";
else if (type == VoidPtrTy)
return "void*";
else if (type->getTypeID() == Type::PointerTyID)
@@ -3794,7 +3924,6 @@
Value *x = args[i];
// check for thunks which must be forced
if (argt[i] != ExprPtrTy) {
-#if 1
// do a quick check on the tag value
Value *idx[2] = { Zero, Zero };
Value *tagv = b.CreateLoad(b.CreateGEP(x, idx, idx+2), "tag");
@@ -3808,9 +3937,6 @@
b.CreateBr(skipbb);
f->getBasicBlockList().push_back(skipbb);
b.SetInsertPoint(skipbb);
-#else
- b.CreateCall(module->getFunction("pure_force"), x);
-#endif
}
if (argt[i] == Type::Int1Ty) {
BasicBlock *okbb = BasicBlock::Create("ok");
@@ -3980,6 +4106,28 @@
idx[1] = ValFldIndex;
Value *ptrv = b.CreateLoad(b.CreateGEP(pv, idx, idx+2), "ptrval");
unboxed[i] = b.CreateBitCast(ptrv, argt[i]);
+ } else if (argt[i] == GSLMatrixPtrTy ||
+ argt[i] == GSLDoubleMatrixPtrTy ||
+ argt[i] == GSLComplexMatrixPtrTy ||
+ argt[i] == GSLIntMatrixPtrTy) {
+ BasicBlock *okbb = BasicBlock::Create("ok");
+ Value *idx[2] = { Zero, Zero };
+ Value *tagv = b.CreateLoad(b.CreateGEP(x, idx, idx+2), "tag");
+ int32_t ttag = -99;
+ if (argt[i] == GSLMatrixPtrTy)
+ ttag = EXPR::MATRIX;
+ else if (argt[i] == GSLDoubleMatrixPtrTy)
+ ttag = EXPR::DMATRIX;
+ else if (argt[i] == GSLComplexMatrixPtrTy)
+ ttag = EXPR::CMATRIX;
+ else if (argt[i] == GSLIntMatrixPtrTy)
+ ttag = EXPR::IMATRIX;
+ b.CreateCondBr
+ (b.CreateICmpEQ(tagv, SInt(ttag), "cmp"), okbb, failedbb);
+ f->getBasicBlockList().push_back(okbb);
+ b.SetInsertPoint(okbb);
+ Value *matv = b.CreateCall(module->getFunction("pure_get_matrix"), x);
+ unboxed[i] = b.CreateBitCast(matv, argt[i]);
} else if (argt[i] == ExprPtrTy) {
// passed through
unboxed[i] = x;
@@ -4061,6 +4209,18 @@
type == PointerType::get(Type::DoubleTy, 0))
u = b.CreateCall(module->getFunction("pure_pointer"),
b.CreateBitCast(u, VoidPtrTy));
+ else if (type == GSLMatrixPtrTy)
+ u = b.CreateCall(module->getFunction("pure_symbolic_matrix"),
+ b.CreateBitCast(u, VoidPtrTy));
+ else if (type == GSLDoubleMatrixPtrTy)
+ u = b.CreateCall(module->getFunction("pure_double_matrix"),
+ b.CreateBitCast(u, VoidPtrTy));
+ else if (type == GSLComplexMatrixPtrTy)
+ u = b.CreateCall(module->getFunction("pure_complex_matrix"),
+ b.CreateBitCast(u, VoidPtrTy));
+ else if (type == GSLIntMatrixPtrTy)
+ u = b.CreateCall(module->getFunction("pure_int_matrix"),
+ b.CreateBitCast(u, VoidPtrTy));
else if (type == ExprPtrTy) {
// check that we actually got a valid pointer; otherwise the call failed
BasicBlock *okbb = BasicBlock::Create("ok");
Modified: pure/trunk/interpreter.hh
===================================================================
--- pure/trunk/interpreter.hh 2008-09-20 20:53:03 UTC (rev 813)
+++ pure/trunk/interpreter.hh 2008-09-21 05:24:34 UTC (rev 814)
@@ -493,9 +493,13 @@
llvm::ExecutionEngine *JIT;
llvm::FunctionPassManager *FPM;
llvm::StructType *ExprTy, *IntExprTy, *DblExprTy, *StrExprTy, *PtrExprTy;
+ llvm::StructType *ComplexTy, *GSLMatrixTy, *GSLDoubleMatrixTy,
+ *GSLComplexMatrixTy, *GSLIntMatrixTy;
llvm::PointerType *ExprPtrTy, *ExprPtrPtrTy;
llvm::PointerType *IntExprPtrTy, *DblExprPtrTy, *StrExprPtrTy, *PtrExprPtrTy;
- llvm::PointerType *VoidPtrTy, *CharPtrTy;
+ llvm::PointerType *VoidPtrTy, *CharPtrTy, *IntPtrTy, *DoublePtrTy;
+ llvm::PointerType *ComplexPtrTy, *GSLMatrixPtrTy, *GSLDoubleMatrixPtrTy,
+ *GSLComplexMatrixPtrTy, *GSLIntMatrixPtrTy;
const llvm::Type *named_type(string name);
const char *type_name(const llvm::Type *type);
map<int32_t,GlobalVar> globalvars;
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-20 20:53:03 UTC (rev 813)
+++ pure/trunk/lib/primitives.pure 2008-09-21 05:24:34 UTC (rev 814)
@@ -117,7 +117,8 @@
pointer x::int |
pointer x::bigint |
pointer x::double |
-pointer x::string = pure_pointerval x;
+pointer x::string |
+pointer x::matrix = pure_pointerval x;
/* Convert signed (8/16/32/64) bit integers to the corresponding unsigned
quantities. These functions behave as if the value was "cast" to the
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-20 20:53:03 UTC (rev 813)
+++ pure/trunk/runtime.cc 2008-09-21 05:24:34 UTC (rev 814)
@@ -2522,6 +2522,14 @@
}
extern "C"
+void *pure_get_matrix(pure_expr *x)
+{
+ assert(x && x->tag == EXPR::MATRIX || x->tag == EXPR::DMATRIX ||
+ x->tag == EXPR::CMATRIX || x->tag == EXPR::IMATRIX);
+ return x->data.mat.p;
+}
+
+extern "C"
pure_expr *pure_matrix_rows(uint32_t n, ...)
{
va_list ap;
@@ -3486,6 +3494,10 @@
#else
return pure_pointer((void*)(uint32_t)pure_get_int(x));
#endif
+ case EXPR::MATRIX:
+ case EXPR::DMATRIX:
+ case EXPR::CMATRIX:
+ case EXPR::IMATRIX: return pure_pointer(x->data.mat.p);
default: return 0;
}
}
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-20 20:53:03 UTC (rev 813)
+++ pure/trunk/runtime.h 2008-09-21 05:24:34 UTC (rev 814)
@@ -438,6 +438,11 @@
int64_t pure_get_long(pure_expr *x);
int32_t pure_get_int(pure_expr *x);
+/* Convert a matrix expression to a pointer to the corresponding GSL matrix
+ struct. This is used to marshall matrix arguments in the C interface. */
+
+void *pure_get_matrix(pure_expr *x);
+
/* Additional matrix constructors. These work like pure_matrix_rowsl and
pure_matrix_columnsl in the public API, but are intended to be called
directly from generated code and raise the appropriate Pure exceptions in
@@ -546,7 +551,8 @@
an expression denoting an int, double, bigint or pointer value. The numeric
value of a pointer is its address, cast to a suitably large integer type,
which can be converted to/from an integer, but not a double value. Strings
- can be converted to pointers as well, but not the other way round. */
+ and matrices can be converted to pointers as well, but not the other way
+ round. */
pure_expr *pure_intval(pure_expr *x);
pure_expr *pure_dblval(pure_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-21 05:29:17
|
Revision: 815
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=815&view=rev
Author: agraef
Date: 2008-09-21 05:29:14 +0000 (Sun, 21 Sep 2008)
Log Message:
-----------
Updated NEWS and ChangeLog.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/NEWS
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-21 05:24:34 UTC (rev 814)
+++ pure/trunk/ChangeLog 2008-09-21 05:29:14 UTC (rev 815)
@@ -3,7 +3,9 @@
* Implemented basic GSL matrix support, including support for
symbolic matrices (which is independent from GSL, so these will
also work when building the interpreter without GSL) and matrix
- comprehensions. This required many additions and changes to the
+ comprehensions. Marshalling of matrices in the C interface is also
+ implemented, so that you can interface to GSL matrix functions
+ without much ado. This required many additions and changes to the
parser, interpreter, compiler, runtime and the prelude; details
can be found in the svn log (see r759 and r769ff.).
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-21 05:24:34 UTC (rev 814)
+++ pure/trunk/NEWS 2008-09-21 05:29:14 UTC (rev 815)
@@ -67,7 +67,10 @@
type-checking and various conversions between the different kinds of
matrices. To these ends, various basic operations to deal with matrix objects
have been added to the runtime, including some public API operations to create
-and inspect Pure matrix values from external C modules.
+and inspect Pure matrix values from external C modules. Moreover, the
+'pointer' function in the prelude can be used to convert matrices to Pure
+pointer values, and marshalling of GSL matrices in the C interface is also
+available.
Here is a brief synopsis of some important operations which are already
implemented in the prelude (you can find the definitions of these and the
@@ -81,6 +84,8 @@
- x!!is, x!!(is,js) slicing (is and js are lists of machine ints)
- x==y, x!=y matrix comparisons
+Adding other operations by interfacing to GSL should be a piece of cake.
+
Other user-visible changes prompted by the introduction of the matrix syntax
are listed below:
@@ -96,11 +101,6 @@
lower precedence than the ',' operator, which makes writing matrix slices like
x!!(i..j,k..l) much more convenient.
-Stuff that's still missing:
-
-- Marshalling of GSL matrices in the C interface, so that it becomes possible
-to access all the advanced functionality available in GSL.
-
** Pure 0.6 2008-09-12
New stuff in this release (please see the ChangeLog and the manual for
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-21 07:44:08
|
Revision: 816
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=816&view=rev
Author: agraef
Date: 2008-09-21 07:44:02 +0000 (Sun, 21 Sep 2008)
Log Message:
-----------
Experimental support for marshalling complex numbers (this doesn't work right now, and so is disabled).
Modified Paths:
--------------
pure/trunk/interpreter.cc
pure/trunk/interpreter.hh
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-21 05:29:14 UTC (rev 815)
+++ pure/trunk/interpreter.cc 2008-09-21 07:44:02 UTC (rev 816)
@@ -133,8 +133,7 @@
// Complex numbers (complex double).
{
std::vector<const Type*> elts;
- elts.push_back(Type::DoubleTy);
- elts.push_back(Type::DoubleTy);
+ elts.push_back(ArrayType::get(Type::DoubleTy, 2));
ComplexTy = StructType::get(elts);
ComplexPtrTy = PointerType::get(ComplexTy, 0);
}
@@ -376,6 +375,16 @@
declare_extern((void*)pure_get_matrix,
"pure_get_matrix", "void*", 1, "expr*");
+#if COMPLEX_NUMBERS
+ /* Marshalling of complex numbers. This doesn't work yet and is disabled. */
+ declare_extern((void*)pure_complex,
+ "pure_complex", "expr*", 1, "complex");
+ declare_extern((void*)pure_is_complex,
+ "pure_is_complex", "bool", 1, "expr*");
+ declare_extern((void*)pure_get_complex,
+ "pure_get_complex", "complex", 1, "expr*");
+#endif
+
declare_extern((void*)pure_catch,
"pure_catch", "expr*", 2, "expr*", "expr*");
declare_extern((void*)pure_throw,
@@ -3664,7 +3673,7 @@
return Type::FloatTy;
else if (name == "double")
return Type::DoubleTy;
-#if 0 // no marshalling available yet, does LLVM support these?
+#if COMPLEX_NUMBERS
else if (name == "complex")
return ComplexTy;
#endif
@@ -3678,7 +3687,7 @@
return PointerType::get(Type::Int64Ty, 0);
else if (name == "double*")
return PointerType::get(Type::DoubleTy, 0);
-#if 0
+#if COMPLEX_NUMBERS
else if (name == "complex*")
return ComplexPtrTy;
#endif
@@ -3721,8 +3730,10 @@
return "float";
else if (type == Type::DoubleTy)
return "double";
+#if COMPLEX_NUMBERS
else if (type == ComplexTy)
return "complex";
+#endif
else if (type == CharPtrTy)
return "char*";
else if (type == PointerType::get(Type::Int16Ty, 0))
@@ -3733,8 +3744,10 @@
return "long*";
else if (type == PointerType::get(Type::DoubleTy, 0))
return "double*";
+#if COMPLEX_NUMBERS
else if (type == ComplexPtrTy)
return "complex*";
+#endif
else if (type == ExprPtrTy)
return "expr*";
else if (type == ExprPtrPtrTy)
@@ -4080,6 +4093,19 @@
idx[1] = ValFldIndex;
Value *dv = b.CreateLoad(b.CreateGEP(pv, idx, idx+2), "dblval");
unboxed[i] = dv;
+#if COMPLEX_NUMBERS
+ } else if (argt[i] == ComplexTy) {
+ BasicBlock *okbb = BasicBlock::Create("ok");
+ // Pure's complex values are a special algebraic data type defined in
+ // math.pure, hence we have to go to some lengths here to get these
+ // values.
+ Value *chk = b.CreateCall(module->getFunction("pure_is_complex"), x);
+ b.CreateCondBr(chk, okbb, failedbb);
+ f->getBasicBlockList().push_back(okbb);
+ b.SetInsertPoint(okbb);
+ Value *cv = b.CreateCall(module->getFunction("pure_get_complex"), x);
+ unboxed[i] = cv;
+#endif
} else if (argt[i] == CharPtrTy) {
BasicBlock *okbb = BasicBlock::Create("ok");
Value *idx[2] = { Zero, Zero };
@@ -4094,7 +4120,11 @@
argt[i] == PointerType::get(Type::Int32Ty, 0) ||
argt[i] == PointerType::get(Type::Int64Ty, 0) ||
argt[i] == PointerType::get(Type::FloatTy, 0) ||
- argt[i] == PointerType::get(Type::DoubleTy, 0)) {
+ argt[i] == PointerType::get(Type::DoubleTy, 0)
+#if COMPLEX_NUMBERS
+ || argt[i] == ComplexPtrTy
+#endif
+ ) {
BasicBlock *okbb = BasicBlock::Create("ok");
Value *idx[2] = { Zero, Zero };
Value *tagv = b.CreateLoad(b.CreateGEP(x, idx, idx+2), "tag");
@@ -4200,13 +4230,21 @@
b.CreateFPExt(u, Type::DoubleTy));
else if (type == Type::DoubleTy)
u = b.CreateCall(module->getFunction("pure_double"), u);
+#if COMPLEX_NUMBERS
+ else if (type == ComplexTy)
+ u = b.CreateCall(module->getFunction("pure_complex"), u);
+#endif
else if (type == CharPtrTy)
u = b.CreateCall(module->getFunction("pure_cstring_dup"), u);
else if (type == PointerType::get(Type::Int16Ty, 0) ||
type == PointerType::get(Type::Int32Ty, 0) ||
type == PointerType::get(Type::Int64Ty, 0) ||
type == PointerType::get(Type::FloatTy, 0) ||
- type == PointerType::get(Type::DoubleTy, 0))
+ type == PointerType::get(Type::DoubleTy, 0)
+#if COMPLEX_NUMBERS
+ || type == ComplexPtrTy
+#endif
+ )
u = b.CreateCall(module->getFunction("pure_pointer"),
b.CreateBitCast(u, VoidPtrTy));
else if (type == GSLMatrixPtrTy)
Modified: pure/trunk/interpreter.hh
===================================================================
--- pure/trunk/interpreter.hh 2008-09-21 05:29:14 UTC (rev 815)
+++ pure/trunk/interpreter.hh 2008-09-21 07:44:02 UTC (rev 816)
@@ -41,6 +41,13 @@
default, use 0 to disable this option). */
#define LIST_KLUDGE 10
+/* Experimental support for marshalling of complex numbers in the C
+ interface. This doesn't work right now and is disabled. LLVM doesn't seem
+ to provide a transparent way to handle complex values in function calls
+ yet, and maybe this isn't even possible because different compilers might
+ specify different ABIs to do that kind of thing. */
+#define COMPLEX_NUMBERS 0
+
using namespace std;
/* The Pure interpreter. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-21 05:29:14 UTC (rev 815)
+++ pure/trunk/runtime.cc 2008-09-21 07:44:02 UTC (rev 816)
@@ -890,6 +890,28 @@
#endif
}
+static inline bool is_complex(pure_expr *x)
+{
+ 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))
+ return false;
+ u = u->data.x[1];
+ if (u->tag != EXPR::INT && u->tag != EXPR::DBL ||
+ v->tag != EXPR::INT && v->tag != EXPR::DBL)
+ return false;
+ else
+ return true;
+ } else
+ return false;
+}
+
static inline bool get_complex(pure_expr *x, double& a, double& b)
{
if (x->tag != EXPR::APP) return false;
@@ -2521,7 +2543,33 @@
return &x->data.z;
}
+#if COMPLEX_NUMBERS
extern "C"
+pure_expr *pure_complex(__complex_double c)
+{
+ return make_complex(c.x[0], c.x[1]);
+}
+
+extern "C"
+bool pure_is_complex(pure_expr *x)
+{
+ return is_complex(x);
+}
+
+extern "C"
+__complex_double pure_get_complex(pure_expr *x)
+{
+ __complex_double res = {{0.0,0.0}};
+ double a, b;
+ if (get_complex(x, a, b)) {
+ res.x[0] = a;
+ res.x[1] = b;
+ }
+ return res;
+}
+#endif
+
+extern "C"
void *pure_get_matrix(pure_expr *x)
{
assert(x && x->tag == EXPR::MATRIX || x->tag == EXPR::DMATRIX ||
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-21 05:29:14 UTC (rev 815)
+++ pure/trunk/runtime.h 2008-09-21 07:44:02 UTC (rev 816)
@@ -438,6 +438,21 @@
int64_t pure_get_long(pure_expr *x);
int32_t pure_get_int(pure_expr *x);
+#if 0
+// This stuff is disabled right now as it doesn't work yet.
+/* Marshall complex numbers. These are actually defined as an algebraic type
+ in math.pure, but we need some basic support for these values in the C
+ interface. */
+
+/* We don't want to rely on ISO complex number support here. The following
+ should do the job on all supported systems. */
+typedef struct { double x[2]; } __complex_double;
+
+pure_expr *pure_complex(__complex_double c);
+bool pure_is_complex(pure_expr *x);
+__complex_double pure_get_complex(pure_expr *x);
+#endif
+
/* Convert a matrix expression to a pointer to the corresponding GSL matrix
struct. This is used to marshall matrix arguments in the C interface. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-21 12:42:44
|
Revision: 817
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=817&view=rev
Author: agraef
Date: 2008-09-21 12:42:34 +0000 (Sun, 21 Sep 2008)
Log Message:
-----------
Eliminate broken complex number marshalling support. No easy way to make this work.
Modified Paths:
--------------
pure/trunk/interpreter.cc
pure/trunk/interpreter.hh
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/interpreter.cc
===================================================================
--- pure/trunk/interpreter.cc 2008-09-21 07:44:02 UTC (rev 816)
+++ pure/trunk/interpreter.cc 2008-09-21 12:42:34 UTC (rev 817)
@@ -133,7 +133,9 @@
// Complex numbers (complex double).
{
std::vector<const Type*> elts;
- elts.push_back(ArrayType::get(Type::DoubleTy, 2));
+ //elts.push_back(ArrayType::get(Type::DoubleTy, 2));
+ elts.push_back(Type::DoubleTy);
+ elts.push_back(Type::DoubleTy);
ComplexTy = StructType::get(elts);
ComplexPtrTy = PointerType::get(ComplexTy, 0);
}
@@ -375,16 +377,6 @@
declare_extern((void*)pure_get_matrix,
"pure_get_matrix", "void*", 1, "expr*");
-#if COMPLEX_NUMBERS
- /* Marshalling of complex numbers. This doesn't work yet and is disabled. */
- declare_extern((void*)pure_complex,
- "pure_complex", "expr*", 1, "complex");
- declare_extern((void*)pure_is_complex,
- "pure_is_complex", "bool", 1, "expr*");
- declare_extern((void*)pure_get_complex,
- "pure_get_complex", "complex", 1, "expr*");
-#endif
-
declare_extern((void*)pure_catch,
"pure_catch", "expr*", 2, "expr*", "expr*");
declare_extern((void*)pure_throw,
@@ -3673,10 +3665,6 @@
return Type::FloatTy;
else if (name == "double")
return Type::DoubleTy;
-#if COMPLEX_NUMBERS
- else if (name == "complex")
- return ComplexTy;
-#endif
else if (name == "char*")
return CharPtrTy;
else if (name == "short*")
@@ -3687,10 +3675,6 @@
return PointerType::get(Type::Int64Ty, 0);
else if (name == "double*")
return PointerType::get(Type::DoubleTy, 0);
-#if COMPLEX_NUMBERS
- else if (name == "complex*")
- return ComplexPtrTy;
-#endif
else if (name == "expr*")
return ExprPtrTy;
else if (name == "expr**")
@@ -3730,10 +3714,6 @@
return "float";
else if (type == Type::DoubleTy)
return "double";
-#if COMPLEX_NUMBERS
- else if (type == ComplexTy)
- return "complex";
-#endif
else if (type == CharPtrTy)
return "char*";
else if (type == PointerType::get(Type::Int16Ty, 0))
@@ -3744,10 +3724,6 @@
return "long*";
else if (type == PointerType::get(Type::DoubleTy, 0))
return "double*";
-#if COMPLEX_NUMBERS
- else if (type == ComplexPtrTy)
- return "complex*";
-#endif
else if (type == ExprPtrTy)
return "expr*";
else if (type == ExprPtrPtrTy)
@@ -4093,19 +4069,6 @@
idx[1] = ValFldIndex;
Value *dv = b.CreateLoad(b.CreateGEP(pv, idx, idx+2), "dblval");
unboxed[i] = dv;
-#if COMPLEX_NUMBERS
- } else if (argt[i] == ComplexTy) {
- BasicBlock *okbb = BasicBlock::Create("ok");
- // Pure's complex values are a special algebraic data type defined in
- // math.pure, hence we have to go to some lengths here to get these
- // values.
- Value *chk = b.CreateCall(module->getFunction("pure_is_complex"), x);
- b.CreateCondBr(chk, okbb, failedbb);
- f->getBasicBlockList().push_back(okbb);
- b.SetInsertPoint(okbb);
- Value *cv = b.CreateCall(module->getFunction("pure_get_complex"), x);
- unboxed[i] = cv;
-#endif
} else if (argt[i] == CharPtrTy) {
BasicBlock *okbb = BasicBlock::Create("ok");
Value *idx[2] = { Zero, Zero };
@@ -4120,11 +4083,7 @@
argt[i] == PointerType::get(Type::Int32Ty, 0) ||
argt[i] == PointerType::get(Type::Int64Ty, 0) ||
argt[i] == PointerType::get(Type::FloatTy, 0) ||
- argt[i] == PointerType::get(Type::DoubleTy, 0)
-#if COMPLEX_NUMBERS
- || argt[i] == ComplexPtrTy
-#endif
- ) {
+ argt[i] == PointerType::get(Type::DoubleTy, 0)) {
BasicBlock *okbb = BasicBlock::Create("ok");
Value *idx[2] = { Zero, Zero };
Value *tagv = b.CreateLoad(b.CreateGEP(x, idx, idx+2), "tag");
@@ -4230,21 +4189,13 @@
b.CreateFPExt(u, Type::DoubleTy));
else if (type == Type::DoubleTy)
u = b.CreateCall(module->getFunction("pure_double"), u);
-#if COMPLEX_NUMBERS
- else if (type == ComplexTy)
- u = b.CreateCall(module->getFunction("pure_complex"), u);
-#endif
else if (type == CharPtrTy)
u = b.CreateCall(module->getFunction("pure_cstring_dup"), u);
else if (type == PointerType::get(Type::Int16Ty, 0) ||
type == PointerType::get(Type::Int32Ty, 0) ||
type == PointerType::get(Type::Int64Ty, 0) ||
type == PointerType::get(Type::FloatTy, 0) ||
- type == PointerType::get(Type::DoubleTy, 0)
-#if COMPLEX_NUMBERS
- || type == ComplexPtrTy
-#endif
- )
+ type == PointerType::get(Type::DoubleTy, 0))
u = b.CreateCall(module->getFunction("pure_pointer"),
b.CreateBitCast(u, VoidPtrTy));
else if (type == GSLMatrixPtrTy)
Modified: pure/trunk/interpreter.hh
===================================================================
--- pure/trunk/interpreter.hh 2008-09-21 07:44:02 UTC (rev 816)
+++ pure/trunk/interpreter.hh 2008-09-21 12:42:34 UTC (rev 817)
@@ -41,13 +41,6 @@
default, use 0 to disable this option). */
#define LIST_KLUDGE 10
-/* Experimental support for marshalling of complex numbers in the C
- interface. This doesn't work right now and is disabled. LLVM doesn't seem
- to provide a transparent way to handle complex values in function calls
- yet, and maybe this isn't even possible because different compilers might
- specify different ABIs to do that kind of thing. */
-#define COMPLEX_NUMBERS 0
-
using namespace std;
/* The Pure interpreter. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-21 07:44:02 UTC (rev 816)
+++ pure/trunk/runtime.cc 2008-09-21 12:42:34 UTC (rev 817)
@@ -2048,6 +2048,26 @@
}
extern "C"
+pure_expr *pure_complex(double c[2])
+{
+ return make_complex(c[0], c[1]);
+}
+
+extern "C"
+bool pure_is_complex(pure_expr *x, double *c)
+{
+ double a, b;
+ if (get_complex(x, a, b)) {
+ if (c) {
+ c[0] = a;
+ c[1] = b;
+ }
+ return true;
+ } else
+ return false;
+}
+
+extern "C"
pure_expr *pure_new(pure_expr *x)
{
return pure_new_internal(x);
@@ -2543,33 +2563,7 @@
return &x->data.z;
}
-#if COMPLEX_NUMBERS
extern "C"
-pure_expr *pure_complex(__complex_double c)
-{
- return make_complex(c.x[0], c.x[1]);
-}
-
-extern "C"
-bool pure_is_complex(pure_expr *x)
-{
- return is_complex(x);
-}
-
-extern "C"
-__complex_double pure_get_complex(pure_expr *x)
-{
- __complex_double res = {{0.0,0.0}};
- double a, b;
- if (get_complex(x, a, b)) {
- res.x[0] = a;
- res.x[1] = b;
- }
- return res;
-}
-#endif
-
-extern "C"
void *pure_get_matrix(pure_expr *x)
{
assert(x && x->tag == EXPR::MATRIX || x->tag == EXPR::DMATRIX ||
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-21 07:44:02 UTC (rev 816)
+++ pure/trunk/runtime.h 2008-09-21 12:42:34 UTC (rev 817)
@@ -290,6 +290,15 @@
bool pure_is_listv(pure_expr *x, size_t *size, pure_expr ***elems);
bool pure_is_tuplev(pure_expr *x, size_t *size, pure_expr ***elems);
+/* Complex number support. These are actually defined as an algebraic type in
+ math.pure, but some applications may require access to these values in the
+ C interface. Note that we don't want to rely on ISOC99 complex number
+ support or any kind of third-party library here, so this API represents
+ complex values simply as a double[2]. */
+
+pure_expr *pure_complex(double c[2]);
+bool pure_is_complex(pure_expr *x, double *c);
+
/* Memory management. */
/* Count a new reference to an expression. This should be called whenever you
@@ -438,21 +447,6 @@
int64_t pure_get_long(pure_expr *x);
int32_t pure_get_int(pure_expr *x);
-#if 0
-// This stuff is disabled right now as it doesn't work yet.
-/* Marshall complex numbers. These are actually defined as an algebraic type
- in math.pure, but we need some basic support for these values in the C
- interface. */
-
-/* We don't want to rely on ISO complex number support here. The following
- should do the job on all supported systems. */
-typedef struct { double x[2]; } __complex_double;
-
-pure_expr *pure_complex(__complex_double c);
-bool pure_is_complex(pure_expr *x);
-__complex_double pure_get_complex(pure_expr *x);
-#endif
-
/* Convert a matrix expression to a pointer to the corresponding GSL matrix
struct. This is used to marshall matrix arguments in the C interface. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-21 22:53:34
|
Revision: 820
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=820&view=rev
Author: agraef
Date: 2008-09-21 22:53:28 +0000 (Sun, 21 Sep 2008)
Log Message:
-----------
Add function to determine the stride of a matrix.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-21 22:32:31 UTC (rev 819)
+++ pure/trunk/lib/primitives.pure 2008-09-21 22:53:28 UTC (rev 820)
@@ -412,11 +412,13 @@
/* Basic matrix operations: size, dimensions and indexing. */
-private matrix_size matrix_dim;
-extern int matrix_size(expr *x), expr* matrix_dim(expr *x);
+private matrix_size matrix_dim matrix_stride;
+extern int matrix_size(expr *x), int matrix_stride(expr *x);
+extern expr* matrix_dim(expr *x);
#x::matrix = matrix_size x;
dim x::matrix = matrix_dim x;
+stride x::matrix = matrix_stride x;
private matrix_elem_at matrix_elem_at2;
extern expr* matrix_elem_at(expr* x, int i);
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-21 22:32:31 UTC (rev 819)
+++ pure/trunk/runtime.cc 2008-09-21 22:53:28 UTC (rev 820)
@@ -4096,6 +4096,33 @@
}
extern "C"
+uint32_t matrix_stride(pure_expr *x)
+{
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ return m->tda;
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return m->tda;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return m->tda;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return m->tda;
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
int32_t matrix_type(pure_expr *x)
{
uint32_t t = (uint32_t)x->tag;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-21 22:32:31 UTC (rev 819)
+++ pure/trunk/runtime.h 2008-09-21 22:53:28 UTC (rev 820)
@@ -660,13 +660,19 @@
/* 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). matrix_type
- determines the exact type of a matrix, returning an integer denoting the
- subtype tag (0 = symbolic, 1 = double, 2 = complex, 3 = integer matrix; -1
- is returned if the given object is not a matrix). */
+ number of rows and columns, which are returned as a pair (n,m).
+ matrix_stride returns the real row length of a matrix, which may be larger
+ than its column count if the matrix is actually a slice of a larger matrix.
+ (This value won't be of much use in Pure programs since there is no way to
+ access the "extra" elements in each row, but may be useful if the data
+ pointer is passed to an external C routine.) matrix_type determines the
+ exact type of a matrix, returning an integer denoting the subtype tag (0 =
+ symbolic, 1 = double, 2 = complex, 3 = integer matrix; -1 is returned if
+ the given object is not a matrix). */
uint32_t matrix_size(pure_expr *x);
pure_expr *matrix_dim(pure_expr *x);
+uint32_t matrix_stride(pure_expr *x);
int32_t matrix_type(pure_expr *x);
/* Matrix elements can be retrieved either by a single index (using row-major
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-22 19:20:32
|
Revision: 828
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=828&view=rev
Author: agraef
Date: 2008-09-22 19:20:24 +0000 (Mon, 22 Sep 2008)
Log Message:
-----------
Add operation to change the dimensions of a matrix without changing its size.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-22 17:46:47 UTC (rev 827)
+++ pure/trunk/lib/primitives.pure 2008-09-22 19:20:24 UTC (rev 828)
@@ -447,7 +447,8 @@
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
+x::matrix!!ns = if packed x then rowvector x!!(1,ns)
+ else colcatmap (nth x) ns with
nth x n = catch (cst {}) {x!n};
end;
@@ -468,11 +469,28 @@
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. */
+/* Extract a submatrix of a given size at a given offset. The result shares
+ the underlying storage with the input matrix (i.e., matrix elements are
+ *not* copied) and so this is a comparatively cheap operation. */
submat x::matrix (i::int,j::int) (n::int,m::int)
= matrix_slice x i j (i+n-1) (j+m-1);
+/* Change the dimensions of a matrix without changing its size. The total
+ number of elements must match that of the input matrix. Reuses the
+ underlying storage of the input matrix if possible. */
+
+private matrix_redim;
+extern expr* matrix_redim(expr* x, int n, int m);
+
+redim x::matrix (n::int,m::int)
+ = matrix_redim x n m if n*m==#x;
+
+/* Convenience functions for converting a matrix to a row or column vector. */
+
+rowvector x::matrix = redim x (1,#x);
+colvector x::matrix = redim x (#x,1);
+
/* 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
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-22 17:46:47 UTC (rev 827)
+++ pure/trunk/runtime.cc 2008-09-22 19:20:24 UTC (rev 828)
@@ -4313,6 +4313,114 @@
}
extern "C"
+pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m)
+{
+ void *p = 0;
+ if (n<0 || m<0) return 0;
+ const size_t n1 = (size_t)n, n2 = (size_t)m;
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix_symbolic *m1 =
+ (gsl_matrix_symbolic*)malloc(sizeof(gsl_matrix_symbolic));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_symbolic_matrix_dup(m);
+ if (y) {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix *m1 = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_double_matrix_dup(m);
+ if (y) {
+ gsl_matrix *m = (gsl_matrix*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix_complex *m1 =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_complex_matrix_dup(m);
+ if (y) {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix_int *m1 =
+ (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_int_matrix_dup(m);
+ if (y) {
+ gsl_matrix_int *m = (gsl_matrix_int*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+#endif
+ default:
+ return 0;
+ }
+ pure_expr *y = new_expr();
+ y->tag = x->tag;
+ y->data.mat.p = p;
+ y->data.mat.refc = x->data.mat.refc;
+ (*y->data.mat.refc)++;
+ MEMDEBUG_NEW(y)
+ return y;
+}
+
+extern "C"
pure_expr *matrix_rows(pure_expr *xs)
{
size_t n;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-22 17:46:47 UTC (rev 827)
+++ pure/trunk/runtime.h 2008-09-22 19:20:24 UTC (rev 828)
@@ -692,6 +692,11 @@
pure_expr *matrix_slice(pure_expr *x, int32_t i1, int32_t j1,
int32_t i2, int32_t j2);
+/* Convert a matrix to a new matrix with same size but different dimensions.
+ Reuse the storage of the original matrix if its data is contiguous. */
+
+pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m);
+
/* Matrix construction. These work like the pure_matrix_rows/
pure_matrix_columns functions in the public API, but take their input from
a Pure list. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-22 20:53:45
|
Revision: 829
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=829&view=rev
Author: agraef
Date: 2008-09-22 20:53:37 +0000 (Mon, 22 Sep 2008)
Log Message:
-----------
Add (sub-,super-) diagonal operations.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-22 19:20:24 UTC (rev 828)
+++ pure/trunk/lib/primitives.pure 2008-09-22 20:53:37 UTC (rev 829)
@@ -469,6 +469,16 @@
cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
+/* Extract (sub-,super-) diagonals from a matrix. Sub- and super-diagonals for
+ k=0 return the main diagonal. Indices for sub- and super-diagonals can also
+ be negative, in which case the corresponding super- or sub-diagonal is
+ returned instead. */
+
+private matrix_diag matrix_subdiag matrix_supdiag;
+extern expr* matrix_diag(expr* x) = diag;
+extern expr* matrix_subdiag(expr* x, int k) = subdiag;
+extern expr* matrix_supdiag(expr* x, int k) = supdiag;
+
/* Extract a submatrix of a given size at a given offset. The result shares
the underlying storage with the input matrix (i.e., matrix elements are
*not* copied) and so this is a comparatively cheap operation. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-22 19:20:24 UTC (rev 828)
+++ pure/trunk/runtime.cc 2008-09-22 20:53:37 UTC (rev 829)
@@ -4313,6 +4313,154 @@
}
extern "C"
+pure_expr *matrix_diag(pure_expr *x)
+{
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix_symbolic *m1 = create_symbolic_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[i*(m->tda+1)];
+ return pure_symbolic_matrix(m1);
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix *m1 = create_double_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[i*(m->tda+1)];
+ return pure_double_matrix(m1);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix_complex *m1 = create_complex_matrix(1, n);
+ for (size_t i = 0; i < n; i++) {
+ const size_t k = 2*i*(m->tda+1);
+ m1->data[2*i] = m->data[k];
+ m1->data[2*i+1] = m->data[k+1];
+ }
+ return pure_complex_matrix(m1);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix_int *m1 = create_int_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[i*(m->tda+1)];
+ return pure_int_matrix(m1);
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
+pure_expr *matrix_subdiag(pure_expr *x, int32_t k)
+{
+ if (k<0) return matrix_supdiag(x, -k);
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix_symbolic *m1 = create_symbolic_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_symbolic_matrix(m1);
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix *m1 = create_double_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_double_matrix(m1);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix_complex *m1 = create_complex_matrix(1, n);
+ for (size_t i = 0; i < n; i++) {
+ const size_t k = 2*k0+2*i*(m->tda+1);
+ m1->data[2*i] = m->data[k];
+ m1->data[2*i+1] = m->data[k+1];
+ }
+ return pure_complex_matrix(m1);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix_int *m1 = create_int_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_int_matrix(m1);
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
+pure_expr *matrix_supdiag(pure_expr *x, int32_t k)
+{
+ if (k<0) return matrix_subdiag(x, -k);
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix_symbolic *m1 = create_symbolic_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_symbolic_matrix(m1);
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix *m1 = create_double_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_double_matrix(m1);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix_complex *m1 = create_complex_matrix(1, n);
+ for (size_t i = 0; i < n; i++) {
+ const size_t k = 2*k0+2*i*(m->tda+1);
+ m1->data[2*i] = m->data[k];
+ m1->data[2*i+1] = m->data[k+1];
+ }
+ return pure_complex_matrix(m1);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix_int *m1 = create_int_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_int_matrix(m1);
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m)
{
void *p = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-22 19:20:24 UTC (rev 828)
+++ pure/trunk/runtime.h 2008-09-22 20:53:37 UTC (rev 829)
@@ -697,6 +697,12 @@
pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m);
+/* Retrieve (sub-,super-)diagonals of a matrix, as a 1xn matrix. */
+
+pure_expr *matrix_diag(pure_expr *x);
+pure_expr *matrix_subdiag(pure_expr *x, int32_t k);
+pure_expr *matrix_supdiag(pure_expr *x, int32_t k);
+
/* Matrix construction. These work like the pure_matrix_rows/
pure_matrix_columns functions in the public API, but take their input from
a Pure list. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 05:54:41
|
Revision: 833
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=833&view=rev
Author: agraef
Date: 2008-09-23 05:54:38 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Make the matrix construction functions in the library throw a 'bad_matrix_value' exception in case of dimension mismatch.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-23 05:31:30 UTC (rev 832)
+++ pure/trunk/lib/primitives.pure 2008-09-23 05:54:38 UTC (rev 833)
@@ -473,7 +473,9 @@
converts it to a list of its elements; otherwise the result is the list of
the rows of the matrix. You can also use list2 to convert a matrix to a
list of lists. Conversely, matrix xs converts a list of lists or matrices
- to the corresponding matrix. Otherwise, the result is a row vector. */
+ to the corresponding matrix; otherwise, the result is a row vector.
+ NOTE: The matrix function may throw a 'bad_matrix_value x' in case of
+ dimension mismatch, where x denotes the offending submatrix. */
list x::matrix = [x!i|i=0..#x-1] if rowvectorp x;
= rows x otherwise;
@@ -520,7 +522,9 @@
/* 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}. */
+ horizontally, like {x,y}. NOTE: Like the built-in matrix constructs, these
+ operations may throw a 'bad_matrix_value x' in case of dimension mismatch,
+ where x denotes the offending submatrix. */
extern expr* matrix_rows(expr *x) = rowcat;
extern expr* matrix_columns(expr *x) = colcat;
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 05:31:30 UTC (rev 832)
+++ pure/trunk/runtime.cc 2008-09-23 05:54:38 UTC (rev 833)
@@ -4568,13 +4568,227 @@
return y;
}
+static pure_expr *matrix_rowsv(uint32_t n, pure_expr **xs)
+{
+ int k = -1;
+ size_t nrows = 0, ncols = 0;
+ int32_t target = 0;
+ bool have_matrix = false;
+ pure_expr *x = 0;
+ for (size_t i = 0; i < n; i++) {
+ x = xs[i];
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *mp = (gsl_matrix_symbolic*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ goto err;
+ nrows += mp->size1; k = mp->size2;
+ set_target_type(target, EXPR::MATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+#ifdef HAVE_GSL
+ case EXPR::DBL:
+ set_target_type(target, EXPR::DMATRIX);
+ if (k >= 0 && k != 1) goto err;
+ nrows++; k = 1;
+ break;
+ case EXPR::INT:
+ set_target_type(target, EXPR::IMATRIX);
+ if (k >= 0 && k != 1) goto err;
+ nrows++; k = 1;
+ break;
+ case EXPR::APP: {
+ double a, b;
+ if (k >= 0 && k != 1) goto err;
+ nrows++; k = 1;
+ if (get_complex(x, a, b))
+ set_target_type(target, EXPR::CMATRIX);
+ else
+ set_target_type(target, EXPR::MATRIX);
+ break;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ goto err;
+ nrows += mp->size1; k = mp->size2;
+ set_target_type(target, EXPR::DMATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ goto err;
+ nrows += mp->size1; k = mp->size2;
+ set_target_type(target, EXPR::CMATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size2 != (size_t)k)
+ goto err;
+ nrows += mp->size1; k = mp->size2;
+ set_target_type(target, EXPR::IMATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+#endif
+ default:
+ if (k >= 0 && k != 1) goto err;
+ nrows++; k = 1;
+ set_target_type(target, EXPR::MATRIX);
+ break;
+ }
+ }
+ if (n == 1 && have_matrix) return xs[0];
+ if (k < 0) k = 0;
+ ncols = k;
+ if (target == 0) target = EXPR::MATRIX;
+ switch (target) {
+ case EXPR::MATRIX:
+ return symbolic_matrix_rows(nrows, ncols, n, xs);
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX:
+ 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);
+#endif
+ default:
+ assert(0 && "this can't happen");
+ return 0;
+ }
+ err:
+ pure_throw(bad_matrix_exception(x));
+ return 0;
+}
+
+static pure_expr *matrix_columnsv(uint32_t n, pure_expr **xs)
+{
+ int k = -1;
+ size_t nrows = 0, ncols = 0;
+ int32_t target = 0;
+ bool have_matrix = false;
+ pure_expr *x = 0;
+ for (size_t i = 0; i < n; i++) {
+ x = xs[i];
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *mp = (gsl_matrix_symbolic*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ goto err;
+ ncols += mp->size2; k = mp->size1;
+ set_target_type(target, EXPR::MATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+#ifdef HAVE_GSL
+ case EXPR::DBL:
+ set_target_type(target, EXPR::DMATRIX);
+ if (k >= 0 && k != 1) goto err;
+ ncols++; k = 1;
+ break;
+ case EXPR::INT:
+ set_target_type(target, EXPR::IMATRIX);
+ if (k >= 0 && k != 1) goto err;
+ ncols++; k = 1;
+ break;
+ case EXPR::APP: {
+ double a, b;
+ if (k >= 0 && k != 1) goto err;
+ ncols++; k = 1;
+ if (get_complex(x, a, b))
+ set_target_type(target, EXPR::CMATRIX);
+ else
+ set_target_type(target, EXPR::MATRIX);
+ break;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *mp = (gsl_matrix*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ goto err;
+ ncols += mp->size2; k = mp->size1;
+ set_target_type(target, EXPR::DMATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *mp = (gsl_matrix_complex*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ goto err;
+ ncols += mp->size2; k = mp->size1;
+ set_target_type(target, EXPR::CMATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *mp = (gsl_matrix_int*)x->data.mat.p;
+ if (mp->size1 > 0 && mp->size2 > 0) {
+ if (k >= 0 && mp->size1 != (size_t)k)
+ goto err;
+ ncols += mp->size2; k = mp->size1;
+ set_target_type(target, EXPR::IMATRIX);
+ have_matrix = true;
+ }
+ break;
+ }
+#endif
+ default:
+ if (k >= 0 && k != 1) goto err;
+ ncols++; k = 1;
+ set_target_type(target, EXPR::MATRIX);
+ break;
+ }
+ }
+ if (n == 1 && have_matrix) return xs[0];
+ if (k < 0) k = 0;
+ nrows = k;
+ if (target == 0) target = EXPR::MATRIX;
+ switch (target) {
+ case EXPR::MATRIX:
+ return symbolic_matrix_columns(nrows, ncols, n, xs);
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX:
+ 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);
+#endif
+ default:
+ assert(0 && "this can't happen");
+ return 0;
+ }
+ err:
+ pure_throw(bad_matrix_exception(x));
+ return 0;
+}
+
extern "C"
pure_expr *matrix_rows(pure_expr *xs)
{
size_t n;
pure_expr **elems;
if (pure_is_listv(xs, &n, &elems)) {
- pure_expr *ret = pure_matrix_rowsv(n, elems);
+ pure_expr *ret = matrix_rowsv(n, elems);
if (elems) free(elems);
return ret;
} else
@@ -4587,7 +4801,7 @@
size_t n;
pure_expr **elems;
if (pure_is_listv(xs, &n, &elems)) {
- pure_expr *ret = pure_matrix_columnsv(n, elems);
+ pure_expr *ret = matrix_columnsv(n, elems);
if (elems) free(elems);
return ret;
} else
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 05:31:30 UTC (rev 832)
+++ pure/trunk/runtime.h 2008-09-23 05:54:38 UTC (rev 833)
@@ -703,9 +703,9 @@
pure_expr *matrix_subdiag(pure_expr *x, int32_t k);
pure_expr *matrix_supdiag(pure_expr *x, int32_t k);
-/* Matrix construction. These work like the pure_matrix_rows/
- pure_matrix_columns functions in the public API, but take their input from
- a Pure list. */
+/* Matrix construction. These work like the corresponding functions in the
+ public API, but take their input from a Pure list and raise the appropriate
+ exception in case of dimension mismatch. */
pure_expr *matrix_rows(pure_expr *xs);
pure_expr *matrix_columns(pure_expr *xs);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 09:15:38
|
Revision: 835
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=835&view=rev
Author: agraef
Date: 2008-09-23 09:15:34 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Bugfixes, make list operations work on matrices.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/lib/primitives.pure
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-23 06:32:29 UTC (rev 834)
+++ pure/trunk/ChangeLog 2008-09-23 09:15:34 UTC (rev 835)
@@ -1,3 +1,9 @@
+2008-09-23 Albert Graef <Dr....@t-...>
+
+ * lib/primitives.pure: Added a bunch of new matrix operations. In
+ particular, list operations like filter and map now work on
+ matrices, too.
+
2008-09-20 Albert Graef <Dr....@t-...>
* Implemented basic GSL matrix support, including support for
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-23 06:32:29 UTC (rev 834)
+++ pure/trunk/lib/primitives.pure 2008-09-23 09:15:34 UTC (rev 835)
@@ -416,6 +416,7 @@
extern int matrix_size(expr *x), int matrix_stride(expr *x);
extern expr* matrix_dim(expr *x);
+null x::matrix = #x==0;
#x::matrix = matrix_size x;
dim x::matrix = matrix_dim x;
stride x::matrix = matrix_stride x;
@@ -447,7 +448,7 @@
nth x n = catch (cst {}) (row x n);
mth x m = catch (cst {}) (col x m);
end;
-x::matrix!!ns = if packed x then rowvector x!!(1,ns)
+x::matrix!!ns = if packed x then rowvector x!!([0],ns)
else colcatmap (nth x) ns with
nth x n = catch (cst {}) {x!n};
end;
@@ -469,16 +470,14 @@
cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
-/* Convert a matrix to a list and vice versa. If x is a row vector then list x
- converts it to a list of its elements; otherwise the result is the list of
- the rows of the matrix. You can also use list2 to convert a matrix to a
- list of lists. Conversely, matrix xs converts a list of lists or matrices
- to the corresponding matrix; otherwise, the result is a row vector.
- NOTE: The matrix function may throw a 'bad_matrix_value x' in case of
- dimension mismatch, where x denotes the offending submatrix. */
+/* Convert a matrix to a list and vice versa. list x converts a matrix x to a
+ flat list of its elements, while list2 converts it to a list of lists.
+ Conversely, matrix xs converts a list of lists or matrices specifying the
+ rows of the matrix to the corresponding matrix; otherwise, the result is a
+ row vector. NOTE: The matrix function may throw a 'bad_matrix_value x' in
+ case of dimension mismatch, where x denotes the offending submatrix. */
-list x::matrix = [x!i|i=0..#x-1] if rowvectorp x;
- = rows x otherwise;
+list x::matrix = [x!i|i=0..#x-1];
list2 x::matrix = [[x!(i,j)|j=0..m-1]|i=0..n-1]
when n::int,m::int = dim x end;
@@ -512,12 +511,17 @@
extern expr* matrix_redim(expr* x, int n, int m);
redim x::matrix (n::int,m::int)
- = matrix_redim x n m if n*m==#x;
+ = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
+/* You can also redim a matrix to a given positive row size. In this case the
+ row size must divide the total size of the matrix, */
+
+redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
+
/* Convenience functions for converting a matrix to a row or column vector. */
-rowvector x::matrix = redim x (1,#x);
-colvector x::matrix = redim x (#x,1);
+rowvector x::matrix = redim x (#x);
+colvector x::matrix = redim x 1;
/* Construct matrices from lists of rows and columns. These take either
scalars or submatrices as inputs; corresponding dimensions must match.
@@ -559,6 +563,78 @@
colrev x::matrix = colcat (reverse (cols x));
reverse x::matrix = rowrev (colrev x);
+/* catmap et al on matrices. This allows list and matrix comprehensions to
+ draw values from matrices just as well as from lists, treating the matrix
+ as a flat list of its elements. */
+
+catmap f x::matrix = catmap f (list x);
+rowcatmap f x::matrix = rowcat (map f (list x));
+colcatmap f x::matrix = colcat (map f (list x));
+
+/* Implementations of the other customary list operations, so that these can
+ be used on matrices, too. These operations treat the matrix essentially as
+ if it was a flat list of its elements. Aggregate results are then converted
+ back to matrices with the appropriate dimensions. (This depends on the
+ operation; operations like map and zip keep the dimensions of the input
+ matrix intact, while other functions like filter, take or scanl always
+ return a flat row vector. Also note that the zip-style operations require
+ that the row sizes of all arguments match.) */
+
+cycle x::matrix = cycle (list x);
+cyclen n::int x::matrix = cyclen n (list x) if not null x;
+
+all p x::matrix = all p (list x);
+any p x::matrix = any p (list x);
+do f x::matrix = do f (list x);
+drop k::int x::matrix = x!!(k..#x-1);
+dropwhile p x::matrix = colcat (dropwhile p (list x));
+filter p x::matrix = colcat (filter p (list x));
+foldl f a x::matrix = foldl f a (list x);
+foldl1 f x::matrix = foldl1 f (list x);
+foldr f a x::matrix = foldr f a (list x);
+foldr1 f x::matrix = foldr1 f (list x);
+head x::matrix = x!0 if not null x;
+init x::matrix = x!!(0..#x-2) if not null x;
+last x::matrix = x!(#x-1) if not null x;
+map f x::matrix = flip redim (dim x) $ colcat (map f (list x));
+scanl f a x::matrix = colcat (scanl f a (list x));
+scanl1 f x::matrix = colcat (scanl1 f (list x));
+scanr f a x::matrix = colcat (scanr f a (list x));
+scanr1 f x::matrix = colcat (scanr1 f (list x));
+take k::int x::matrix = x!!(0..k-1);
+takewhile p x::matrix = colcat (takewhile p (list x));
+tail x::matrix = x!!(1..#x-1) if not null x;
+
+zip x::matrix y::matrix = flip redim (dim x!1) $
+ colcat (zip (list x) (list y))
+ if dim x!1==dim y!1;
+zip3 x::matrix y::matrix z::matrix
+ = flip redim (dim x!1) $
+ colcat (zip3 (list x) (list y) (list z))
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+zipwith f x::matrix y::matrix
+ = flip redim (dim x!1) $
+ colcat (zipwith f (list x) (list y))
+ if dim x!1==dim y!1;
+zipwith3 f x::matrix y::matrix z::matrix
+ = flip redim (dim x!1) $
+ colcat (zipwith3 f (list x) (list y) (list z))
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+dowith f x::matrix y::matrix
+ = dowith f (list x) (list y)
+ if dim x!1==dim y!1;
+dowith3 f x::matrix y::matrix z::matrix
+ = dowith3 f (list x) (list y) (list z)
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+
+unzip x::matrix = flip redim (dim x) (colcat u),
+ flip redim (dim x) (colcat v)
+ when u,v = unzip (list x) end;
+unzip3 x::matrix = flip redim (dim x) (colcat u),
+ flip redim (dim x) (colcat v),
+ flip redim (dim x) (colcat w)
+ when u,v,w = unzip3 (list x) end;
+
/* Matrix conversions. */
private matrix_double matrix_complex matrix_int;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 10:32:03
|
Revision: 836
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=836&view=rev
Author: agraef
Date: 2008-09-23 10:31:58 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Moved the matrix operations from primitives.pure to their own module.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/lib/prelude.pure
pure/trunk/lib/primitives.pure
Added Paths:
-----------
pure/trunk/lib/matrices.pure
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-23 09:15:34 UTC (rev 835)
+++ pure/trunk/ChangeLog 2008-09-23 10:31:58 UTC (rev 836)
@@ -1,5 +1,8 @@
2008-09-23 Albert Graef <Dr....@t-...>
+ * lib/matrices.pure: Moved the matrix operations from
+ primitives.pure to their own module.
+
* lib/primitives.pure: Added a bunch of new matrix operations. In
particular, list operations like filter and map now work on
matrices, too.
Added: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure (rev 0)
+++ pure/trunk/lib/matrices.pure 2008-09-23 10:31:58 UTC (rev 836)
@@ -0,0 +1,339 @@
+
+/* Basic matrix functions. */
+
+/* Copyright (c) 2008 by Albert Graef <Dr....@t-...>.
+
+ This file is part of the Pure programming language and system.
+
+ Pure is free software: you can redistribute it and/or modify it under the
+ terms of the GNU General Public License as published by the Free Software
+ Foundation, either version 3 of the License, or (at your option) any later
+ version.
+
+ Pure is distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+/* Additional matrix predicates. The numeric matrix types will only be
+ available if the interpreter was built with support for the GNU Scientific
+ Library (GSL). This is the case iff the global gsl_version variable is
+ predefined by the interpreter. */
+
+extern int matrix_type(expr* x);
+
+dmatrixp x = case x of _::matrix = matrix_type x==1; _ = 0 end;
+cmatrixp x = case x of _::matrix = matrix_type x==2; _ = 0 end;
+imatrixp x = case x of _::matrix = matrix_type x==3; _ = 0 end;
+
+/* The nmatrixp predicate checks for any kind of numeric (double, complex or
+ int) matrix, smatrix for symbolic matrices. */
+
+nmatrixp x = case x of _::matrix = matrix_type x>=1; _ = 0 end;
+smatrixp x = case x of _::matrix = matrix_type x==0; _ = 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;
+
+/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
+ of rows, m the number of columns). Note that Pure supports empty matrices,
+ thus the total size may well be zero, meaning that either the number of
+ rows or the number of columns is zero, or both. The null predicate checks
+ for empty matrices. */
+
+private matrix_size matrix_dim;
+extern int matrix_size(expr *x), expr* matrix_dim(expr *x);
+
+#x::matrix = matrix_size x;
+dim x::matrix = matrix_dim x;
+null x::matrix = #x==0;
+
+/* The stride of a matrix denotes the real row size of the underlying C array,
+ see the description of the 'pack' function below for further details.
+ There's little use for this value in Pure, but it may be needed when
+ interfacing to C. */
+
+private matrix_stride;
+extern int matrix_stride(expr *x);
+
+stride x::matrix = matrix_stride x;
+
+/* Indexing. Note that in difference to Octave and MATLAB, all indices are
+ zero-based, and you *must* use Pure's indexing operator '!' to retrieve an
+ element. You can either get an element by specifying its row-major index in
+ the range 0..#x-1, or with a two-dimensional subscript of the form of a
+ pair (i,j), where i denotes the row and j the column index. Both operations
+ take only constant time, and an 'out_of_bounds' exception is thrown if an
+ index falls outside the valid range. */
+
+private matrix_elem_at matrix_elem_at2;
+extern expr* matrix_elem_at(expr* x, int i);
+extern expr* matrix_elem_at2(expr* x, int i, int j);
+
+x::matrix!i::int = matrix_elem_at x i if i>=0 && i<#x;
+ = throw out_of_bounds otherwise;
+
+x::matrix!(i::int,j::int)
+ = matrix_elem_at2 x i j
+ if (i>=0 && i<n && j>=0 && j<m
+ when n::int,m::int = dim x end);
+ = throw out_of_bounds otherwise;
+
+/* Matrix slices (x!!ns). As with simple indexing, elements can addressed
+ using either singleton (row-major) indices or index pair (row,column). The
+ former is specified as a list of int values, the latter as a pair of lists
+ of int values. As with list slicing, index ranges may be non-contiguous
+ and/or non-monotonous. However, the case of contiguous and monotonous
+ ranges is optimized by making good use of the 'submat' operation below. */
+
+x::matrix!!(ns,ms) = case ns,ms of
+ ns@(n:_),ms@(m:_) = submat x (n,m) (#ns,#ms)
+ if cont ns && cont ms;
+ _ = colcatmap (mth (rowcatmap (nth x) ns)) ms;
+ end with
+ cont [n] = 1;
+ cont (n::int:ns@(m::int:_)) = cont ns if m==n+1;
+ cont _ = 0 otherwise;
+ nth x n = catch (cst {}) (row x n);
+ mth x m = catch (cst {}) (col x m);
+ end;
+x::matrix!!ns = if packed x then rowvector x!!([0],ns)
+ else colcatmap (nth x) ns with
+ nth x n = catch (cst {}) {x!n};
+ end;
+
+/* 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;
+
+/* Convert a matrix to a list and vice versa. list x converts a matrix x to a
+ flat list of its elements, while list2 converts it to a list of lists.
+ Conversely, matrix xs converts a list of lists or matrices specifying the
+ rows of the matrix to the corresponding matrix; otherwise, the result is a
+ row vector. NOTE: The matrix function may throw a 'bad_matrix_value x' in
+ case of dimension mismatch, where x denotes the offending submatrix. */
+
+list x::matrix = [x!i|i=0..#x-1];
+list2 x::matrix = [[x!(i,j)|j=0..m-1]|i=0..n-1]
+ when n::int,m::int = dim x end;
+
+matrix [] = {};
+matrix xs@(x:_) = rowcatmap colcat xs if all listp xs;
+ = rowcat xs if any matrixp xs;
+ = colcat xs otherwise;
+
+/* Extract (sub-,super-) diagonals from a matrix. Sub- and super-diagonals for
+ k=0 return the main diagonal. Indices for sub- and super-diagonals can also
+ be negative, in which case the corresponding super- or sub-diagonal is
+ returned instead. In each case the result is a row vector. */
+
+private matrix_diag matrix_subdiag matrix_supdiag;
+extern expr* matrix_diag(expr* x) = diag;
+extern expr* matrix_subdiag(expr* x, int k) = subdiag;
+extern expr* matrix_supdiag(expr* x, int k) = supdiag;
+
+/* Extract a submatrix of a given size at a given offset. The result shares
+ the underlying storage with the input matrix (i.e., matrix elements are
+ *not* copied) and so this is a comparatively cheap operation. */
+
+submat x::matrix (i::int,j::int) (n::int,m::int)
+ = matrix_slice x i j (i+n-1) (j+m-1);
+
+/* Change the dimensions of a matrix without changing its size. The total
+ number of elements must match that of the input matrix. Reuses the
+ underlying storage of the input matrix if possible. */
+
+private matrix_redim;
+extern expr* matrix_redim(expr* x, int n, int m);
+
+redim x::matrix (n::int,m::int)
+ = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
+
+/* You can also redim a matrix to a given row size. In this case the row size
+ must be positive and divide the total size of the matrix, */
+
+redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
+
+/* Convenience functions for converting a matrix to a row or column vector. */
+
+rowvector x::matrix = redim x (#x);
+colvector x::matrix = redim x 1;
+
+/* 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}. NOTE: Like the built-in matrix constructs, these
+ operations may throw a 'bad_matrix_value x' in case of dimension mismatch,
+ where x denotes the offending submatrix. */
+
+extern expr* matrix_rows(expr *x) = rowcat;
+extern expr* matrix_columns(expr *x) = colcat;
+
+/* Combinations of rowcat/colcat and map. */
+
+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;
+
+/* Reverse a matrix. 'rowrev' reverses the rows, 'colrev' the columns,
+ 'reverse' both dimensions. */
+
+rowrev x::matrix = rowcat (reverse (rows x));
+colrev x::matrix = colcat (reverse (cols x));
+reverse x::matrix = rowrev (colrev x);
+
+/* catmap et al on matrices. This allows list and matrix comprehensions to
+ draw values from matrices as well as from lists, treating the matrix as a
+ flat list of its elements. */
+
+catmap f x::matrix = catmap f (list x);
+rowcatmap f x::matrix = rowcat (map f (list x));
+colcatmap f x::matrix = colcat (map f (list x));
+
+/* Implementations of the other customary list operations, so that these can
+ be used on matrices, too. These operations treat the matrix essentially as
+ if it was a flat list of its elements. Aggregate results are then converted
+ back to matrices with the appropriate dimensions. (This depends on the
+ particular operation; functions like map and zip keep the dimensions of the
+ input matrix intact, while other functions like filter, take or scanl
+ always return a flat row vector. Also note that the zip-style operations
+ require that the row sizes of all arguments match.) */
+
+cycle x::matrix = cycle (list x);
+cyclen n::int x::matrix = cyclen n (list x) if not null x;
+
+all p x::matrix = all p (list x);
+any p x::matrix = any p (list x);
+do f x::matrix = do f (list x);
+drop k::int x::matrix = x!!(k..#x-1);
+dropwhile p x::matrix = colcat (dropwhile p (list x));
+filter p x::matrix = colcat (filter p (list x));
+foldl f a x::matrix = foldl f a (list x);
+foldl1 f x::matrix = foldl1 f (list x);
+foldr f a x::matrix = foldr f a (list x);
+foldr1 f x::matrix = foldr1 f (list x);
+head x::matrix = x!0 if not null x;
+init x::matrix = x!!(0..#x-2) if not null x;
+last x::matrix = x!(#x-1) if not null x;
+map f x::matrix = flip redim (dim x) $ colcat (map f (list x));
+scanl f a x::matrix = colcat (scanl f a (list x));
+scanl1 f x::matrix = colcat (scanl1 f (list x));
+scanr f a x::matrix = colcat (scanr f a (list x));
+scanr1 f x::matrix = colcat (scanr1 f (list x));
+take k::int x::matrix = x!!(0..k-1);
+takewhile p x::matrix = colcat (takewhile p (list x));
+tail x::matrix = x!!(1..#x-1) if not null x;
+
+zip x::matrix y::matrix = flip redim (dim x!1) $
+ colcat (zip (list x) (list y))
+ if dim x!1==dim y!1;
+zip3 x::matrix y::matrix z::matrix
+ = flip redim (dim x!1) $
+ colcat (zip3 (list x) (list y) (list z))
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+zipwith f x::matrix y::matrix
+ = flip redim (dim x!1) $
+ colcat (zipwith f (list x) (list y))
+ if dim x!1==dim y!1;
+zipwith3 f x::matrix y::matrix z::matrix
+ = flip redim (dim x!1) $
+ colcat (zipwith3 f (list x) (list y) (list z))
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+dowith f x::matrix y::matrix
+ = dowith f (list x) (list y)
+ if dim x!1==dim y!1;
+dowith3 f x::matrix y::matrix z::matrix
+ = dowith3 f (list x) (list y) (list z)
+ if dim x!1==dim y!1 && dim x!1==dim z!1;
+
+unzip x::matrix = flip redim (dim x) (colcat u),
+ flip redim (dim x) (colcat v)
+ when u,v = unzip (list x) end;
+unzip3 x::matrix = flip redim (dim x) (colcat u),
+ flip redim (dim x) (colcat v),
+ flip redim (dim x) (colcat w)
+ when u,v,w = unzip3 (list x) end;
+
+/* Matrix conversions. These convert between different types of numeric
+ matrices. You can also extract the real and imaginary parts of a (complex)
+ matrix. */
+
+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 = matrix_double x if imatrixp x || dmatrixp x;
+imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
+cmatrix x::matrix = matrix_complex x if nmatrixp x;
+
+private matrix_re matrix_im;
+extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
+
+re x::matrix = matrix_re x if nmatrixp x;
+im x::matrix = matrix_im x if nmatrixp x;
+
+/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
+ contiguous storage. If possible, it also frees up extra memory if the
+ matrix was created as a slice from a bigger matrix. The 'packed' predicate
+ can be used to verify whether a matrix is already packed. */
+
+pack x::matrix = colcat [x,{}];
+packed x::matrix = stride x==dim x!1;
+
+/* Convert a matrix to a pointer. This returns a pointer to the underlying C
+ array, which allows you to change the data in-place by shelling out to C
+ functions, so beware. You also need to be careful when passing such a
+ pointer to C functions if the underlying data is non-contiguous; when in
+ doubt use the 'pack' function from above to place the data in contiguous
+ storage first. */
+
+pointer x::matrix = pure_pointerval x;
+
+/* Matrix comparisons. */
+
+x::matrix == y::matrix = x === y
+ if nmatrixp x && matrix_type x == matrix_type y;
+// mixed numeric cases
+ = cmatrix x === y if nmatrixp x && cmatrixp y;
+ = x === cmatrix y if cmatrixp x && nmatrixp y;
+ = dmatrix x === y if imatrixp x && dmatrixp y;
+ = x === dmatrix y if dmatrixp x && imatrixp y;
+// comparisons with symbolic matrices
+ = 0 if dim x != dim y;
+ = compare 0 with
+ compare i::int = 1 if i>=n;
+ = 0 if x!i != y!i;
+ = compare (i+1);
+ end when n::int = #x end;
+
+x::matrix != y::matrix = not x == y;
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-23 09:15:34 UTC (rev 835)
+++ pure/trunk/lib/prelude.pure 2008-09-23 10:31:58 UTC (rev 836)
@@ -74,7 +74,7 @@
Note that the math and system modules are *not* included here, so you have
to do that yourself if your program requires any of those operations. */
-using primitives, strings;
+using primitives, matrices, strings;
/* Basic combinators. */
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-23 09:15:34 UTC (rev 835)
+++ pure/trunk/lib/primitives.pure 2008-09-23 10:31:58 UTC (rev 836)
@@ -41,33 +41,8 @@
doublep x = case x of _::double = 1; _ = 0 end;
stringp x = case x of _::string = 1; _ = 0 end;
pointerp x = case x of _::pointer = 1; _ = 0 end;
-
-/* Matrix predicates. The numeric matrix types will only be available if the
- interpreter was built with support for the GNU Scientific Library (GSL).
- This is the case iff the global gsl_version variable is predefined by the
- interpreter. */
-
-extern int matrix_type(expr* x);
-
matrixp x = case x of _::matrix = 1; _ = 0 end;
-dmatrixp x = case x of _::matrix = matrix_type x==1; _ = 0 end;
-cmatrixp x = case x of _::matrix = matrix_type x==2; _ = 0 end;
-imatrixp x = case x of _::matrix = matrix_type x==3; _ = 0 end;
-/* The nmatrixp predicate checks for any kind of numeric (double, complex or
- int) matrix, smatrix for symbolic matrices. */
-
-nmatrixp x = case x of _::matrix = matrix_type x>=1; _ = 0 end;
-smatrixp x = case x of _::matrix = matrix_type x==0; _ = 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;
-
/* Predicates to check for function objects, global (unbound) variables,
function applications, proper lists, list nodes and tuples. */
@@ -117,8 +92,7 @@
pointer x::int |
pointer x::bigint |
pointer x::double |
-pointer x::string |
-pointer x::matrix = pure_pointerval x;
+pointer x::string = pure_pointerval x;
/* Convert signed (8/16/32/64) bit integers to the corresponding unsigned
quantities. These functions behave as if the value was "cast" to the
@@ -410,266 +384,6 @@
x::pointer==y::pointer = bigint x == bigint y;
x::pointer!=y::pointer = bigint x != bigint y;
-/* Basic matrix operations: size, dimensions and indexing. */
-
-private matrix_size matrix_dim matrix_stride;
-extern int matrix_size(expr *x), int matrix_stride(expr *x);
-extern expr* matrix_dim(expr *x);
-
-null x::matrix = #x==0;
-#x::matrix = matrix_size x;
-dim x::matrix = matrix_dim x;
-stride x::matrix = matrix_stride x;
-
-private matrix_elem_at matrix_elem_at2;
-extern expr* matrix_elem_at(expr* x, int i);
-extern expr* matrix_elem_at2(expr* x, int i, int j);
-
-x::matrix!i::int = matrix_elem_at x i if i>=0 && i<#x;
- = throw out_of_bounds otherwise;
-
-x::matrix!(i::int,j::int)
- = matrix_elem_at2 x i j
- if (i>=0 && i<n && j>=0 && j<m
- when n::int,m::int = dim x end);
- = throw out_of_bounds otherwise;
-
-/* Matrix slices. */
-
-x::matrix!!(ns,ms) = case ns,ms of
- // optimize the case of contiguous slices
- ns@(n:_),ms@(m:_) = submat x (n,m) (#ns,#ms)
- if cont ns && cont ms;
- _ = colcatmap (mth (rowcatmap (nth x) ns)) ms;
- end with
- cont [n] = 1;
- cont (n::int:ns@(m::int:_)) = cont ns if m==n+1;
- cont _ = 0 otherwise;
- nth x n = catch (cst {}) (row x n);
- mth x m = catch (cst {}) (col x m);
- end;
-x::matrix!!ns = if packed x then rowvector x!!([0],ns)
- else colcatmap (nth x) ns with
- nth x n = catch (cst {}) {x!n};
- end;
-
-/* 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;
-
-/* Convert a matrix to a list and vice versa. list x converts a matrix x to a
- flat list of its elements, while list2 converts it to a list of lists.
- Conversely, matrix xs converts a list of lists or matrices specifying the
- rows of the matrix to the corresponding matrix; otherwise, the result is a
- row vector. NOTE: The matrix function may throw a 'bad_matrix_value x' in
- case of dimension mismatch, where x denotes the offending submatrix. */
-
-list x::matrix = [x!i|i=0..#x-1];
-list2 x::matrix = [[x!(i,j)|j=0..m-1]|i=0..n-1]
- when n::int,m::int = dim x end;
-
-matrix [] = {};
-matrix xs@(x:_) = rowcatmap colcat xs if all listp xs;
- = rowcat xs if any matrixp xs;
- = colcat xs otherwise;
-
-/* Extract (sub-,super-) diagonals from a matrix. Sub- and super-diagonals for
- k=0 return the main diagonal. Indices for sub- and super-diagonals can also
- be negative, in which case the corresponding super- or sub-diagonal is
- returned instead. In each case the result is a row vector. */
-
-private matrix_diag matrix_subdiag matrix_supdiag;
-extern expr* matrix_diag(expr* x) = diag;
-extern expr* matrix_subdiag(expr* x, int k) = subdiag;
-extern expr* matrix_supdiag(expr* x, int k) = supdiag;
-
-/* Extract a submatrix of a given size at a given offset. The result shares
- the underlying storage with the input matrix (i.e., matrix elements are
- *not* copied) and so this is a comparatively cheap operation. */
-
-submat x::matrix (i::int,j::int) (n::int,m::int)
- = matrix_slice x i j (i+n-1) (j+m-1);
-
-/* Change the dimensions of a matrix without changing its size. The total
- number of elements must match that of the input matrix. Reuses the
- underlying storage of the input matrix if possible. */
-
-private matrix_redim;
-extern expr* matrix_redim(expr* x, int n, int m);
-
-redim x::matrix (n::int,m::int)
- = matrix_redim x n m if n>=0 && m>=0 && n*m==#x;
-
-/* You can also redim a matrix to a given positive row size. In this case the
- row size must divide the total size of the matrix, */
-
-redim x::matrix m::int = redim x (#x div m,m) if m>0 && #x mod m==0;
-
-/* Convenience functions for converting a matrix to a row or column vector. */
-
-rowvector x::matrix = redim x (#x);
-colvector x::matrix = redim x 1;
-
-/* 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}. NOTE: Like the built-in matrix constructs, these
- operations may throw a 'bad_matrix_value x' in case of dimension mismatch,
- where x denotes the offending submatrix. */
-
-extern expr* matrix_rows(expr *x) = rowcat;
-extern expr* matrix_columns(expr *x) = colcat;
-
-/* Combinations of rowcat/colcat and map. */
-
-rowcatmap f [] = {};
-rowcatmap f xs@(_:_) = rowcat (map f xs);
-
-colcatmap f [] = {};
-colcatmap f xs@(_:_) = colcat (map f xs);
-
-/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
- contiguous storage. If possible, it also frees up extra memory if the
- matrix was created as a slice from a bigger matrix. The 'packed' predicate
- can be used to verify whether a matrix is already packed. */
-
-pack x::matrix = colcat [x,{}];
-packed x::matrix = stride x==dim x!1;
-
-/* Transpose a matrix. */
-
-private matrix_transpose;
-extern expr* matrix_transpose(expr *x);
-
-x::matrix' = matrix_transpose x;
-
-/* Reverse a matrix. rowrev reverses the rows, colrev the columns, reverse
- both dimensions. */
-
-rowrev x::matrix = rowcat (reverse (rows x));
-colrev x::matrix = colcat (reverse (cols x));
-reverse x::matrix = rowrev (colrev x);
-
-/* catmap et al on matrices. This allows list and matrix comprehensions to
- draw values from matrices just as well as from lists, treating the matrix
- as a flat list of its elements. */
-
-catmap f x::matrix = catmap f (list x);
-rowcatmap f x::matrix = rowcat (map f (list x));
-colcatmap f x::matrix = colcat (map f (list x));
-
-/* Implementations of the other customary list operations, so that these can
- be used on matrices, too. These operations treat the matrix essentially as
- if it was a flat list of its elements. Aggregate results are then converted
- back to matrices with the appropriate dimensions. (This depends on the
- operation; operations like map and zip keep the dimensions of the input
- matrix intact, while other functions like filter, take or scanl always
- return a flat row vector. Also note that the zip-style operations require
- that the row sizes of all arguments match.) */
-
-cycle x::matrix = cycle (list x);
-cyclen n::int x::matrix = cyclen n (list x) if not null x;
-
-all p x::matrix = all p (list x);
-any p x::matrix = any p (list x);
-do f x::matrix = do f (list x);
-drop k::int x::matrix = x!!(k..#x-1);
-dropwhile p x::matrix = colcat (dropwhile p (list x));
-filter p x::matrix = colcat (filter p (list x));
-foldl f a x::matrix = foldl f a (list x);
-foldl1 f x::matrix = foldl1 f (list x);
-foldr f a x::matrix = foldr f a (list x);
-foldr1 f x::matrix = foldr1 f (list x);
-head x::matrix = x!0 if not null x;
-init x::matrix = x!!(0..#x-2) if not null x;
-last x::matrix = x!(#x-1) if not null x;
-map f x::matrix = flip redim (dim x) $ colcat (map f (list x));
-scanl f a x::matrix = colcat (scanl f a (list x));
-scanl1 f x::matrix = colcat (scanl1 f (list x));
-scanr f a x::matrix = colcat (scanr f a (list x));
-scanr1 f x::matrix = colcat (scanr1 f (list x));
-take k::int x::matrix = x!!(0..k-1);
-takewhile p x::matrix = colcat (takewhile p (list x));
-tail x::matrix = x!!(1..#x-1) if not null x;
-
-zip x::matrix y::matrix = flip redim (dim x!1) $
- colcat (zip (list x) (list y))
- if dim x!1==dim y!1;
-zip3 x::matrix y::matrix z::matrix
- = flip redim (dim x!1) $
- colcat (zip3 (list x) (list y) (list z))
- if dim x!1==dim y!1 && dim x!1==dim z!1;
-zipwith f x::matrix y::matrix
- = flip redim (dim x!1) $
- colcat (zipwith f (list x) (list y))
- if dim x!1==dim y!1;
-zipwith3 f x::matrix y::matrix z::matrix
- = flip redim (dim x!1) $
- colcat (zipwith3 f (list x) (list y) (list z))
- if dim x!1==dim y!1 && dim x!1==dim z!1;
-dowith f x::matrix y::matrix
- = dowith f (list x) (list y)
- if dim x!1==dim y!1;
-dowith3 f x::matrix y::matrix z::matrix
- = dowith3 f (list x) (list y) (list z)
- if dim x!1==dim y!1 && dim x!1==dim z!1;
-
-unzip x::matrix = flip redim (dim x) (colcat u),
- flip redim (dim x) (colcat v)
- when u,v = unzip (list x) end;
-unzip3 x::matrix = flip redim (dim x) (colcat u),
- flip redim (dim x) (colcat v),
- flip redim (dim x) (colcat w)
- when u,v,w = unzip3 (list x) end;
-
-/* 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);
-
-dmatrix x::matrix = matrix_double x if imatrixp x || dmatrixp x;
-imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
-cmatrix x::matrix = matrix_complex x if nmatrixp x;
-
-private matrix_re matrix_im;
-extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
-
-re x::matrix = matrix_re x if nmatrixp x;
-im x::matrix = matrix_im x if nmatrixp x;
-
-/* Matrix comparisons. */
-
-x::matrix == y::matrix = x === y
- if nmatrixp x && matrix_type x == matrix_type y;
-// mixed numeric cases
- = cmatrix x === y if nmatrixp x && cmatrixp y;
- = x === cmatrix y if cmatrixp x && nmatrixp y;
- = dmatrix x === y if imatrixp x && dmatrixp y;
- = x === dmatrix y if dmatrixp x && imatrixp y;
-// comparisons with symbolic matrices
- = 0 if dim x != dim y;
- = compare 0 with
- compare i::int = 1 if i>=n;
- = 0 if x!i != y!i;
- = compare (i+1);
- end when n::int = #x end;
-
-x::matrix != y::matrix = not x == y;
-
/* IEEE floating point infinities and NaNs. Place these after the definitions
of the built-in operators so that the double arithmetic works. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 12:56:57
|
Revision: 839
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=839&view=rev
Author: agraef
Date: 2008-09-23 12:56:30 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add more matrix creation and pointer->matrix conversion functions.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 10:58:30 UTC (rev 838)
+++ pure/trunk/lib/matrices.pure 2008-09-23 12:56:30 UTC (rev 839)
@@ -43,6 +43,18 @@
rowvectorp x = matrixp x && dim x!0==1;
colvectorp x = matrixp x && dim x!1==1;
+/* Convenience functions to create numeric matrices with the given dimensions
+ (either a pair denoting the number of rows and columns, or just the row
+ size in order to create a row vector), and all elements zero. */
+
+dmatrix (n::int,m::int) = cdmatrix (pointer 0) (n,m);
+cmatrix (n::int,m::int) = ccmatrix (pointer 0) (n,m);
+imatrix (n::int,m::int) = cimatrix (pointer 0) (n,m);
+
+dmatrix n::int = cdmatrix (pointer 0) n;
+cmatrix n::int = ccmatrix (pointer 0) n;
+imatrix n::int = cimatrix (pointer 0) n;
+
/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
of rows, m the number of columns). Note that Pure supports empty matrices,
thus the total size may well be zero, meaning that either the number of
@@ -332,6 +344,37 @@
pointer x::matrix = pure_pointerval x;
+/* Create a numeric matrix from a pointer. The pointer to be converted must
+ point to a malloc'ed, properly initialized memory area big enough to
+ accommodate the requested dimensions. The pointer may also be NULL in which
+ case a suitable pointer is allocated automatically. This memory is taken
+ over by Pure and will be freed automatically when the matrix object is
+ garbage-collected. CAVEAT: These are low-level operations and hence should
+ be used with care. They are useful, in particular, if you need to handle
+ massive amounts of matrix or vector data generated or processed by external
+ C software, such as graphics and sound libraries. */
+
+private matrix_from_double_array;
+private matrix_from_complex_array;
+private matrix_from_int_array;
+extern expr* matrix_from_double_array(void* p, int n, int m);
+extern expr* matrix_from_complex_array(void* p, int n, int m);
+extern expr* matrix_from_int_array(void* p, int n, int m);
+
+cdmatrix p::pointer (n::int,m::int)
+ = matrix_from_double_array p n m;
+ccmatrix p::pointer (n::int,m::int)
+ = matrix_from_complex_array p n m;
+cimatrix p::pointer (n::int,m::int)
+ = matrix_from_int_array p n m;
+
+cdmatrix p::pointer n::int
+ = cdmatrix p (1,n);
+ccmatrix p::pointer n::int
+ = ccmatrix p (1,n);
+cimatrix p::pointer n::int
+ = cimatrix p (1,n);
+
/* Matrix comparisons. */
x::matrix == y::matrix = x === y
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 10:58:30 UTC (rev 838)
+++ pure/trunk/runtime.cc 2008-09-23 12:56:30 UTC (rev 839)
@@ -105,7 +105,7 @@
s.tda = m->tda;
s.block = m->block;
s.owner = 0;
- view.matrix = s;
+ view.matrix = s;
return view;
}
}
@@ -5010,6 +5010,95 @@
#endif
}
+extern "C"
+pure_expr *matrix_from_double_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ if (!p) p = calloc(n1*n2, sizeof(double));
+ if (!p) return 0;
+ gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::DMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ if (!p) p = calloc(2*n1*n2, sizeof(double));
+ if (!p) return 0;
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ gsl_block_complex *b = (gsl_block_complex*)malloc(sizeof(gsl_block_complex));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p) p = calloc(n1*n2, sizeof(int));
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#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-23 10:58:30 UTC (rev 838)
+++ pure/trunk/runtime.h 2008-09-23 12:56:30 UTC (rev 839)
@@ -733,6 +733,18 @@
pure_expr *matrix_re(pure_expr *x);
pure_expr *matrix_im(pure_expr *x);
+/* Create a matrix object of the given dimensions which uses the given pointer
+ p as its underlying C array. There are no checks whatsoever, so the caller
+ is responsible for making sure that the memory has the right size and is
+ properly initialized. p must point to a malloc'ed memory area, which is
+ taken over by Pure and will be freed automatically when the matrix is
+ garbage-collected. p may also be NULL in which case a suitable pointer is
+ allocated automatically. */
+
+pure_expr *matrix_from_double_array(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_complex_array(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_int_array(void *p, uint32_t n, uint32_t m);
+
/* 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-23 21:56:57
|
Revision: 842
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=842&view=rev
Author: agraef
Date: 2008-09-23 21:56:51 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add some missing memory operations.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-23 18:56:29 UTC (rev 841)
+++ pure/trunk/lib/primitives.pure 2008-09-23 21:56:51 UTC (rev 842)
@@ -396,30 +396,40 @@
/* Direct memory accesses. Use with care ... or else! */
-private pointer_get_byte pointer_get_int pointer_get_double
+private pointer_get_byte pointer_get_short pointer_get_int
+ pointer_get_float pointer_get_double
pointer_get_string pointer_get_pointer;
extern int pointer_get_byte(void *ptr);
+extern int pointer_get_short(void *ptr);
extern int pointer_get_int(void *ptr);
+extern double pointer_get_float(void *ptr);
extern double pointer_get_double(void *ptr);
extern char *pointer_get_string(void *ptr);
extern void *pointer_get_pointer(void *ptr);
get_byte x::pointer = pointer_get_byte x;
+get_short x::pointer = pointer_get_short x;
get_int x::pointer = pointer_get_int x;
+get_float x::pointer = pointer_get_float x;
get_double x::pointer = pointer_get_double x;
get_string x::pointer = pointer_get_string x;
get_pointer x::pointer = pointer_get_pointer x;
-private pointer_put_byte pointer_put_int pointer_put_double
+private pointer_put_byte pointer_put_short pointer_put_int
+ pointer_put_float pointer_put_double
pointer_put_string pointer_put_pointer;
extern void pointer_put_byte(void *ptr, int x); // IMPURE!
+extern void pointer_put_short(void *ptr, int x); // IMPURE!
extern void pointer_put_int(void *ptr, int x); // IMPURE!
+extern void pointer_put_float(void *ptr, double x); // IMPURE!
extern void pointer_put_double(void *ptr, double x); // IMPURE!
extern void pointer_put_string(void *ptr, char *x); // IMPURE!
extern void pointer_put_pointer(void *ptr, void *x); // IMPURE!
put_byte x::pointer y::int = pointer_put_byte x y;
+put_short x::pointer y::int = pointer_put_short x y;
put_int x::pointer y::int = pointer_put_int x y;
+put_float x::pointer y::double = pointer_put_float x y;
put_double x::pointer y::double = pointer_put_double x y;
put_string x::pointer y::string = pointer_put_string x y;
put_pointer x::pointer y::string = pointer_put_pointer x y;
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 18:56:29 UTC (rev 841)
+++ pure/trunk/runtime.cc 2008-09-23 21:56:51 UTC (rev 842)
@@ -5301,11 +5301,18 @@
extern "C"
int32_t pointer_get_byte(void *ptr)
{
- uint8_t *p = (uint8_t*)ptr;
+ int8_t *p = (int8_t*)ptr;
return *p;
}
extern "C"
+int32_t pointer_get_short(void *ptr)
+{
+ int16_t *p = (int16_t*)ptr;
+ return *p;
+}
+
+extern "C"
int32_t pointer_get_int(void *ptr)
{
int32_t *p = (int32_t*)ptr;
@@ -5313,6 +5320,13 @@
}
extern "C"
+double pointer_get_float(void *ptr)
+{
+ float *p = (float*)ptr;
+ return *p;
+}
+
+extern "C"
double pointer_get_double(void *ptr)
{
double *p = (double*)ptr;
@@ -5343,11 +5357,18 @@
extern "C"
void pointer_put_byte(void *ptr, int32_t x)
{
- uint8_t *p = (uint8_t*)ptr;
+ int8_t *p = (int8_t*)ptr;
*p = x;
}
extern "C"
+void pointer_put_short(void *ptr, int32_t x)
+{
+ int16_t *p = (int16_t*)ptr;
+ *p = x;
+}
+
+extern "C"
void pointer_put_int(void *ptr, int32_t x)
{
int32_t *p = (int32_t*)ptr;
@@ -5355,6 +5376,13 @@
}
extern "C"
+void pointer_put_float(void *ptr, double x)
+{
+ float *p = (float*)ptr;
+ *p = x;
+}
+
+extern "C"
void pointer_put_double(void *ptr, double x)
{
double *p = (double*)ptr;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 18:56:29 UTC (rev 841)
+++ pure/trunk/runtime.h 2008-09-23 21:56:51 UTC (rev 842)
@@ -767,14 +767,18 @@
you'll have to use the memory management routines above to do that. */
int32_t pointer_get_byte(void *ptr);
+int32_t pointer_get_short(void *ptr);
int32_t pointer_get_int(void *ptr);
+double pointer_get_float(void *ptr);
double pointer_get_double(void *ptr);
char *pointer_get_string(void *ptr);
void *pointer_get_pointer(void *ptr);
pure_expr *pointer_get_expr(void *ptr);
void pointer_put_byte(void *ptr, int32_t x);
+void pointer_put_short(void *ptr, int32_t x);
void pointer_put_int(void *ptr, int32_t x);
+void pointer_put_float(void *ptr, double x);
void pointer_put_double(void *ptr, double x);
void pointer_put_string(void *ptr, const char *x);
void pointer_put_pointer(void *ptr, void *x);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-23 23:25:32
|
Revision: 844
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=844&view=rev
Author: agraef
Date: 2008-09-23 23:25:21 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Overhaul of matrix-pointer operations.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 22:03:18 UTC (rev 843)
+++ pure/trunk/lib/matrices.pure 2008-09-23 23:25:21 UTC (rev 844)
@@ -43,17 +43,36 @@
rowvectorp x = matrixp x && dim x!0==1;
colvectorp x = matrixp x && dim x!1==1;
-/* Convenience functions to create numeric matrices with the given dimensions
+/* Matrix comparisons. */
+
+x::matrix == y::matrix = x === y
+ if nmatrixp x && matrix_type x == matrix_type y;
+// mixed numeric cases
+ = cmatrix x === y if nmatrixp x && cmatrixp y;
+ = x === cmatrix y if cmatrixp x && nmatrixp y;
+ = dmatrix x === y if imatrixp x && dmatrixp y;
+ = x === dmatrix y if dmatrixp x && imatrixp y;
+// comparisons with symbolic matrices
+ = 0 if dim x != dim y;
+ = compare 0 with
+ compare i::int = 1 if i>=n;
+ = 0 if x!i != y!i;
+ = compare (i+1);
+ end when n::int = #x end;
+
+x::matrix != y::matrix = not x == y;
+
+/* Convenience functions to create zero matrices with the given dimensions
(either a pair denoting the number of rows and columns, or just the row
- size in order to create a row vector), and all elements zero. */
+ size in order to create a row vector). */
-dmatrix (n::int,m::int) = cdmatrix (pointer 0) (n,m);
-cmatrix (n::int,m::int) = ccmatrix (pointer 0) (n,m);
-imatrix (n::int,m::int) = cimatrix (pointer 0) (n,m);
+dmatrix (n::int,m::int) = dmatrix_from_array_dup (pointer 0) (n,m);
+cmatrix (n::int,m::int) = cmatrix_from_array_dup (pointer 0) (n,m);
+imatrix (n::int,m::int) = imatrix_from_array_dup (pointer 0) (n,m);
-dmatrix n::int = cdmatrix (pointer 0) n;
-cmatrix n::int = ccmatrix (pointer 0) n;
-imatrix n::int = cimatrix (pointer 0) n;
+dmatrix n::int = dmatrix (1,n);
+cmatrix n::int = cmatrix (1,n);
+imatrix n::int = imatrix (1,n);
/* Size of a matrix (#x) and its dimensions (dim x=n,m where n is the number
of rows, m the number of columns). Note that Pure supports empty matrices,
@@ -328,34 +347,34 @@
im x::matrix = matrix_im x if nmatrixp x;
/* 'Pack' a matrix. This creates a copy of the matrix which has the data in
- contiguous storage. If possible, it also frees up extra memory if the
- matrix was created as a slice from a bigger matrix. The 'packed' predicate
- can be used to verify whether a matrix is already packed. */
+ contiguous storage. It also frees up extra memory if the matrix was created
+ as a slice from a bigger matrix. The 'packed' predicate can be used to
+ verify whether a matrix is already packed. */
pack x::matrix = colcat [x,{}];
packed x::matrix = stride x==dim x!1;
-/* Convert a matrix to a pointer. This returns a pointer to the underlying C
- array, which allows you to change the data in-place by shelling out to C
- functions, so beware. You also need to be careful when passing such a
- pointer to C functions if the underlying data is non-contiguous; when in
- doubt use the 'pack' function from above to place the data in contiguous
- storage first. */
+/* Low-level operations for converting between matrices and raw pointers.
+ These are typically used to shovel around massive amounts of numeric data
+ between Pure and external C routines, when performance and throughput is an
+ important consideration (e.g., graphics, video and audio applications). The
+ usual caveats apply. */
+/* Convert a matrix to a pointer. Returns a pointer pointing directly to the
+ underlying C array. You have to be careful when passing such a pointer to C
+ functions if the underlying data is non-contiguous; when in doubt, first
+ use the 'pack' function from above to place the data in contiguous
+ storage. */
+
private pure_pointerval;
extern expr* pure_pointerval(expr*);
pointer x::matrix = pure_pointerval x;
-/* Create a numeric matrix from a pointer. The pointer to be converted must
- point to a malloc'ed, properly initialized memory area big enough to
- accommodate the requested dimensions. The pointer may also be NULL in which
- case a suitable pointer is allocated automatically. This memory is taken
- over by Pure and will be freed automatically when the matrix object is
- garbage-collected. CAVEAT: These are low-level operations and hence should
- be used with care. They are useful, in particular, if you need to handle
- massive amounts of matrix or vector data generated or processed by external
- C functions, such as routines dealing with graphics and sound data. */
+/* Create a numeric matrix from a pointer, without copying the data. The
+ caller must ensure that the pointer points to properly initialized memory
+ big enough to accommodate the requested dimensions, which persists for the
+ entire lifetime of the matrix object. */
private matrix_from_double_array;
private matrix_from_complex_array;
@@ -364,35 +383,41 @@
extern expr* matrix_from_complex_array(void* p, int n, int m);
extern expr* matrix_from_int_array(void* p, int n, int m);
-cdmatrix p::pointer (n::int,m::int)
+dmatrix_from_array p::pointer (n::int,m::int)
= matrix_from_double_array p n m;
-ccmatrix p::pointer (n::int,m::int)
+cmatrix_from_array p::pointer (n::int,m::int)
= matrix_from_complex_array p n m;
-cimatrix p::pointer (n::int,m::int)
+imatrix_from_array p::pointer (n::int,m::int)
= matrix_from_int_array p n m;
-cdmatrix p::pointer n::int
- = cdmatrix p (1,n);
-ccmatrix p::pointer n::int
- = ccmatrix p (1,n);
-cimatrix p::pointer n::int
- = cimatrix p (1,n);
+dmatrix_from_array p::pointer n::int
+ = dmatrix_from_array p (1,n);
+cmatrix_from_array p::pointer n::int
+ = cmatrix_from_array p (1,n);
+imatrix_from_array p::pointer n::int
+ = imatrix_from_array p (1,n);
-/* Matrix comparisons. */
+/* Create a numeric matrix from a pointer, copying the data to fresh memory.
+ The source pointer p may also be NULL, in which case the new matrix is
+ filled with zeros instead. */
-x::matrix == y::matrix = x === y
- if nmatrixp x && matrix_type x == matrix_type y;
-// mixed numeric cases
- = cmatrix x === y if nmatrixp x && cmatrixp y;
- = x === cmatrix y if cmatrixp x && nmatrixp y;
- = dmatrix x === y if imatrixp x && dmatrixp y;
- = x === dmatrix y if dmatrixp x && imatrixp y;
-// comparisons with symbolic matrices
- = 0 if dim x != dim y;
- = compare 0 with
- compare i::int = 1 if i>=n;
- = 0 if x!i != y!i;
- = compare (i+1);
- end when n::int = #x end;
+private matrix_from_double_array_dup;
+private matrix_from_complex_array_dup;
+private matrix_from_int_array_dup;
+extern expr* matrix_from_double_array_dup(void* p, int n, int m);
+extern expr* matrix_from_complex_array_dup(void* p, int n, int m);
+extern expr* matrix_from_int_array_dup(void* p, int n, int m);
-x::matrix != y::matrix = not x == y;
+dmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_double_array_dup p n m;
+cmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_complex_array_dup p n m;
+imatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_int_array_dup p n m;
+
+dmatrix_from_array_dup p::pointer n::int
+ = dmatrix_from_array_dup p (1,n);
+cmatrix_from_array_dup p::pointer n::int
+ = cmatrix_from_array_dup p (1,n);
+imatrix_from_array_dup p::pointer n::int
+ = imatrix_from_array_dup p (1,n);
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 22:03:18 UTC (rev 843)
+++ pure/trunk/runtime.cc 2008-09-23 23:25:21 UTC (rev 844)
@@ -388,19 +388,19 @@
#ifdef HAVE_GSL
case EXPR::DMATRIX: {
gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
- m->owner = owner;
+ m->owner = owner && m->block;
gsl_matrix_free(m);
break;
}
case EXPR::CMATRIX: {
gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
- m->owner = owner;
+ m->owner = owner && m->block;
gsl_matrix_complex_free(m);
break;
}
case EXPR::IMATRIX: {
gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
- m->owner = owner;
+ m->owner = owner && m->block;
gsl_matrix_int_free(m);
break;
}
@@ -5016,14 +5016,94 @@
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_double_matrix(create_double_matrix(n1, n2));
- if (!p) p = calloc(n1*n2, sizeof(double));
if (!p) return 0;
gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
// take a copy of the view matrix
gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
- gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
assert(m && v.matrix.data);
*m = v.matrix;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::DMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ if (!p) return 0;
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ assert(m && v.matrix.data);
+ *m = v.matrix;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_double_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(double));
+ else {
+ void *q = malloc(n1*n2*sizeof(double));
+ memcpy(q, p, n1*n2*sizeof(double));
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
b->size = n1*n2;
b->data = m->data;
m->block = b;
@@ -5040,12 +5120,18 @@
}
extern "C"
-pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n1, uint32_t n2)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_complex_matrix(create_complex_matrix(n1, n2));
- if (!p) p = calloc(2*n1*n2, sizeof(double));
+ if (!p)
+ p = calloc(2*n1*n2, sizeof(double));
+ else {
+ void *q = malloc(2*n1*n2*sizeof(double));
+ memcpy(q, p, 2*n1*n2*sizeof(double));
+ p = q;
+ }
if (!p) return 0;
gsl_matrix_complex_view v =
gsl_matrix_complex_view_array((double*)p, n1, n2);
@@ -5053,7 +5139,7 @@
gsl_matrix_complex *m =
(gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
gsl_block_complex *b = (gsl_block_complex*)malloc(sizeof(gsl_block_complex));
- assert(m && v.matrix.data);
+ assert(m && b && v.matrix.data);
*m = v.matrix;
b->size = n1*n2;
b->data = m->data;
@@ -5071,18 +5157,24 @@
}
extern "C"
-pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_int_array_dup(void *p, uint32_t n1, uint32_t n2)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_int_matrix(create_int_matrix(n1, n2));
- if (!p) p = calloc(n1*n2, sizeof(int));
+ if (!p)
+ p = calloc(n1*n2, sizeof(int));
+ else {
+ void *q = malloc(n1*n2*sizeof(int));
+ memcpy(q, p, n1*n2*sizeof(int));
+ p = q;
+ }
if (!p) return 0;
gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
// take a copy of the view matrix
gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
- assert(m && v.matrix.data);
+ assert(m && b && v.matrix.data);
*m = v.matrix;
b->size = n1*n2;
b->data = m->data;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-23 22:03:18 UTC (rev 843)
+++ pure/trunk/runtime.h 2008-09-23 23:25:21 UTC (rev 844)
@@ -734,17 +734,25 @@
pure_expr *matrix_im(pure_expr *x);
/* Create a matrix object of the given dimensions which uses the given pointer
- p as its underlying C array. There are no checks whatsoever, so the caller
- is responsible for making sure that the memory has the right size and is
- properly initialized. p must point to a malloc'ed memory area, which is
- taken over by Pure and will be freed automatically when the matrix is
- garbage-collected. p may also be NULL in which case a suitable pointer is
- allocated automatically. */
+ p as its underlying storage. There are no checks whatsoever and the data is
+ *not* copied, so the caller is responsible for making sure that the memory
+ has the right size, is properly initialized and is not freed while the
+ matrix is still in use. The memory is *not* freed when the matrix is
+ garbage-collected but remains in the ownership of the caller. */
pure_expr *matrix_from_double_array(void *p, uint32_t n, uint32_t m);
pure_expr *matrix_from_complex_array(void *p, uint32_t n, uint32_t m);
pure_expr *matrix_from_int_array(void *p, uint32_t n, uint32_t m);
+/* The following routines work like the above, but copy the data to newly
+ allocated memory, so the original data can be freed after the call.
+ Moreover, the source pointer p may also be NULL in which case the new
+ matrix is filled with zeros instead. */
+
+pure_expr *matrix_from_double_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_int_array_dup(void *p, uint32_t n, uint32_t m);
+
/* 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-23 23:57:04
|
Revision: 845
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=845&view=rev
Author: agraef
Date: 2008-09-23 23:56:53 +0000 (Tue, 23 Sep 2008)
Log Message:
-----------
Add some more matrix-pointer operations for alternative base types.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-23 23:25:21 UTC (rev 844)
+++ pure/trunk/lib/matrices.pure 2008-09-23 23:56:53 UTC (rev 845)
@@ -421,3 +421,34 @@
= cmatrix_from_array_dup p (1,n);
imatrix_from_array_dup p::pointer n::int
= imatrix_from_array_dup p (1,n);
+
+/* Some helper routines for dealing with alternate base types. These work like
+ the routines above, but initialize the data from float, complex float,
+ short and byte arrays, respectively. */
+
+private matrix_from_float_array_dup;
+private matrix_from_complex_float_array_dup;
+private matrix_from_short_array_dup;
+private matrix_from_byte_array_dup;
+extern expr* matrix_from_float_array_dup(void* p, int n, int m);
+extern expr* matrix_from_complex_float_array_dup(void* p, int n, int m);
+extern expr* matrix_from_short_array_dup(void* p, int n, int m);
+extern expr* matrix_from_byte_array_dup(void* p, int n, int m);
+
+fmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_float_array_dup p n m;
+cfmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_complex_float_array_dup p n m;
+smatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_short_array_dup p n m;
+bmatrix_from_array_dup p::pointer (n::int,m::int)
+ = matrix_from_byte_array_dup p n m;
+
+fmatrix_from_array_dup p::pointer n::int
+ = fmatrix_from_array_dup p (1,n);
+cfmatrix_from_array_dup p::pointer n::int
+ = cfmatrix_from_array_dup p (1,n);
+smatrix_from_array_dup p::pointer n::int
+ = smatrix_from_array_dup p (1,n);
+bmatrix_from_array_dup p::pointer n::int
+ = bmatrix_from_array_dup p (1,n);
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-23 23:25:21 UTC (rev 844)
+++ pure/trunk/runtime.cc 2008-09-23 23:56:53 UTC (rev 845)
@@ -5191,6 +5191,152 @@
#endif
}
+extern "C"
+pure_expr *matrix_from_float_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(double));
+ else {
+ void *q = malloc(n1*n2*sizeof(double));
+ float *p1 = (float*)p; double *q1 = (double*)q;
+ for (size_t i = 0; i < n1*n2; i++) q1[i] = (double)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_view v = gsl_matrix_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ gsl_block *b = (gsl_block*)malloc(sizeof(gsl_block));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::DMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_complex_float_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ if (!p)
+ p = calloc(2*n1*n2, sizeof(double));
+ else {
+ void *q = malloc(2*n1*n2*sizeof(double));
+ float *p1 = (float*)p; double *q1 = (double*)q;
+ for (size_t i = 0; i < 2*n1*n2; i++) q1[i] = (double)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_view_array((double*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ gsl_block_complex *b = (gsl_block_complex*)malloc(sizeof(gsl_block_complex));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::CMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_short_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(int));
+ else {
+ void *q = malloc(n1*n2*sizeof(int));
+ short *p1 = (short*)p; int *q1 = (int*)q;
+ for (size_t i = 0; i < n1*n2; i++) q1[i] = (int)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_byte_array_dup(void *p, uint32_t n1, uint32_t n2)
+{
+#ifdef HAVE_GSL
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ if (!p)
+ p = calloc(n1*n2, sizeof(int));
+ else {
+ void *q = malloc(n1*n2*sizeof(int));
+ int8_t *p1 = (int8_t*)p; int *q1 = (int*)q;
+ for (size_t i = 0; i < n1*n2; i++) q1[i] = (int)p1[i];
+ p = q;
+ }
+ if (!p) return 0;
+ gsl_matrix_int_view v = gsl_matrix_int_view_array((int*)p, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ gsl_block_int *b = (gsl_block_int*)malloc(sizeof(gsl_block_int));
+ assert(m && b && v.matrix.data);
+ *m = v.matrix;
+ b->size = n1*n2;
+ b->data = m->data;
+ m->block = b;
+ pure_expr *x = new_expr();
+ x->tag = EXPR::IMATRIX;
+ x->data.mat.p = m;
+ x->data.mat.refc = new uint32_t;
+ *x->data.mat.refc = 1;
+ MEMDEBUG_NEW(x)
+ return x;
+#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-23 23:25:21 UTC (rev 844)
+++ pure/trunk/runtime.h 2008-09-23 23:56:53 UTC (rev 845)
@@ -746,6 +746,7 @@
/* The following routines work like the above, but copy the data to newly
allocated memory, so the original data can be freed after the call.
+ Additional routines are provided for various alternate data sizes.
Moreover, the source pointer p may also be NULL in which case the new
matrix is filled with zeros instead. */
@@ -753,6 +754,11 @@
pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n, uint32_t m);
pure_expr *matrix_from_int_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_float_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_complex_float_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_short_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_byte_array_dup(void *p, uint32_t n, uint32_t m);
+
/* 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-25 00:55:21
|
Revision: 851
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=851&view=rev
Author: agraef
Date: 2008-09-25 00:55:11 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Add Gaussian elimination example.
Modified Paths:
--------------
pure/trunk/NEWS
Added Paths:
-----------
pure/trunk/examples/gauss.pure
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-25 00:54:30 UTC (rev 850)
+++ pure/trunk/NEWS 2008-09-25 00:55:11 UTC (rev 851)
@@ -78,7 +78,7 @@
- #x total number of elements
- dim x number of rows and columns (as a pair)
-- x' transposed matrix
+- x' transpose of a matrix
- x!i ith matrix element in row-major order (zero-based)
- x!(i,j) jth element in ith row of the matrix (zero-based)
- x!!is, x!!(is,js) slicing (is and js are lists of machine ints)
Added: pure/trunk/examples/gauss.pure
===================================================================
--- pure/trunk/examples/gauss.pure (rev 0)
+++ pure/trunk/examples/gauss.pure 2008-09-25 00:55:11 UTC (rev 851)
@@ -0,0 +1,69 @@
+
+/* Pure matrix example: Gaussian elimination algorithm. 2008-09-25 AG */
+
+/* The Gaussian elimination algorithm brings a matrix into row echelon form,
+ which makes it possible to solve the original system of linear equations
+ using back substitution. Our version of the algorithm just returns the row
+ echelon form of the matrix, along with the corresponding permutation of the
+ row indices performed during pivoting. */
+
+gauss_elimination x::matrix = p,x
+when n,m = dim x; p,_,x = foldl step (0..n-1,0,x) (0..m-1) end;
+
+/* The elimination step. x is our matrix, i the current row index, j the
+ current column index, and p keeps track of the current permutation of the
+ row indices performed during pivoting. The algorithm returns the updated
+ matrix x, row index i and row permutation p. */
+
+/* Here is a brief rundown of the algorithm: First we find the pivot element
+ in column j of the matrix. (We're doing partial pivoting here, i.e., we
+ only look for the pivot in column j, starting at row i.) If the pivot is
+ zero then we're done (the pivot column is already zeroed out). Otherwise,
+ we bring it into the pivot position (swapping row i and the pivot row),
+ divide the picot row by the pivot, and subtract suitable multiples of the
+ pivot row to eliminate the pivot column in all subsequent rows. Finally we
+ update i and p accordingly and return the result. */
+
+step (p,i,x) j
+= if max_x>0 then
+ // updated row permutation and index:
+ transp i max_i p, i+1,
+ {// the top rows of the matrix remain unchanged:
+ x!!(0..i-1,0..m-1);
+ // the pivot row, divided by the pivot:
+ {x!(i,l)/x!(i,j) | l=0..m-1};
+ // subtract suitable multiples of the pivot row:
+ {x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
+ else p,i,x
+when
+ n,m = dim x; max_i, max_x = pivot i (col x j);
+ x = if max_x>0 then swap x i max_i else x;
+end with
+ pivot i x = foldl max (0,0) [j,abs (x!j)|j=i..#x-1];
+ max (i,x) (j,y) = if x<y then j,y else i,x;
+end;
+
+/* Swap rows i and j of the matrix x. This operation is used in the partial
+ pivoting step. (This is not really needed because we could just index the
+ rows indirectly through the row permutation readily available in our
+ implementation of the algorithm, but we omitted this here for clarity.) */
+
+swap x i j = x!!(transp i j (0..n-1),0..m-1) when n,m = dim x end;
+
+/* A little helper function to apply a transposition to a permutation. */
+
+transp i j p = [p!tr k | k=0..#p-1]
+with tr k = if k==i then j else if k==j then i else k end;
+
+/* For convenience, print a double matrix in "short" format a la Octave. */
+
+using system;
+__show__ x::matrix
+= strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\n"
+with printd 0 = sprintf "\n%10.5f"; printd _ = sprintf "%10.5f" end
+when n,m = dim x end if dmatrixp x;
+
+/* Example: */
+
+let x = dmatrix {2,1,-1,8; -3,-1,2,-11; -2,1,2,-3};
+x; gauss_elimination x;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-25 10:53:00
|
Revision: 863
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=863&view=rev
Author: agraef
Date: 2008-09-25 10:52:57 +0000 (Thu, 25 Sep 2008)
Log Message:
-----------
Cosmetic changes.
Modified Paths:
--------------
pure/trunk/examples/gauss.pure
pure/trunk/pure.1.in
Modified: pure/trunk/examples/gauss.pure
===================================================================
--- pure/trunk/examples/gauss.pure 2008-09-25 10:48:21 UTC (rev 862)
+++ pure/trunk/examples/gauss.pure 2008-09-25 10:52:57 UTC (rev 863)
@@ -26,7 +26,8 @@
i and p accordingly and return the result. */
step (p,i,x) j
-= if max_x>0 then
+= if max_x==0 then p,i,x
+ else
// updated row permutation and index:
transp i max_i p, i+1,
{// the top rows of the matrix remain unchanged:
@@ -35,7 +36,6 @@
{x!(i,l)/x!(i,j) | l=0..m-1};
// subtract suitable multiples of the pivot row:
{x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
- else p,i,x
when
n,m = dim x; max_i, max_x = pivot i (col x j);
x = if max_x>0 then swap x i max_i else x;
Modified: pure/trunk/pure.1.in
===================================================================
--- pure/trunk/pure.1.in 2008-09-25 10:48:21 UTC (rev 862)
+++ pure/trunk/pure.1.in 2008-09-25 10:52:57 UTC (rev 863)
@@ -1261,7 +1261,8 @@
// One pivoting and elimination step in column j of the matrix:
step (p,i,x) j
-= \fBif\fP max_x>0 \fBthen\fP
+= \fBif\fP max_x==0 \fBthen\fP p,i,x
+ \fBelse\fP
// updated row permutation and index:
transp i max_i p, i+1,
{// the top rows of the matrix remain unchanged:
@@ -1270,7 +1271,6 @@
{x!(i,l)/x!(i,j) | l=0..m-1};
// subtract suitable multiples of the pivot row:
{x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}}
- \fBelse\fP p,i,x
\fBwhen\fP
n,m = dim x; max_i, max_x = pivot i (col x j);
x = \fBif\fP max_x>0 \fBthen\fP swap x i max_i \fBelse\fP x;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-27 16:03:53
|
Revision: 884
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=884&view=rev
Author: agraef
Date: 2008-09-27 16:03:46 +0000 (Sat, 27 Sep 2008)
Log Message:
-----------
Added missing complex->double/int matrix conversions, overhaul of matrix<->pointer conversions.
Modified Paths:
--------------
pure/trunk/lib/matrices.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/matrices.pure
===================================================================
--- pure/trunk/lib/matrices.pure 2008-09-27 13:11:09 UTC (rev 883)
+++ pure/trunk/lib/matrices.pure 2008-09-27 16:03:46 UTC (rev 884)
@@ -217,9 +217,9 @@
(either a pair denoting the number of rows and columns, or just the row
size in order to create a row vector). */
-dmatrix (n::int,m::int) = dmatrix_dup (pointer 0,n,m);
-cmatrix (n::int,m::int) = cmatrix_dup (pointer 0,n,m);
-imatrix (n::int,m::int) = imatrix_dup (pointer 0,n,m);
+dmatrix (n::int,m::int) = double_matrix (n,m) (pointer 0);
+cmatrix (n::int,m::int) = complex_matrix (n,m) (pointer 0);
+imatrix (n::int,m::int) = int_matrix (n,m) (pointer 0);
dmatrix n::int = dmatrix (1,n);
cmatrix n::int = cmatrix (1,n);
@@ -233,8 +233,8 @@
extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
expr* matrix_int(expr *x);
-dmatrix x::matrix = matrix_double x if imatrixp x || dmatrixp x;
-imatrix x::matrix = matrix_int x if imatrixp x || dmatrixp x;
+dmatrix x::matrix = matrix_double x if nmatrixp x;
+imatrix x::matrix = matrix_int x if nmatrixp x;
cmatrix x::matrix = matrix_complex x if nmatrixp x;
private matrix_re matrix_im;
@@ -374,94 +374,126 @@
important consideration (e.g., graphics, video and audio applications). The
usual caveats apply. */
-/* Convert a matrix to a pointer, which points directly to the underlying C
- array. You have to be careful when passing such a pointer to C functions if
- the underlying data is non-contiguous; when in doubt, first use the 'pack'
- function from above to place the data in contiguous storage. */
+/* Get a pointer to the underlying C array of a matrix. The data is *not*
+ copied. Hence you have to be careful when passing such a pointer to C
+ functions if the underlying data is non-contiguous; when in doubt, first
+ use the 'pack' function from above to place the data in contiguous storage,
+ or use one of the matrix-pointer conversion routines below. */
private pure_pointerval;
extern expr* pure_pointerval(expr*);
pointer x::matrix = pure_pointerval x;
-/* Create a numeric matrix from a pointer, without copying the data. The
- caller must ensure that the pointer points to properly initialized memory
- big enough to accommodate the requested dimensions, which persists for the
- entire lifetime of the matrix object. */
+/* Copy the contents of a matrix to a given pointer and return that pointer,
+ converting to the target data type on the fly if necessary. The given
+ pointer may also be NULL, in which case suitable memory is malloc'ed and
+ returned; otherwise the caller must ensure that the memory pointed to by p
+ is big enough for the contents of the given matrix. */
+private matrix_to_double_array;
+private matrix_to_float_array;
+private matrix_to_complex_array;
+private matrix_to_complex_float_array;
+private matrix_to_int_array;
+private matrix_to_short_array;
+private matrix_to_byte_array;
+extern void* matrix_to_double_array(void* p, expr* x);
+extern void* matrix_to_float_array(void* p, expr* x);
+extern void* matrix_to_complex_array(void* p, expr* x);
+extern void* matrix_to_complex_float_array(void* p, expr* x);
+extern void* matrix_to_int_array(void* p, expr* x);
+extern void* matrix_to_short_array(void* p, expr* x);
+extern void* matrix_to_byte_array(void* p, expr* x);
+
+double_pointer p::pointer x::matrix
+ = matrix_to_double_array p x if nmatrixp x;
+float_pointer p::pointer x::matrix
+ = matrix_to_float_array p x if nmatrixp x;
+complex_pointer p::pointer x::matrix
+ = matrix_to_complex_array p x if nmatrixp x;
+complex_float_pointer p::pointer x::matrix
+ = matrix_to_complex_float_array p x if nmatrixp x;
+int_pointer p::pointer x::matrix
+ = matrix_to_int_array p x if nmatrixp x;
+short_pointer p::pointer x::matrix
+ = matrix_to_short_array p x if nmatrixp x;
+byte_pointer p::pointer x::matrix
+ = matrix_to_byte_array p x if nmatrixp x;
+
+/* Create a numeric matrix from a pointer, copying the data and converting it
+ from the source type on the fly if necessary. The source pointer p may also
+ be NULL, in which case the new matrix is filled with zeros instead.
+ Otherwise the caller must ensure that the pointer points to properly
+ initialized memory big enough for the requested dimensions. */
+
private matrix_from_double_array;
+private matrix_from_float_array;
private matrix_from_complex_array;
+private matrix_from_complex_float_array;
private matrix_from_int_array;
-extern expr* matrix_from_double_array(void* p, int n, int m);
-extern expr* matrix_from_complex_array(void* p, int n, int m);
-extern expr* matrix_from_int_array(void* p, int n, int m);
+private matrix_from_short_array;
+private matrix_from_byte_array;
+extern expr* matrix_from_double_array(int n, int m, void* p);
+extern expr* matrix_from_float_array(int n, int m, void* p);
+extern expr* matrix_from_complex_array(int n, int m, void* p);
+extern expr* matrix_from_complex_float_array(int n, int m, void* p);
+extern expr* matrix_from_int_array(int n, int m, void* p);
+extern expr* matrix_from_short_array(int n, int m, void* p);
+extern expr* matrix_from_byte_array(int n, int m, void* p);
-dmatrix (p::pointer,n::int,m::int)
- = matrix_from_double_array p n m;
-cmatrix (p::pointer,n::int,m::int)
- = matrix_from_complex_array p n m;
-imatrix (p::pointer,n::int,m::int)
- = matrix_from_int_array p n m;
+double_matrix (n::int,m::int) p::pointer
+ = matrix_from_double_array n m p;
+float_matrix (n::int,m::int) p::pointer
+ = matrix_from_float_array n m p;
+complex_matrix (n::int,m::int) p::pointer
+ = matrix_from_complex_array n m p;
+complex_float_matrix (n::int,m::int) p::pointer
+ = matrix_from_complex_float_array n m p;
+int_matrix (n::int,m::int) p::pointer
+ = matrix_from_int_array n m p;
+short_matrix (n::int,m::int) p::pointer
+ = matrix_from_short_array n m p;
+byte_matrix (n::int,m::int) p::pointer
+ = matrix_from_byte_array n m p;
-dmatrix (p::pointer,n::int)
- = dmatrix (p,1,n);
-cmatrix (p::pointer,n::int)
- = cmatrix (p,1,n);
-imatrix (p::pointer,n::int)
- = imatrix (p,1,n);
+double_matrix n::int p::pointer
+ = double_matrix (1,n) p;
+float_matrix n::int p::pointer
+ = float_matrix (1,n) p;
+complex_matrix n::int p::pointer
+ = complex_matrix (1,n) p;
+complex_float_matrix n::int p::pointer
+ = complex_float_matrix (1,n) p;
+int_matrix n::int p::pointer
+ = int_matrix (1,n) p;
+short_matrix n::int p::pointer
+ = short_matrix (1,n) p;
+byte_matrix n::int p::pointer
+ = byte_matrix (1,n) p;
-/* Create a numeric matrix from a pointer, copying the data to fresh memory.
- The source pointer p may also be NULL, in which case the new matrix is
- filled with zeros instead. */
+/* Create a numeric matrix view of existing data, without copying the data.
+ The data must be double, complex or int, the pointer must not be NULL and
+ the caller must also ensure that the memory persists for the entire
+ lifetime of the matrix object. */
-private matrix_from_double_array_dup;
-private matrix_from_complex_array_dup;
-private matrix_from_int_array_dup;
-extern expr* matrix_from_double_array_dup(void* p, int n, int m);
-extern expr* matrix_from_complex_array_dup(void* p, int n, int m);
-extern expr* matrix_from_int_array_dup(void* p, int n, int m);
+private matrix_from_double_array_nodup;
+private matrix_from_complex_array_nodup;
+private matrix_from_int_array_nodup;
+extern expr* matrix_from_double_array_nodup(int n, int m, void* p);
+extern expr* matrix_from_complex_array_nodup(int n, int m, void* p);
+extern expr* matrix_from_int_array_nodup(int n, int m, void* p);
-dmatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_double_array_dup p n m;
-cmatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_complex_array_dup p n m;
-imatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_int_array_dup p n m;
+double_matrix_view (n::int,m::int) p::pointer
+ = matrix_from_double_array_nodup n m p;
+complex_matrix_view (n::int,m::int) p::pointer
+ = matrix_from_complex_array_nodup n m p;
+int_matrix_view (n::int,m::int) p::pointer
+ = matrix_from_int_array_nodup n m p;
-dmatrix_dup (p::pointer,n::int)
- = dmatrix_dup (p,1,n);
-cmatrix_dup (p::pointer,n::int)
- = cmatrix_dup (p,1,n);
-imatrix_dup (p::pointer,n::int)
- = imatrix_dup (p,1,n);
-
-/* Some additional functions for alternate base types. These work like the
- routines above, but initialize the data from float, complex float, short
- and byte arrays, respectively. */
-
-private matrix_from_float_array_dup;
-private matrix_from_complex_float_array_dup;
-private matrix_from_short_array_dup;
-private matrix_from_byte_array_dup;
-extern expr* matrix_from_float_array_dup(void* p, int n, int m);
-extern expr* matrix_from_complex_float_array_dup(void* p, int n, int m);
-extern expr* matrix_from_short_array_dup(void* p, int n, int m);
-extern expr* matrix_from_byte_array_dup(void* p, int n, int m);
-
-float_dmatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_float_array_dup p n m;
-float_cmatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_complex_float_array_dup p n m;
-short_imatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_short_array_dup p n m;
-byte_imatrix_dup (p::pointer,n::int,m::int)
- = matrix_from_byte_array_dup p n m;
-
-float_dmatrix_dup (p::pointer,n::int)
- = float_dmatrix_dup (p,1,n);
-float_cmatrix_dup (p::pointer,n::int)
- = float_cmatrix_dup (p,1,n);
-short_imatrix_dup (p::pointer,n::int)
- = short_imatrix_dup (p,1,n);
-byte_imatrix_dup (p::pointer,n::int)
- = byte_imatrix_dup (p,1,n);
+double_matrix_view n::int p::pointer
+ = double_matrix_view (1,n) p;
+complex_matrix_view n::int p::pointer
+ = complex_matrix_view (1,n) p;
+int_matrix_view n::int p::pointer
+ = int_matrix_view (1,n) p;
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-27 13:11:09 UTC (rev 883)
+++ pure/trunk/runtime.cc 2008-09-27 16:03:46 UTC (rev 884)
@@ -4882,6 +4882,19 @@
m2->data[i*m2->tda+j] = (double)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 *m2 = create_double_matrix(n, 2*m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = i*m2->tda+2*j;
+ size_t l = 2*(i*m1->tda+j);
+ m2->data[k] = m1->data[l];
+ m2->data[k+1] = m1->data[l+1];
+ }
+ return pure_double_matrix(m2);
+ }
default:
return 0;
}
@@ -4945,6 +4958,19 @@
}
case EXPR::IMATRIX:
return x;
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(n, 2*m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = i*m2->tda+2*j;
+ size_t l = 2*(i*m1->tda+j);
+ m2->data[k] = (int)m1->data[l];
+ m2->data[k+1] = (int)m1->data[l+1];
+ }
+ return pure_int_matrix(m2);
+ }
default:
return 0;
}
@@ -5019,7 +5045,7 @@
}
extern "C"
-pure_expr *matrix_from_double_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_double_array_nodup(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5043,7 +5069,7 @@
}
extern "C"
-pure_expr *matrix_from_complex_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_complex_array_nodup(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5069,7 +5095,7 @@
}
extern "C"
-pure_expr *matrix_from_int_array(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_int_array_nodup(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5093,7 +5119,7 @@
}
extern "C"
-pure_expr *matrix_from_double_array_dup(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_double_array(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5128,7 +5154,7 @@
}
extern "C"
-pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_complex_array(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5165,7 +5191,7 @@
}
extern "C"
-pure_expr *matrix_from_int_array_dup(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_int_array(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5200,9 +5226,178 @@
}
extern "C"
-pure_expr *matrix_from_float_array_dup(void *p, uint32_t n1, uint32_t n2)
+void *matrix_to_double_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(double));
+ if (!p) return 0;
+ double *q = (double*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = m1->data[l];
+ q[k++] = m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(double));
+ if (!p) return 0;
+ double *q = (double*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = m1->data[i*m1->tda+j];
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(double));
+ if (!p) return 0;
+ double *q = (double*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (double)m1->data[i*m1->tda+j];
+ return p;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+void *matrix_to_complex_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(double));
+ if (!p) return 0;
+ double *q = (double*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = m1->data[l];
+ q[k++] = m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(double));
+ if (!p) return 0;
+ double *q = (double*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ q[k++] = m1->data[i*m1->tda+j];
+ q[k++] = 0.0;
+ }
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(double));
+ if (!p) return 0;
+ double *q = (double*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ q[k++] = (double)m1->data[i*m1->tda+j];
+ q[k++] = 0.0;
+ }
+ return p;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+void *matrix_to_int_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(int));
+ if (!p) return 0;
+ int *q = (int*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = (int)m1->data[l];
+ q[k++] = (int)m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(int));
+ if (!p) return 0;
+ int *q = (int*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (int)m1->data[i*m1->tda+j];
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(int));
+ if (!p) return 0;
+ int *q = (int*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = m1->data[i*m1->tda+j];
+ return p;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_from_float_array(uint32_t n1, uint32_t n2, void *p)
+{
+#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
return pure_double_matrix(create_double_matrix(n1, n2));
if (!p)
@@ -5236,7 +5431,7 @@
}
extern "C"
-pure_expr *matrix_from_complex_float_array_dup(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_complex_float_array(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5274,7 +5469,7 @@
}
extern "C"
-pure_expr *matrix_from_short_array_dup(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_short_array(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5310,7 +5505,7 @@
}
extern "C"
-pure_expr *matrix_from_byte_array_dup(void *p, uint32_t n1, uint32_t n2)
+pure_expr *matrix_from_byte_array(uint32_t n1, uint32_t n2, void *p)
{
#ifdef HAVE_GSL
if (n1 == 0 || n2 == 0) // empty matrix
@@ -5345,6 +5540,230 @@
#endif
}
+extern "C"
+void *matrix_to_float_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(float));
+ if (!p) return 0;
+ float *q = (float*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = (float)m1->data[l];
+ q[k++] = (float)m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(float));
+ if (!p) return 0;
+ float *q = (float*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (float)m1->data[i*m1->tda+j];
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(float));
+ if (!p) return 0;
+ float *q = (float*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (float)m1->data[i*m1->tda+j];
+ return p;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+void *matrix_to_complex_float_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(float));
+ if (!p) return 0;
+ float *q = (float*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = (float)m1->data[l];
+ q[k++] = (float)m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(float));
+ if (!p) return 0;
+ float *q = (float*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ q[k++] = (float)m1->data[i*m1->tda+j];
+ q[k++] = 0.0;
+ }
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(float));
+ if (!p) return 0;
+ float *q = (float*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ q[k++] = (float)m1->data[i*m1->tda+j];
+ q[k++] = 0.0;
+ }
+ return p;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+void *matrix_to_short_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(short));
+ if (!p) return 0;
+ short *q = (short*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = (short)m1->data[l];
+ q[k++] = (short)m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(short));
+ if (!p) return 0;
+ short *q = (short*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (short)m1->data[i*m1->tda+j];
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(short));
+ if (!p) return 0;
+ short *q = (short*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (short)m1->data[i*m1->tda+j];
+ return p;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+void *matrix_to_byte_array(void *p, 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;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(2*n*m*sizeof(int8_t));
+ if (!p) return 0;
+ int8_t *q = (int8_t*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t l = 2*(i*m1->tda+j);
+ q[k++] = (int8_t)m1->data[l];
+ q[k++] = (int8_t)m1->data[l+1];
+ }
+ return p;
+ }
+ case EXPR::DMATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(int8_t));
+ if (!p) return 0;
+ int8_t *q = (int8_t*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (int8_t)m1->data[i*m1->tda+j];
+ return p;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ if (n==0 || m==0) return p;
+ if (!p) p = malloc(n*m*sizeof(int8_t));
+ if (!p) return 0;
+ int8_t *q = (int8_t*)p;
+ size_t k = 0;
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ q[k++] = (int8_t)m1->data[i*m1->tda+j];
+ return p;
+ }
+ 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-27 13:11:09 UTC (rev 883)
+++ pure/trunk/runtime.h 2008-09-27 16:03:46 UTC (rev 884)
@@ -715,10 +715,7 @@
pure_expr *matrix_transpose(pure_expr *x);
-/* Convert between different types of numeric matrices. Note that any numeric
- matrix can be converted to a complex matrix, but for the other conversions
- the input must be a a double or integer matrix (see matrix_re and matrix_im
- below to handle the complex->double case). */
+/* Convert between different types of numeric matrices. */
pure_expr *matrix_double(pure_expr *x);
pure_expr *matrix_complex(pure_expr *x);
@@ -740,25 +737,42 @@
matrix is still in use. The memory is *not* freed when the matrix is
garbage-collected but remains in the ownership of the caller. */
-pure_expr *matrix_from_double_array(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_complex_array(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_int_array(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_double_array_nodup(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_complex_array_nodup(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_int_array_nodup(uint32_t n, uint32_t m, void *p);
/* The following routines work like the above, but copy the data to newly
allocated memory, so the original data can be freed after the call.
- Additional routines are provided for various alternate data sizes.
Moreover, the source pointer p may also be NULL in which case the new
matrix is filled with zeros instead. */
-pure_expr *matrix_from_double_array_dup(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_complex_array_dup(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_int_array_dup(void *p, uint32_t n, uint32_t m);
+pure_expr *matrix_from_double_array(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_complex_array(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_int_array(uint32_t n, uint32_t m, void *p);
-pure_expr *matrix_from_float_array_dup(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_complex_float_array_dup(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_short_array_dup(void *p, uint32_t n, uint32_t m);
-pure_expr *matrix_from_byte_array_dup(void *p, uint32_t n, uint32_t m);
+/* Copy data from the given matrix to the given data pointer, which is then
+ returned. If p is NULL then memory of sufficient size is malloc'ed;
+ otherwise p must point to a memory area of sufficient size. */
+void *matrix_to_double_array(void *p, pure_expr *x);
+void *matrix_to_complex_array(void *p, pure_expr *x);
+void *matrix_to_int_array(void *p, pure_expr *x);
+
+/* Additional routines for alternative base types. These work like the
+ routines above but take data consisting of base types which are not
+ directly supported by Pure GSL matrices: float, complex float, short,
+ byte. */
+
+pure_expr *matrix_from_float_array(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_complex_float_array(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_short_array(uint32_t n, uint32_t m, void *p);
+pure_expr *matrix_from_byte_array(uint32_t n, uint32_t m, void *p);
+
+void *matrix_to_float_array(void *p, pure_expr *x);
+void *matrix_to_complex_float_array(void *p, pure_expr *x);
+void *matrix_to_short_array(void *p, pure_expr *x);
+void *matrix_to_byte_array(void *p, 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-27 16:18:49
|
Revision: 885
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=885&view=rev
Author: agraef
Date: 2008-09-27 16:18:43 +0000 (Sat, 27 Sep 2008)
Log Message:
-----------
Bump version number, update ChangeLog and NEWS.
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/NEWS
pure/trunk/configure
pure/trunk/configure.ac
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-27 16:03:46 UTC (rev 884)
+++ pure/trunk/ChangeLog 2008-09-27 16:18:43 UTC (rev 885)
@@ -1,3 +1,16 @@
+2008-09-27 Albert Graef <Dr....@t-...>
+
+ * configure.ac: Bump version number.
+
+ * lib/matrices.pure, interpreter.cc/h: Added missing
+ complex->double/int matrix conversions. Thorough overhaul of
+ matrix<->pointer conversions, which now provide a complete set of
+ conversions from/to any reasonable C array representation,
+ including float, complex float, short and byte arrays.
+
+ * lib/prelude.pure: Fix a bug in cat which broke catmap on
+ strings. Reported by Eddie Rucker.
+
2008-09-25 Albert Graef <Dr....@t-...>
* 0.7 release.
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-27 16:03:46 UTC (rev 884)
+++ pure/trunk/NEWS 2008-09-27 16:18:43 UTC (rev 885)
@@ -1,6 +1,14 @@
-** Pure 0.7 (2008-09-26)
+** Pure 0.8 (in progress)
+This is a maintenance release. It fixes a rather annoying bug in catmap
+(reported by Eddie Rucker) and some minor glitches in the Makefile and the
+documentation. Also, the matrix conversion operations in matrices.pure have
+been overhauled (in particular, the matrix pointer conversions are much more
+complete now). Details can be found in the ChangeLog.
+
+** Pure 0.7 2008-09-26
+
This release brings an important new feature: GSL (GNU Scientific Library)
matrix support. Here's a brief overview of the new stuff. For more
information on GSL please refer to http://www.gnu.org/software/gsl.
Modified: pure/trunk/configure
===================================================================
--- pure/trunk/configure 2008-09-27 16:03:46 UTC (rev 884)
+++ pure/trunk/configure 2008-09-27 16:18:43 UTC (rev 885)
@@ -1,6 +1,6 @@
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.61 for pure 0.7.
+# Generated by GNU Autoconf 2.61 for pure 0.8.
#
# Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
# 2002, 2003, 2004, 2005, 2006 Free Software Foundation, Inc.
@@ -572,8 +572,8 @@
# Identity of this package.
PACKAGE_NAME='pure'
PACKAGE_TARNAME='pure'
-PACKAGE_VERSION='0.7'
-PACKAGE_STRING='pure 0.7'
+PACKAGE_VERSION='0.8'
+PACKAGE_STRING='pure 0.8'
PACKAGE_BUGREPORT=''
# Factoring default headers for most tests.
@@ -1199,7 +1199,7 @@
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
-\`configure' configures pure 0.7 to adapt to many kinds of systems.
+\`configure' configures pure 0.8 to adapt to many kinds of systems.
Usage: $0 [OPTION]... [VAR=VALUE]...
@@ -1264,7 +1264,7 @@
if test -n "$ac_init_help"; then
case $ac_init_help in
- short | recursive ) echo "Configuration of pure 0.7:";;
+ short | recursive ) echo "Configuration of pure 0.8:";;
esac
cat <<\_ACEOF
@@ -1358,7 +1358,7 @@
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
-pure configure 0.7
+pure configure 0.8
generated by GNU Autoconf 2.61
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1372,7 +1372,7 @@
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.
-It was created by pure $as_me 0.7, which was
+It was created by pure $as_me 0.8, which was
generated by GNU Autoconf 2.61. Invocation command line was
$ $0 $@
@@ -6716,7 +6716,7 @@
# report actual input values of CONFIG_FILES etc. instead of their
# values after options handling.
ac_log="
-This file was extended by pure $as_me 0.7, which was
+This file was extended by pure $as_me 0.8, which was
generated by GNU Autoconf 2.61. Invocation command line was
CONFIG_FILES = $CONFIG_FILES
@@ -6765,7 +6765,7 @@
_ACEOF
cat >>$CONFIG_STATUS <<_ACEOF
ac_cs_version="\\
-pure config.status 0.7
+pure config.status 0.8
configured by $0, generated by GNU Autoconf 2.61,
with options \\"`echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"
Modified: pure/trunk/configure.ac
===================================================================
--- pure/trunk/configure.ac 2008-09-27 16:03:46 UTC (rev 884)
+++ pure/trunk/configure.ac 2008-09-27 16:18:43 UTC (rev 885)
@@ -2,7 +2,7 @@
dnl To regenerate the configury after changes:
dnl autoconf -I config && autoheader -I config
-AC_INIT(pure, 0.7)
+AC_INIT(pure, 0.8)
AC_CONFIG_AUX_DIR(config)
dnl AC_CONFIG_MACRO_DIR(config)
AC_CONFIG_HEADERS(config.h)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ag...@us...> - 2008-09-28 07:28:06
|
Revision: 890
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=890&view=rev
Author: agraef
Date: 2008-09-28 07:28:01 +0000 (Sun, 28 Sep 2008)
Log Message:
-----------
Final polish (0.8 release).
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/NEWS
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-28 06:37:13 UTC (rev 889)
+++ pure/trunk/ChangeLog 2008-09-28 07:28:01 UTC (rev 890)
@@ -1,5 +1,7 @@
2008-09-28 Albert Graef <Dr....@t-...>
+ * 0.8 release.
+
* test/test025.pure: matrix tests. NOTE: This test is expected to
fail if Pure was built without GSL support.
Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS 2008-09-28 06:37:13 UTC (rev 889)
+++ pure/trunk/NEWS 2008-09-28 07:28:01 UTC (rev 890)
@@ -1,11 +1,11 @@
-** Pure 0.8 (in progress)
+** Pure 0.8 2008-09-28
-This is a maintenance release. It fixes a rather annoying bug in catmap
-(reported by Eddie Rucker) and some minor glitches in the Makefile and the
-documentation. Also, the matrix conversion operations in matrices.pure have
-been overhauled (in particular, the matrix pointer conversions are much more
-complete now). Details can be found in the ChangeLog.
+This is a maintenance release. It fixes a bug in catmap which slipped into the
+0.6 release (reported by Eddie Rucker) and some minor glitches in the Makefile
+and the documentation. Also, the matrix conversion operations in matrices.pure
+have been overhauled (in particular, the matrix-pointer conversions are much
+more complete now). Details can be found in the ChangeLog.
** Pure 0.7 2008-09-26
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X