[pure-lang-svn] SF.net SVN: pure-lang:[844] pure/trunk
Status: Beta
Brought to you by:
agraef
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. |