diff --git a/setup.nqp b/setup.nqp index 1f836b0..d4394b5 100755 --- a/setup.nqp +++ b/setup.nqp @@ -32,8 +32,8 @@ sub MAIN(@argv) { $mode := @argv[0]; } if $mode eq "build" { - probe_for_cblas(%PLA); find_blas(%PLA); + probe_for_cblas_h(%PLA); find_lapack(%PLA); } if $mode eq "test" { @@ -45,6 +45,7 @@ sub MAIN(@argv) { setup_testlib(%PLA); setup_nqp_bootstrapper(%PLA); setup_dynpmc_flags(%PLA); + setup_docs(%PLA); setup(@argv, %PLA); } @@ -70,15 +71,18 @@ sub setup_PLA_keys(%PLA) { %PLA{'inst_lib'} := new_array(); %PLA{'pir_pir'} := new_hash(); %PLA{'pir_pir'}{'t/testlib/pla_test.pir'} := new_array(); + %PLA{'need_cblas_h'} := 0; } # -sub probe_for_cblas(%PLA) { - if probe_include("cblas.h", :verbose(1)) { - pir::say("Cannot find cblas.h\nPlease install libatlas-base-dev"); - pir::exit__vI(1); - } else { - %PLA{'dynpmc_cflags_list'}.push("-D_PLA_HAVE_CBLAS_H"); +sub probe_for_cblas_h(%PLA) { + if %PLA{'need_cblas_h'} { + if probe_include("cblas.h", :verbose(1)) { + pir::say("Cannot find cblas.h. This is required if you are using CBLAS or ATLAS"); + pir::exit(1); + } else { + %PLA{'dynpmc_cflags_list'}.push("-D_PLA_HAVE_CBLAS_H"); + } } } @@ -87,19 +91,7 @@ sub find_blas(%PLA) { my $osname := %config{'osname'}; my $found_blas := 0; if $osname eq 'linux' { - my %searches; - %searches{'/usr/lib/libblas.so'} := ['-lblas', '-D_PLA_HAVE_ATLAS']; - %searches{'/usr/lib/atlas/libcblas.so'} := ['-L/usr/lib/atlas -lcblas', '-D_PLA_HAVE_ATLAS']; - for %searches { - my $searchloc := $_; - my $test_ldd := pir::spawnw__IS('ldd ' ~ $searchloc); - if $test_ldd == 0 { - $found_blas := 1; - %PLA{'dynpmc_ldflags_list'}.push(%searches{$searchloc}[0]); - %PLA{'dynpmc_cflags_list'}.push(%searches{$searchloc}[1]); - return; - } - } + $found_blas := find_blas_linux(%PLA); } else { pir::say("Only Linux is currently supported"); @@ -111,6 +103,29 @@ sub find_blas(%PLA) { } } +sub find_blas_linux(%PLA) { + my $found_blas := 0; + my %searches; + # TODO: We should search in /usr/lib and /usr/local/lib for each + %searches{'/usr/lib/libblas-3.so'} := ['-lblas-3', '-D_PLA_HAVE_BLAS', 0]; + %searches{'/usr/lib/libblas.so'} := ['-lblas', '-D_PLA_HAVE_ATLAS', 1]; + %searches{'/usr/lib/atlas/libcblas.so'} := ['-L/usr/lib/atlas -lcblas', '-D_PLA_HAVE_ATLAS', 1]; + for %searches { + my $searchloc := $_; + my $test_ldd := pir::spawnw__IS('ldd ' ~ $searchloc); + if $test_ldd == 0 { + $found_blas := 1; + my @options := %searches{$searchloc}; + %PLA{'dynpmc_ldflags_list'}.push(@options[0]); + %PLA{'dynpmc_cflags_list'}.push(@options[1]); + %PLA{'need_cblas_h'} := @options[2]; + pir::say("=== PLA: Using BLAS library $searchloc"); + return $found_blas; + } + } + return $found_blas; +} + sub find_lapack(%PLA) { my %config := get_config(); my $osname := %config{'osname'}; @@ -236,4 +251,10 @@ sub setup_nqp_bootstrapper(%PLA) { %PLA{'inst_lib'}.push('pla_nqp.pbc'); } +sub setup_docs(%PLA) { + %PLA{'html_pod'}{'docs/nummatrix2d.html'} := 'src/pmc/nummatrix2d.pmc'; + %PLA{'html_pod'}{'docs/pmcmatrix2d.html'} := 'src/pmc/pmcmatrix2d.pmc'; + %PLA{'html_pod'}{'docs/complexmatrix2d.html'} := 'src/pmc/complexmatrix2d.pmc'; +} + diff --git a/src/include/pla_blas.h b/src/include/pla_blas.h index 36c359e..e1a8813 100644 --- a/src/include/pla_blas.h +++ b/src/include/pla_blas.h @@ -5,60 +5,67 @@ need to manually write bindings for that */ #ifdef _PLA_HAVE_CBLAS_H # include -# define dgemm cblas_dgemm -# define zgemm cblas_zgemm +# define PLA_HAVE_CBLAS +# define IS_TRANSPOSED_BLAS(flags) (IS_TRANSPOSED(flags) ? CblasTrans : CblasNoTrans) #else +# define IS_TRANSPOSED_BLAS(flags) (IS_TRANSPOSED(flags) ? "T" : "N") +# ifdef PLA_HAVE_CBLAS +# undef PLA_HAVE_CBLAS +# endif /* Define manual mappings to the BLAS library here, and map them to the same function names that would be used by ATLAS/cblas.h */ -/* Using the same names for these enums as cblas.h does, so that we can be - transparent between the two */ -typedef enum CBLAS_ORDER_ { - CblasRowMajor = 101, - CblasColMajor=102 -} CBLAS_ORDER; - -typedef enum CBLAS_TRANSPOSE_ { - CblasNoTrans = 111, - CblasTrans = 112, - CblasConjTrans = 113 -} CBLAS_TRANSPOSE; - extern void dgemm_( - const CBLAS_ORDER Order, - const CBLAS_TRANSPOSE TransA, - const CBLAS_TRANSPOSE TransB, - const int M, - const int N, - const int K, - const double alpha, + const void *TransA, + const void *TransB, + const void *M, + const void *N, + const void *K, + const double *alpha, const double *A, - const int lda, + const void *lda, const double *B, - const int ldb, - const double beta, + const void *ldb, + const double *beta, double *C, - const int ldc + const void *ldc ); -#define dgemm dgemm_ extern void zgemm_( - const CBLAS_ORDER Order, - const CBLAS_TRANSPOSE TransA, - const CBLAS_TRANSPOSE TransB, - const int M, - const int N, - const int K, + const void *TransA, + const void *TransB, + const void *M, + const void *N, + const void *K, const void *alpha, const void *A, - const int lda, + const void *lda, const void *B, - const int ldb, + const void *ldb, const void *beta, void *C, - const int ldc + const void *ldc ); # endif /* _PLA_HAVE_CBLAS_H */ + +void call_dgemm( + FLOATVAL alpha, + INTVAL flags_A, FLOATVAL * A, INTVAL rows_A, INTVAL cols_A, + INTVAL flags_B, FLOATVAL * B, INTVAL cols_B, + FLOATVAL beta, + FLOATVAL * C +); + +void call_zgemm( + FLOATVAL alpha_r, FLOATVAL alpha_i, + INTVAL flags_a, FLOATVAL * A, INTVAL rows_a, INTVAL cols_a, + INTVAL flags_b, FLOATVAL * B, INTVAL cols_b, + FLOATVAL beta_r, FLOATVAL beta_i, + FLOATVAL * C +); + + + #endif /* _PLA_BLAS_H_ */ diff --git a/src/include/pla_matrix_types.h b/src/include/pla_matrix_types.h index 2089b06..12b08c7 100644 --- a/src/include/pla_matrix_types.h +++ b/src/include/pla_matrix_types.h @@ -5,6 +5,20 @@ extern INTVAL __PLA_Type_NumMatrix2D; extern INTVAL __PLA_Type_ComplexMatrix2D; extern INTVAL __PLA_Type_PMCMatrix2D; +#define ALLOCATE_STORAGE_NumMatrix2D(s) \ + (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL)) +#define DECLATTRS_NumMatrix2D(p, a) Parrot_NumMatrix2D_attributes * const (a) = \ + (Parrot_NumMatrix2D_attributes *)((p)->data) + +#define ALLOCATE_STORAGE_ComplexMatrix2D(s) \ + (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL) * 2) +#define DECLATTRS_ComplexMatrix2D(p, a) Parrot_ComplexMatrix2D_attributes * const (a) = \ + (Parrot_ComplexMatrix2D_attributes *)((p)->data) + +#define ALLOCATE_STORAGE_PMCMatrix2D(s) (PMC **)mem_sys_allocate_zeroed(s * sizeof (PMC *)) +#define DECLATTRS_PMCMatrix2D(p, a) Parrot_PMCMatrix2D_attributes * const (a) = \ + (Parrot_PMCMatrix2D_attributes *)((p)->data) + #define SWAP_XY(a) do { \ const INTVAL __temp_val = a->rows; \ a->rows = a->cols; \ @@ -164,6 +178,5 @@ do { \ #define IS_DIAGONAL(flags) ((((flags) & (FLAG_DIAGONAL)) == FLAG_DIAGONAL)) #define IS_TRIDIAGONAL(flags) (((flags) & (FLAG_TRIDIAGONAL))) #define IS_TRANSPOSED(flags) (((flags) & (FLAG_TRANSPOSED))) -#define IS_TRANSPOSED_BLAS(flags) (IS_TRANSPOSED(flags) ? CblasTrans : CblasNoTrans) #endif /* _PLA_MATRIX_TYPES_H_ */ diff --git a/src/lib/pla_blas.c b/src/lib/pla_blas.c index 64d2387..0d69fa8 100644 --- a/src/lib/pla_blas.c +++ b/src/lib/pla_blas.c @@ -3,81 +3,100 @@ /* Wrapper to call the dgemm function from BLAS with PMC arguments. Assumes A, B, and C are all NumMatrix2D. */ void -call_dgemm(PARROT_INTERP, FLOATVAL alpha, PMC * A, PMC *B, FLOATVAL beta, PMC *C) +call_dgemm(FLOATVAL alpha, + INTVAL flags_a, FLOATVAL * A, INTVAL rows_a, INTVAL cols_a, + INTVAL flags_b, FLOATVAL * B, INTVAL cols_b, + FLOATVAL beta, FLOATVAL * C) { - DECLATTRS(A, attrs_a); - DECLATTRS(B, attrs_b); - DECLATTRS(C, attrs_c); - const INTVAL M = attrs_a->rows; - const INTVAL N = attrs_b->cols; - const INTVAL K = attrs_a->cols; - if (attrs_c->rows != M) - Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, - "PLA DGEMM: A, C indices do not match in gemm"); - if (attrs_c->cols != N) - Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, - "PLA DGEMM: B, C indices do not match in gemm"); - if (attrs_b->rows != K) - Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, - "PLA DGEMM: A, B indeces do not match in gemm"); - dgemm(CblasRowMajor, - IS_TRANSPOSED_BLAS(attrs_a->flags), - IS_TRANSPOSED_BLAS(attrs_b->flags), + const INTVAL M = rows_a; + const INTVAL N = cols_b; + const INTVAL K = cols_a; +#ifdef PLA_HAVE_CBLAS + cblas_dgemm(CblasRowMajor, + IS_TRANSPOSED_BLAS(flags_a), + IS_TRANSPOSED_BLAS(flags_b), M, N, K, alpha, - attrs_a->storage, + A, M, - attrs_b->storage, + B, N, beta, - attrs_c->storage, + C, M ); +#else + dgemm_( + IS_TRANSPOSED_BLAS(flags_a), + IS_TRANSPOSED_BLAS(flags_b), + &M, + &N, + &K, + &alpha, + A, + &M, + B, + &N, + &beta, + C, + &M + ); +#endif } /* Wrapper to call the zdgemm function from BLAS with PMC arguments. Assumes A, B, and C are all ComplexMatrix2D. */ -static void -call_zgemm(PARROT_INTERP, FLOATVAL alpha_r, FLOATVAL alpha_i, PMC * A, PMC *B, - FLOATVAL beta_r, FLOATVAL beta_i, PMC *C) +void +call_zgemm(FLOATVAL alpha_r, FLOATVAL alpha_i, + INTVAL flags_a, FLOATVAL * A, INTVAL rows_a, INTVAL cols_a, + INTVAL flags_b, FLOATVAL * B, INTVAL cols_b, + FLOATVAL beta_r, FLOATVAL beta_i, FLOATVAL * C) { - DECLATTRS(A, attrs_a); - DECLATTRS(B, attrs_b); - DECLATTRS(C, attrs_c); - const INTVAL M = attrs_a->rows; - const INTVAL N = attrs_b->cols; - const INTVAL K = attrs_a->cols; + const INTVAL M = rows_a; + const INTVAL N = cols_b; + const INTVAL K = cols_a; FLOATVAL alpha_p[2]; FLOATVAL beta_p[2]; - if (attrs_c->rows != M) - Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, - "PLA ZGEMM: A, C indices do not match in gemm"); - if (attrs_c->cols != N) - Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, - "PLA ZGEMM: B, C indices do not match in gemm"); - if (attrs_b->rows != K) - Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, - "PLA ZGEMM: A, B indeces do not match in gemm"); + alpha_p[0] = alpha_r; alpha_p[1] = alpha_i; beta_p[0] = beta_r; beta_p[1] = beta_i; - zgemm(CblasRowMajor, - IS_TRANSPOSED_BLAS(attrs_a->flags), - IS_TRANSPOSED_BLAS(attrs_b->flags), + +#ifdef PLA_HAVE_CBLAS + cblas_zgemm(CblasRowMajor, + IS_TRANSPOSED_BLAS(flags_a), + IS_TRANSPOSED_BLAS(flags_b), M, N, K, alpha_p, - attrs_a->storage, + A, M, - attrs_b->storage, + B, N, beta_p, - attrs_c->storage, + C, M ); +#else + zgemm_( + IS_TRANSPOSED_BLAS(flags_a), + IS_TRANSPOSED_BLAS(flags_b), + &M, + &N, + &K, + alpha_p, + A, + &M, + B, + &N, + beta_p, + C, + &M + ); +#endif } diff --git a/src/pmc/complexmatrix2d.pmc b/src/pmc/complexmatrix2d.pmc index 16c912b..6c266cc 100644 --- a/src/pmc/complexmatrix2d.pmc +++ b/src/pmc/complexmatrix2d.pmc @@ -1,9 +1,19 @@ #include "pla.h" - -#define ALLOCATE_STORAGE(s) (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL) * 2) #define PLATYPENAME "ComplexMatrix2D" -#define DECLATTRS(p, a) Parrot_ComplexMatrix2D_attributes * const (a) = \ - (Parrot_ComplexMatrix2D_attributes *)((p)->data) + +/* + +=head1 ComplexMatrix2D + +=head2 Description + +ComplexMatrix2D is a 2-dimensional complex-valued matrix type for the Parrot +Virtual Machine. It supports a variety of element-wise and matrix operations, +including bindings to the BLAS and LAPACK libraries. + +=cut + +*/ INTVAL __PLA_Type_ComplexMatrix2D; @@ -13,12 +23,12 @@ INTVAL __PLA_Type_ComplexMatrix2D; static void normalize_lazy_transpose(PARROT_INTERP, PMC * self) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); if (IS_TRANSPOSED(attrs->flags)) { const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL size = rows_size * cols_size; - FLOATVAL * const new_s = ALLOCATE_STORAGE(size); + FLOATVAL * const new_s = ALLOCATE_STORAGE_ComplexMatrix2D(size); FLOATVAL * const old_s = attrs->storage; INTVAL i, j; @@ -49,7 +59,8 @@ convert_to_ComplexMatrix2D(PARROT_INTERP, PMC * p, INTVAL forcecopy) Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, PLATYPENAME ": cannot convert unknown PMC type"); else { - PMC * const meth = VTABLE_find_method(interp, p, CONST_STRING(interp, "convert_to_complex_matrix")); + STRING * const converter = CONST_STRING(interp, "convert_to_complex_matrix"); + PMC * const meth = VTABLE_find_method(interp, p, converter); PMC * dest = PMCNULL; Parrot_ext_call(interp, meth, "Pi->P", p, &dest); return dest; @@ -71,7 +82,7 @@ convert_to_ComplexMatrix2D(PARROT_INTERP, PMC * p, INTVAL forcecopy) static void resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); /* Store the old values */ const INTVAL old_rows = attrs->rows; const INTVAL old_cols = attrs->cols; @@ -82,7 +93,7 @@ resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) const INTVAL new_rows = INDEX_MAX(old_rows, row + 1); const INTVAL new_cols = INDEX_MAX(old_cols, col + 1); const INTVAL newsize = new_rows * new_cols; - FLOATVAL * new_s = ALLOCATE_STORAGE(newsize); + FLOATVAL * new_s = ALLOCATE_STORAGE_ComplexMatrix2D(newsize); INTVAL i, j; @@ -105,7 +116,7 @@ static void init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, INTVAL cols_size, PMC * values) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); INTVAL num = 0; const INTVAL init_elems = VTABLE_elements(interp, values); const INTVAL total_elems = rows_size * cols_size; @@ -120,7 +131,7 @@ init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, static PMC * get_complex_pmc_at_xy(PARROT_INTERP, PMC *self, INTVAL rows, INTVAL cols) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; PMC * newcomplex; @@ -143,7 +154,7 @@ static void set_complex_pmc_at_xy(PARROT_INTERP, PMC * self, PMC * value, INTVAL row, INTVAL col) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); INTVAL rows_size = attrs->rows; INTVAL cols_size = attrs->cols; FLOATVAL real, imag; @@ -165,7 +176,7 @@ static void set_scalar_at_xy(PARROT_INTERP, PMC * self, FLOATVAL value, INTVAL row, INTVAL col) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); INTVAL rows_size = attrs->rows; INTVAL cols_size = attrs->cols; const INTVAL flags = attrs->flags; @@ -185,7 +196,7 @@ set_scalar_at_xy(PARROT_INTERP, PMC * self, FLOATVAL value, INTVAL row, static void add_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; FLOATVAL * const s = attrs->storage; @@ -202,7 +213,7 @@ add_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) static void add_scalar_complex(PARROT_INTERP, PMC * self, PMC *v, INTVAL sub) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; FLOATVAL * const s = attrs->storage; @@ -232,7 +243,7 @@ add_scalar_complex(PARROT_INTERP, PMC * self, PMC *v, INTVAL sub) static void multiply_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; FLOATVAL * const s = attrs->storage; @@ -250,7 +261,7 @@ multiply_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) static void multiply_scalar_complex(PARROT_INTERP, PMC * self, PMC * v) { - DECLATTRS(self, attrs); + DECLATTRS_ComplexMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; FLOATVAL * const s = attrs->storage; @@ -270,13 +281,34 @@ multiply_scalar_complex(PARROT_INTERP, PMC * self, PMC * v) } } +static void +multiply_matrices(PARROT_INTERP, FLOATVAL alpha_r, FLOATVAL alpha_i, PMC * A, + PMC * B, FLOATVAL beta_r, FLOATVAL beta_i, PMC * C) +{ + DECLATTRS_ComplexMatrix2D(A, attrs_a); + DECLATTRS_ComplexMatrix2D(B, attrs_b); + DECLATTRS_ComplexMatrix2D(C, attrs_c); + if (attrs_c->rows != attrs_a->rows) + Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": A, C indices do not match in gemm"); + if (attrs_c->cols != attrs_b->cols) + Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": B, C indices do not match in gemm"); + if (attrs_b->rows != attrs_a->cols) + Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": A, B indeces do not match in gemm"); + call_zgemm(alpha_r, alpha_i, attrs_a->flags, attrs_a->storage, attrs_a->rows, attrs_a->cols, + attrs_b->flags, attrs_b->storage, attrs_b->cols, beta_r, beta_i, attrs_c->storage); +} + /* item-by-item addition or subtraction A = A + B */ +// TODO: Support adding other types of matrices here, also. static void add_matrices(PARROT_INTERP, PMC * A, PMC * B, INTVAL sub) { - DECLATTRS(A, attrs_a); - DECLATTRS(B, attrs_b); + DECLATTRS_ComplexMatrix2D(A, attrs_a); + DECLATTRS_ComplexMatrix2D(B, attrs_b); const INTVAL rows = attrs_a->rows; const INTVAL cols = attrs_a->cols; FLOATVAL * const s_a = attrs_a->storage; @@ -321,16 +353,26 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix =head1 VTABLEs +=head2 System VTABLEs + =over 4 -=item* init +=item * init + +Initialize the new PMC + +=item * destroy + +Destroy the PMC and free all associated memory + +=back =cut */ VTABLE void init() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); attrs->storage = NULL; attrs->rows = 0; attrs->cols = 0; @@ -339,7 +381,7 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix } VTABLE void destroy() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; if (s) mem_sys_free(s); @@ -347,20 +389,50 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix /* -=item* get_number_keyed +=head2 Keyed lookup VTABLEs + +In each of these cases, the specified Key PMC must have exactly two elements +to specify a location in the matrix. + +Attempting to retrieve a value outside the boundaries of the matrix will throw +an OUT_OF_BOUNDS exception. + +=over 4 + +=item * get_number_keyed + +Get a floating point value at the location specified by the key. -=item* get_integer_keyed +The floating point value is calculated as the magnitude of the complex value +at the point in the matrix. -=item* get_string_keyed +=item * get_integer_keyed -=item* get_pmc_keyed +Get an integer value at the location specified by the key. + +The integer value is the magnitude of the complex value at the specified +location, cast from a floating point to an integer value. + +=item * get_string_keyed + +Get the string representation of the complex value at the point specified by +the key. + +The returned string has the same format as is used by the Parrot built-in +Complex PMC type. + +=item * get_pmc_keyed + +Get a Complex PMC from the point specified by the key. + +=back =cut */ VTABLE FLOATVAL get_number_keyed(PMC * key) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; INTVAL rows, cols; @@ -374,31 +446,13 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix return sqrt(real * real + imag + imag); } - VTABLE FLOATVAL get_number_keyed_int(INTVAL key) { - DECLATTRS(SELF, attrs); - const INTVAL totalsize = attrs->rows * attrs->cols; - FLOATVAL real, imag; - if (key >= totalsize || key < 0) { - Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, - PLATYPENAME ": Index out of bounds."); - } - real = attrs->storage[key * 2]; - imag = attrs->storage[key * 2 + 1]; - return sqrt(real * real + imag * imag); - } - VTABLE INTVAL get_integer_keyed(PMC * key) { const FLOATVAL f = VTABLE_get_number_keyed(INTERP, SELF, key); return (INTVAL)f; } - VTABLE INTVAL get_integer_keyed_int(INTVAL key) { - const FLOATVAL f = VTABLE_get_number_keyed_int(INTERP, SELF, key); - return (INTVAL)f; - } - VTABLE STRING * get_string_keyed(PMC * key) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; INTVAL rows, cols; @@ -414,8 +468,69 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix return Parrot_sprintf_c(INTERP, "%vg%+vgi", real, imag); } + VTABLE PMC * get_pmc_keyed(PMC * key) { + INTVAL rows, cols; + GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols); + return get_complex_pmc_at_xy(INTERP, SELF, rows, cols); + } + +/* + +=head2 Integer-Keyed Lookup VTABLES + +These VTABLEs treat the matrix, which is a contiguous region in memory, as a +linear array of values. The matrix data is stored by rows. + +These routines are used for low-level access. Attempting to access a value +outside the bounds of the matrix will throw an OUT_OF_BOUNDS exception. + +=over 4 + +=item * get_number_keyed_int + +Get a floating point number, the magnitude of the complex value, at the +specified location in the array + +=item * get_integer_keyed_int + +Get an integer, cast from the magnitude of the complex value, at the specified +location in the array + +=item * get_string_keyed_int + +Get a string representation of the value at the specified point. Stringification +is done the same as Parrot's Complex PMC type. + +=item * get_pmc_keyed_int + +Get a Complex PMC from the value at the specified point. + +=back + +=cut + +*/ + + VTABLE FLOATVAL get_number_keyed_int(INTVAL key) { + DECLATTRS_ComplexMatrix2D(SELF, attrs); + const INTVAL totalsize = attrs->rows * attrs->cols; + FLOATVAL real, imag; + if (key >= totalsize || key < 0) { + Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": Index out of bounds."); + } + real = attrs->storage[key * 2]; + imag = attrs->storage[key * 2 + 1]; + return sqrt(real * real + imag * imag); + } + + VTABLE INTVAL get_integer_keyed_int(INTVAL key) { + const FLOATVAL f = VTABLE_get_number_keyed_int(INTERP, SELF, key); + return (INTVAL)f; + } + VTABLE STRING * get_string_keyed_int(INTVAL key) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL real, imag; const INTVAL totalsize = attrs->rows * attrs->cols; if (key >= totalsize || totalsize < 0) { @@ -427,14 +542,8 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix return Parrot_sprintf_c(INTERP, "%vg%+vgi", real, imag); } - VTABLE PMC * get_pmc_keyed(PMC * key) { - INTVAL rows, cols; - GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols); - return get_complex_pmc_at_xy(INTERP, SELF, rows, cols); - } - VTABLE PMC * get_pmc_keyed_int(INTVAL key) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); INTVAL row, col; intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col); return get_complex_pmc_at_xy(INTERP, SELF, row, col); @@ -442,9 +551,38 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix /* -=item* set_pmc_keyed +=head2 Keyed Setter VTABLES + +These VTABLEs insert new values into the matrix at a point specified by the +Key PMC. The Key PMC must have exactly two elements. If the matrix is not large +enough to accomodate the specified location, it will be grown with zero-padding +so that it is at least large enough to hold the specified point and all existing +data. + +=over 4 + +=item * set_pmc_keyed + +Set the value at the specified location to the complex value represented by the +given PMC. Different PMC types have different behaviors in this operation. + +=item * set_number_keyed + +Set the value at the specified location to have a real value given by the +number, and 0 for the complex value. + +=item * set_string_keyed -=item* set_string_keyed +Set the value at the specified location to the value represented by the string. +The string is parsed by creating a temporary Complex PMC type, so all the same +rules for that type apply to the input string format. + +=item * set_integer_keyed + +Set the value at the specified location to have a real value given by the +integer, and 0 for the complex value. + +=back =cut @@ -456,26 +594,12 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix set_complex_pmc_at_xy(INTERP, SELF, value, rows, cols); } - VTABLE void set_pmc_keyed_int(INTVAL key, PMC * value) { - DECLATTRS(SELF, attrs); - INTVAL row, col; - intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col); - set_complex_pmc_at_xy(INTERP, SELF, value, row, col); - } - VTABLE void set_number_keyed(PMC * key, FLOATVAL value) { INTVAL rows, cols; GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols); set_scalar_at_xy(INTERP, SELF, value, rows, cols); } - VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL value) { - DECLATTRS(SELF, attrs); - INTVAL row, col; - intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col); - set_scalar_at_xy(INTERP, SELF, value, row, col); - } - VTABLE void set_string_keyed(PMC * key, STRING * value) { PMC * const item = Parrot_pmc_new(INTERP, enum_class_Complex); VTABLE_set_string_native(INTERP, item, value); @@ -486,12 +610,65 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix VTABLE_set_number_keyed(INTERP, SELF, key, (FLOATVAL)value); } +/* + +=head2 Integer-Keyed Setter VTABLEs + +These VTABLEs treat the matrix as a linear array in memory and allow fast +lookup based on the integer offset of values in the array. These are low-level +routines and are not intended for general use. + +Unlike the PMC-keyed VTABLEs, these routines will not automatically grow the +matrix if an index is provided which is outside the boundaries of the matrix. +In that case, an OUT_OF_BOUNDS exception will be thrown. + +=over 4 + +=item * set_pmc_keyed_int + +Set a PMC at the specified location. The rules for extracting a complex value +out of the input PMC are the same as used for set_pmc_keyed() + +=item * set_number_keyed_int + +Set the complex value at the specified location to the given real value. The +complex part of the value will be 0. + +=item * set_integer_keyed_int + +Set the complex value at the specified location to the given real value. The +complex part of the value will be 0 + +=item * set_string_keyed_int + +Convert the string to a Complex PMC, and set that at the specified location. + +=back + +=cut + +*/ + + VTABLE void set_pmc_keyed_int(INTVAL key, PMC * value) { + DECLATTRS_ComplexMatrix2D(SELF, attrs); + INTVAL row, col; + intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col); + set_complex_pmc_at_xy(INTERP, SELF, value, row, col); + } + + VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL value) { + DECLATTRS_ComplexMatrix2D(SELF, attrs); + INTVAL row, col; + intkey_to_coords(INTERP, attrs->rows, attrs->cols, key, &row, &col); + set_scalar_at_xy(INTERP, SELF, value, row, col); + } + VTABLE void set_integer_keyed_int(INTVAL key, INTVAL value) { VTABLE_set_number_keyed_int(INTERP, SELF, key, (FLOATVAL)value); } VTABLE void set_string_keyed_int(INTVAL key, STRING * value) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); INTVAL row, col; PMC * const item = Parrot_pmc_new(INTERP, enum_class_Complex); VTABLE_set_string_native(INTERP, item, value); @@ -500,14 +677,34 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix /* -=item* get_string +=head2 Miscellaneous VTABLEs + +=over 4 + +=item * get_string + +Get a string representation of the matrix, suitable for printing to the console + +=item * get_attr_string + +Get a named attribute. The name can be one of "rows", "cols", or "size". + +=item * clone + +Clone the matrix + +=item * is_equal + +Determine if two matrices are equal in size and composition. + +=back =cut */ VTABLE STRING *get_string() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); INTVAL rows, cols; PMC * const builder = Parrot_pmc_new(INTERP, enum_class_StringBuilder); STRING * const newline = Parrot_str_new(INTERP, "\n", 1); @@ -541,7 +738,7 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix } VTABLE PMC * get_attr_str(STRING * idx) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "rows"))) { PMC * const rows = Parrot_pmc_new(INTERP, enum_class_Integer); VTABLE_set_integer_native(INTERP, rows, attrs->rows); @@ -562,14 +759,14 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix VTABLE PMC * clone() { PMC * const c = Parrot_pmc_new(INTERP, SELF->vtable->base_type); - DECLATTRS(SELF, old_atts); - DECLATTRS(c, new_atts); + DECLATTRS_ComplexMatrix2D(SELF, old_atts); + DECLATTRS_ComplexMatrix2D(c, new_atts); INTVAL rows, cols; INTVAL const rows_size = old_atts->rows; INTVAL const cols_size = old_atts->cols; INTVAL const newsize = rows_size * cols_size; FLOATVAL * const old_s = old_atts->storage; - FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize); + FLOATVAL * const new_s = ALLOCATE_STORAGE_ComplexMatrix2D(newsize); for (rows = 0; rows < rows_size; ++rows) { for (cols = 0; cols < cols_size; ++cols) { R_ITEM_XY_ROWMAJOR(new_s, rows_size, cols_size, rows, cols) = @@ -587,8 +784,8 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix VTABLE INTVAL is_equal(PMC * other) { if (other->vtable->base_type == SELF->vtable->base_type) { - DECLATTRS(SELF, self_attrs); - DECLATTRS(other, other_attrs); + DECLATTRS_ComplexMatrix2D(SELF, self_attrs); + DECLATTRS_ComplexMatrix2D(other, other_attrs); const INTVAL self_rows = self_attrs->rows; const INTVAL self_cols = self_attrs->cols; const INTVAL self_flags = self_attrs->flags; @@ -623,16 +820,27 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix /* -=item* freeze +=head2 Serialization/Deserialization VTABLEs + +=over 4 + +=item * freeze + +Freeze the PMC for serialization to a string suitable for long-term storage in +a file. -=item* thaw +=item * thaw + +Thaw a serialized PMC + +=back =cut */ VTABLE void freeze(PMC *info) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); INTVAL const rows = attrs->rows; INTVAL const cols = attrs->cols; INTVAL const flags = attrs->flags; @@ -652,7 +860,7 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix } VTABLE void thaw(PMC *info) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); INTVAL const rows = VTABLE_shift_integer(INTERP, info); INTVAL const cols = VTABLE_shift_integer(INTERP, info); INTVAL const flags = VTABLE_shift_integer(INTERP, info); @@ -677,23 +885,33 @@ pmclass ComplexMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix /* -=item* add(ComplexMatrix2D) +=head2 Addition VTABLEs + +=over 4 + +=item * add(ComplexMatrix2D) Add two matrices together, element-by-element. -=item* add(DEFAULT) +=item * add(Complex) + +Add a complex value to each element of the matrix + +=item * add(DEFAULT) Take the number representation of the argument PMC and add it to every element in the matrix. -=item* add_float +=item * add_float Add the float value to every element in the matrix. -=item* add_int +=item * add_int Add the integer value to every element in the matrix. +=back + =cut */ @@ -731,13 +949,31 @@ Add the integer value to every element in the matrix. /* -=item* i_add(ComplexMatrix2D) +=head2 In-Place Addition VTABLEs + +=over 4 + +=item * i_add(ComplexMatrix2D) + +Add a matrix into SELF + +=item * i_add(Complex) + +Add a Complex value to every element of SELF + +=item * i_add(DEFAULT) + +Add the numeric value of the given PMC to every element of SELF + +=item * i_add_int -=item* i_add(DEFAULT) +Add the integer value to every element of SELF -=item* i_add_int +=item * i_add_float -=item* i_add_float +Add the floating-point value to every element of SELF + +=back =cut @@ -766,22 +1002,34 @@ Add the integer value to every element in the matrix. /* -=item* subtract(ComplexMatrix2D) +=head2 Subtraction VTABLEs -Add two matrices together, element-by-element. +=over 4 -=item* subtract(DEFAULT) +=item * subtract(ComplexMatrix2D) -Take the number representation of the argument PMC and add it to every element -in the matrix. +Perform the matrix operation -=item* subtract_float + C = SELF - A -Add the float value to every element in the matrix. +=item * subtract(Complex) -=item* subtract_int +Subtract a complex value from each element of the matrix -Add the integer value to every element in the matrix. +=item * subtract(DEFAULT) + +Take the number representation of the argument PMC and subtract it from every +element in the matrix. + +=item * subtract_float + +Subtract the float value from every element in the matrix. + +=item * subtract_int + +Subtract the integer value from every element in the matrix. + +=back =cut @@ -821,13 +1069,35 @@ Add the integer value to every element in the matrix. /* -=item* i_subtract(ComplexMatrix2D) +=head2 In-Place Subtraction VTABLEs + +=over 4 + +=item * i_subtract(ComplexMatrix2D) -=item* i_subtract(DEFAULT) +Perform the matrix operation -=item* i_subtract_int + SELF = SELF - A -=item* i_subtract_float +=item * i_subtract(Complex) + +Substract a complex value from each element of the matrix + +=item * i_subtract(DEFAULT) + +Subtract the numeric value of the PMC from the real part of each element in the +matrix + +=item * i_subtract_int + +Subtract the integer value from the real part of each element in the matrix + +=item * i_subtract_float + +Subtract the floating point value from the real part of each element in the +matrix + +=back =cut @@ -856,17 +1126,41 @@ Add the integer value to every element in the matrix. /* -=item* multiply(ComplexMatrix2D) +=head2 Multiplication VTABLEs -=item* multiply(DEFAULT) +=over 4 + +=item * multiply(ComplexMatrix2D) + +Perform the matrix operation + + C = SELF × A + +=item * multiply(Complex) + +Multiply each element in the matrix by the complex value + +=item * multiply(DEFAULT) + +Multiply each element in the matrix by the numeric value of the PMC + +=item * multiply_int + +Multiply each element in the matrix by the integer + +=item * multiply_float + +Multiply each element in the matrix by the floating point number + +=back =cut */ MULTI PMC *multiply(ComplexMatrix2D *value, PMC *dest) { - DECLATTRS(SELF, selfattr); - DECLATTRS(value, valattr); + DECLATTRS_ComplexMatrix2D(SELF, selfattr); + DECLATTRS_ComplexMatrix2D(value, valattr); const INTVAL new_rows = selfattr->rows; const INTVAL new_cols = valattr->cols; @@ -877,7 +1171,7 @@ Add the integer value to every element in the matrix. dest = Parrot_pmc_new(INTERP, VTABLE_type(INTERP, SELF)); resize_matrix(INTERP, dest, new_rows - 1, new_cols - 1); - call_zgemm(INTERP, 1.0, 0.0, SELF, value, 0.0, 0.0, dest); + multiply_matrices(INTERP, 1.0, 0.0, SELF, value, 0.0, 0.0, dest); return dest; } @@ -894,16 +1188,6 @@ Add the integer value to every element in the matrix. return dest; } -/* - -=item* multiply_int - -=item* multiply_float - -=cut - -*/ - VTABLE PMC *multiply_int(INTVAL v, PMC * dest) { dest = VTABLE_clone(INTERP, SELF); multiply_scalar_float(INTERP, dest, (FLOATVAL)v); @@ -918,9 +1202,33 @@ Add the integer value to every element in the matrix. /* -=item* i_multiply(NumMatrix2D) +=head2 In-Place Multiplication VTABLEs + +=over 4 + +=item * i_multiply(ComplexMatrix2D) + +Perform the matrix operation + + SELF = SELF × A + +=item * i_multiply(Complex) -=item* i_multiply(DEFAULT) +Multiply every element in SELF by the complex value + +=item * i_multiply(DEFAULT) + +Multiply every element in SELF by the numeric value of the PMC + +=item * i_multiply_int + +Multiply every element in SELF by the integer + +=item * i_multiply_float + +Multiply every element in SELF by the floating point number + +=back =cut @@ -928,7 +1236,7 @@ Add the integer value to every element in the matrix. MULTI void i_multiply(ComplexMatrix2D* value) { PMC * const temp = VTABLE_clone(INTERP, SELF); - call_zgemm(INTERP, 1.0, 0.0, temp, value, 0.0, 0.0, SELF); + multiply_matrices(INTERP, 1.0, 0.0, temp, value, 0.0, 0.0, SELF); } MULTI void i_multiply(Complex *value) { @@ -940,16 +1248,6 @@ Add the integer value to every element in the matrix. multiply_scalar_float(INTERP, SELF, v); } -/* - -=item i_multiply_int - -=item i_multiply_float - -=cut - -*/ - VTABLE void i_multiply_int(INTVAL v) { multiply_scalar_float(INTERP, SELF, (FLOATVAL)v); } @@ -961,13 +1259,16 @@ Add the integer value to every element in the matrix. /* -=back - =head1 METHODS =over 4 -=item resize() +=item * resize() + +Resize the matrix to include at least the specified number of rows and columns. + +Resizing the matrix never causes the matrix to shrink. If you need a subset of +the matrix, use get_block instead. =cut @@ -979,10 +1280,11 @@ Add the integer value to every element in the matrix. /* -=item fill() +=item * fill() Fill the matrix with a single value. if sizes are provided, fill to those -sizes, growing the matrix if needed. +sizes, growing the matrix if needed. Elements outside the specified area are +unaffected. Calling fill() never causes the matrix to shrink. =cut @@ -993,7 +1295,7 @@ sizes, growing the matrix if needed. INTVAL cols_size :optional, INTVAL has_cols_size :opt_flag ) { /* TODO: Value here is going to be a Complex PMC. Handle that */ - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL * s = attrs->storage; INTVAL const curr_rows_size = attrs->rows; INTVAL const curr_cols_size = attrs->cols; @@ -1019,7 +1321,9 @@ sizes, growing the matrix if needed. /* -=item item_at() +=item * item_at() + +Return a single Complex PMC from the item at the specified coordinates =cut @@ -1027,7 +1331,7 @@ sizes, growing the matrix if needed. METHOD item_at(INTVAL row, INTVAL col, PMC * value :optional, INTVAL has_value :opt_flag) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; if (row < 0 || col < 0 || row >= rows || col >= cols) { @@ -1046,16 +1350,19 @@ sizes, growing the matrix if needed. /* -=item transpose() +=item * transpose() -Transposes the matrix. +Transposes the matrix lazily. This operation is O(1). Some operations, such as +mathematical operations do not work on a matrix which has been lazily +transposed, so those operations will force the matrix memory to be eagerly +transposed. =cut */ METHOD transpose() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); INTVAL tmp = 0; INTVAL transposed = IS_TRANSPOSED(attrs->flags); @@ -1071,21 +1378,22 @@ Transposes the matrix. /* -=item mem_transpose() +=item * mem_transpose() -Transposes the actual data storage of the matrix. More expensive up-front -than the transpose() method. +Transposes the actual data storage of the matrix. More expensive O(n) up-front +than the transpose() method, but the resulting memory structure is more suitable +for use in certain mathematical operations. =cut */ METHOD mem_transpose() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL newsize = rows_size * cols_size; - FLOATVAL * new_s = ALLOCATE_STORAGE(newsize); + FLOATVAL * new_s = ALLOCATE_STORAGE_ComplexMatrix2D(newsize); FLOATVAL * old_s = attrs->storage; INTVAL i, j; @@ -1105,14 +1413,16 @@ than the transpose() method. /* -=item conjugate +=item * conjugate + +Convert the matrix to the complex conjugate of itself. =cut */ METHOD conjugate() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL newsize = rows_size * cols_size; @@ -1129,17 +1439,22 @@ than the transpose() method. /* -=item iterate_function_inplace() +=item * iterate_function_inplace() -Calls a function for every element in the array, replacing the current +Calls a function for every element in the matrix, replacing the current value with the return value of the called function. +=item * iterate_function_external() + +Calls a function for every element in the matrix, creating a new matrix with +the returned values of the called function. + =cut */ METHOD iterate_function_inplace(PMC * func, PMC * args :slurpy) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; INTVAL i, j; @@ -1155,7 +1470,7 @@ value with the return value of the called function. } METHOD iterate_function_external(PMC * func, PMC * args :slurpy) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; @@ -1175,6 +1490,20 @@ value with the return value of the called function. RETURN(PMC * new_matrix); } +/* + +=item * initialize_from_array + +Initialize the matrix using a list of values from an array. + +=item * initialize_from_args + +Initialize the matrix using values from a variadic (slurpy) argument list. + +=cut + +*/ + METHOD initialize_from_array(INTVAL rows_size, INTVAL cols_size, PMC *values) { init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values); @@ -1185,9 +1514,24 @@ value with the return value of the called function. init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values); } +/* + +=item * get_block + +Get a specified sub-block of the matrix. If the bounds of the sub-block are +outside the bounds of the matrix, an OUT_OF_BOUNDS exception is thrown. + +=item * set_block + +Set a block in the matrix, growing it if needed. + +=cut + +*/ + METHOD get_block(INTVAL rows_idx, INTVAL cols_idx, INTVAL rows_size, INTVAL cols_size) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1221,8 +1565,8 @@ value with the return value of the called function. } METHOD set_block(INTVAL rows_idx, INTVAL cols_idx, PMC * blck) { - DECLATTRS(SELF, self_attrs); - DECLATTRS(blck, blck_attrs); + DECLATTRS_ComplexMatrix2D(SELF, self_attrs); + DECLATTRS_ComplexMatrix2D(blck, blck_attrs); FLOATVAL * self_s = self_attrs->storage; FLOATVAL * const blck_s = blck_attrs->storage; INTVAL self_rows = self_attrs->rows; @@ -1259,12 +1603,16 @@ value with the return value of the called function. /* -=item* gemm +=item * gemm Calculates the matrix equation: Z = aAB + bC +The matrices must all be ComplexMatrix2D, or must be convertable to it. The +matrix SELF is not used in the calculation, but the result matrix will have the +same type as SELF. + =cut */ @@ -1276,30 +1624,34 @@ Calculates the matrix equation: C = convert_to_ComplexMatrix2D(interp, C, 1); get_complex_value_from_pmc(interp, alpha, &alpha_r, &alpha_i); get_complex_value_from_pmc(interp, beta, &beta_r, &beta_i); - call_zgemm(INTERP, alpha_r, alpha_i, A, B, beta_r, beta_i, C); + multiply_matrices(INTERP, alpha_r, alpha_i, A, B, beta_r, beta_i, C); RETURN(PMC* C); } /* -=item row_combine(srcidx, destidx, gain) +=item * row_combine(srcidx, destidx, gain) -add a multiple of the source row to the destination row. +add a multiple of the source row to the destination row. If either of the row +indices are outside the bounds of the matrix, an OUT_OF_BOUNDS exception is +thrown. -=item row_scale(idx, gain) +=item * row_scale(idx, gain) -Multiply all elements in the row by a gain factor. +Multiply all elements in the row by a gain factor. If the row index is outside +the bounds of the matrix and OUT_OF_BOUNDS exception is thrown. -=item row_swap(idx_a, idx_b) +=item * row_swap(idx_a, idx_b) -Swap two rows +Swap two rows. If either of the row indices are outside the bounds of the +matrix, an OUT_OF_BOUNDS exception is thrown. =cut */ METHOD row_combine(INTVAL srcidx, INTVAL destidx, FLOATVAL gain) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1317,7 +1669,7 @@ Swap two rows } METHOD row_scale(INTVAL idx, FLOATVAL gain) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1333,7 +1685,7 @@ Swap two rows } METHOD row_swap(INTVAL idx_a, INTVAL idx_b) { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1356,24 +1708,27 @@ Swap two rows /* -=item convert_to_number_matrix +=item * convert_to_number_matrix -Get a NumMatrix2D from the current matrix +Get a NumMatrix2D from the current matrix. If the matrix is already a +NumMatrix2D, return a clone. -=item convert_to_complex_matrix +=item * convert_to_complex_matrix -Get a ComplexMatrix2D from the current matrix +Get a ComplexMatrix2D from the current matrix. If the matrix is already a +ComplexMatrix2D, return a clone. -=item convert_to_pmc_matrix +=item * convert_to_pmc_matrix -Get a PMCMatrix2D from the current matrix +Get a PMCMatrix2D from the current matrix. If the matrix is already a +PMCMatrix2D, return a clone. =cut */ METHOD convert_to_number_matrix() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_NumMatrix2D); const INTVAL totalsize = attrs->rows * attrs->cols; PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize")); @@ -1394,7 +1749,7 @@ Get a PMCMatrix2D from the current matrix } METHOD convert_to_pmc_matrix() { - DECLATTRS(SELF, attrs); + DECLATTRS_ComplexMatrix2D(SELF, attrs); PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_PMCMatrix2D); const INTVAL totalsize = attrs->rows * attrs->cols; PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize")); @@ -1413,7 +1768,7 @@ Get a PMCMatrix2D from the current matrix =back -=end +=cut */ } diff --git a/src/pmc/nummatrix2d.pmc b/src/pmc/nummatrix2d.pmc index b3a1649..11e59a2 100644 --- a/src/pmc/nummatrix2d.pmc +++ b/src/pmc/nummatrix2d.pmc @@ -1,9 +1,20 @@ #include "pla.h" - -#define ALLOCATE_STORAGE(s) (FLOATVAL *)mem_sys_allocate_zeroed(s * sizeof (FLOATVAL)) #define PLATYPENAME "NumMatrix2D" -#define DECLATTRS(p, a) Parrot_NumMatrix2D_attributes * const (a) = \ - (Parrot_NumMatrix2D_attributes *)((p)->data) + +/* + +=head1 NumMatrix2D + +=head2 Description + +NumMatrix2D is a 2-dimensional real-valued matrix type for the Parrot Virtual +Machine. It supports a variety of element-wise and matrix operations, including +bindings to the BLAS and LAPACK libraries. + +=cut + +*/ + INTVAL __PLA_Type_NumMatrix2D; @@ -21,7 +32,7 @@ INTVAL __PLA_Type_NumMatrix2D; static void resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) { - DECLATTRS(self, attrs); + DECLATTRS_NumMatrix2D(self, attrs); /* Store the old values */ const INTVAL old_rows = attrs->rows; const INTVAL old_cols = attrs->cols; @@ -32,7 +43,7 @@ resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) const INTVAL new_rows = INDEX_MAX(old_rows, row + 1); const INTVAL new_cols = INDEX_MAX(old_cols, col + 1); const INTVAL newsize = new_rows * new_cols; - FLOATVAL * new_s = ALLOCATE_STORAGE(newsize); + FLOATVAL * new_s = ALLOCATE_STORAGE_NumMatrix2D(newsize); INTVAL i, j; for (i = 0; i < old_rows; i++) { @@ -48,18 +59,37 @@ resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) mem_sys_free(old_s); } +static void +multiply_matrices(PARROT_INTERP, FLOATVAL alpha, PMC * A, PMC * B, FLOATVAL beta, PMC * C) +{ + DECLATTRS_NumMatrix2D(A, attrs_a); + DECLATTRS_NumMatrix2D(B, attrs_b); + DECLATTRS_NumMatrix2D(C, attrs_c); + if (attrs_c->rows != attrs_a->rows) + Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": A, C indices do not match in gemm"); + if (attrs_c->cols != attrs_b->cols) + Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": B, C indices do not match in gemm"); + if (attrs_b->rows != attrs_a->cols) + Parrot_ex_throw_from_c_args(interp, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": A, B indeces do not match in gemm"); + call_dgemm(alpha, attrs_a->flags, attrs_a->storage, attrs_a->rows, attrs_a->cols, + attrs_b->flags, attrs_b->storage, attrs_b->cols, beta, attrs_c->storage); +} + /* If the matrix is lazily transposed, actually transpose the physical memory layout. This is necessary for calculations, especially BLAS calculations, which aren't lazy-transpose-aware. */ static void normalize_lazy_transpose(PARROT_INTERP, PMC * self) { - DECLATTRS(self, attrs); + DECLATTRS_NumMatrix2D(self, attrs); if (IS_TRANSPOSED(attrs->flags)) { const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL size = rows_size * cols_size; - FLOATVAL * const new_s = ALLOCATE_STORAGE(size); + FLOATVAL * const new_s = ALLOCATE_STORAGE_NumMatrix2D(size); FLOATVAL * const old_s = attrs->storage; INTVAL i, j; @@ -100,7 +130,7 @@ static void init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, INTVAL cols_size, PMC * values) { - DECLATTRS(self, attrs); + DECLATTRS_NumMatrix2D(self, attrs); FLOATVAL * s; INTVAL self_rows, self_cols, i, j, num = 0; const INTVAL init_elems = VTABLE_elements(interp, values); @@ -125,7 +155,7 @@ init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, static void add_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) { - DECLATTRS(self, attrs); + DECLATTRS_NumMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; FLOATVAL * const s = attrs->storage; @@ -142,7 +172,7 @@ add_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) static void multiply_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) { - DECLATTRS(self, attrs); + DECLATTRS_NumMatrix2D(self, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; FLOATVAL * const s = attrs->storage; @@ -160,8 +190,8 @@ multiply_scalar_float(PARROT_INTERP, PMC * self, FLOATVAL v) static void add_matrices(PARROT_INTERP, PMC * A, PMC * B, INTVAL sub) { - DECLATTRS(A, attrs_a); - DECLATTRS(B, attrs_b); + DECLATTRS_NumMatrix2D(A, attrs_a); + DECLATTRS_NumMatrix2D(B, attrs_b); const INTVAL rows = attrs_a->rows; const INTVAL cols = attrs_a->cols; FLOATVAL * const s_a = attrs_a->storage; @@ -200,24 +230,28 @@ pmclass NumMatrix2D dynpmc auto_attrs provides matrix provides numericmatrix { /* -=head1 VTABLEs +=head2 VTABLEs + +=head3 System VTABLEs =over 4 -=item* init +=item * init Create a new NumMatrix2D -=item* destroy +=item * destroy Destrow the matrix and free it's allocated storage +=back + =cut */ VTABLE void init() { - DECLATTRS(SELF, a); + DECLATTRS_NumMatrix2D(SELF, a); a->storage = NULL; a->rows = 0; a->cols = 0; @@ -226,7 +260,7 @@ Destrow the matrix and free it's allocated storage } VTABLE void destroy() { - DECLATTRS(SELF, a); + DECLATTRS_NumMatrix2D(SELF, a); FLOATVAL * const s = a->storage; if (s) mem_sys_free(s); @@ -234,30 +268,42 @@ Destrow the matrix and free it's allocated storage /* -=item* get_number_keyed +=head3 Keyed Lookup VTABLEs + +In each of these cases, the specified Key PMC must have exactly two elements +to specify a location in the matrix. + +Attempting to retrieve a value outside the boundaries of the matrix will throw +an OUT_OF_BOUNDS exception. + +=over 4 + +=item * get_number_keyed Get the number at the location X, Y. The key must have two elements. -=item* get_integer_keyed +=item * get_integer_keyed Get the integer at the location X, Y. The key must have two elements. -=item* get_string_keyed +=item * get_string_keyed Get a string representation of the number at the location X, Y. The key must have two elements. -=item* get_pmc_keyed +=item * get_pmc_keyed Get a Float PMC of the number at the location X, Y. The key must have two elements. +=back + =cut */ VTABLE FLOATVAL get_number_keyed(PMC * key) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; INTVAL rows, cols; @@ -287,32 +333,44 @@ elements. /* -=item* get_number_keyed_int +=head3 Integer-Keyed Lookup VTABLES + +These VTABLEs treat the matrix, which is a contiguous region in memory, as a +linear array of values. The matrix data is stored by rows. + +These routines are used for low-level access. Attempting to access a value +outside the bounds of the matrix will throw an OUT_OF_BOUNDS exception. + +=over 4 + +=item * get_number_keyed_int Treating the matrix memory block like a C array, get the Nth number. The array is arranged in memory row by row. -=item* get_integer_keyed_int +=item * get_integer_keyed_int Treating the matrix memory block like a C array, get the Nth integer. The array is arranged in memory row by row. -=item* get_string_keyed_int +=item * get_string_keyed_int Treating the matrix memory block like a C array, get the string representation of the Nth number. The array is arranged in memory row by row. -=item* get_pmc_keyed_int +=item * get_pmc_keyed_int Treating the matrix memory block like a C array, get a Float PMC for the Nth number. The array is arranged in memory row by row. +=back + =cut */ VTABLE FLOATVAL get_number_keyed_int(INTVAL key) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); const INTVAL total_size = attrs->rows * attrs->cols; if (key >= total_size) { Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, @@ -340,28 +398,40 @@ Nth number. The array is arranged in memory row by row. /* -=item* set_number_keyed +=head3 Keyed Setter VTABLES + +These VTABLEs insert new values into the matrix at a point specified by the +Key PMC. The Key PMC must have exactly two elements. If the matrix is not large +enough to accomodate the specified location, it will be grown with zero-padding +so that it is at least large enough to hold the specified point and all existing +data. + +=over 4 + +=item * set_number_keyed Set the number at position (X, Y), growing the matrix if necessary and padding empty spaces with 0.0. The key must have two elements. -=item* set_integer_keyed +=item * set_integer_keyed Set the integer at position (X, Y), growing the matrix if necessary and padding empty spaces with 0.0. The key must have two elements. -=item* set_pmc_keyed +=item * set_pmc_keyed Get the numeric value from the PMC argument and set it at position (X, Y), growing the matrix if necessary and padding empty spaces with 0.0. The key must have two elements. +=back + =cut */ VTABLE void set_number_keyed(PMC * key, FLOATVAL value) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); INTVAL rows, cols, rows_size = attrs->rows, cols_size = attrs->cols; GET_KEY_INDICES_ROWMAJOR(INTERP, key, rows, cols); if (rows >= rows_size || cols >= cols_size) { @@ -383,8 +453,47 @@ must have two elements. VTABLE_set_number_keyed(INTERP, SELF, key, v); } +/* + +=head3 Integer-Keyed Setter VTABLEs + +These VTABLEs treat the matrix as a linear array in memory and allow fast +lookup based on the integer offset of values in the array. These are low-level +routines and are not intended for general use. + +Unlike the PMC-keyed VTABLEs, these routines will not automatically grow the +matrix if an index is provided which is outside the boundaries of the matrix. +In that case, an OUT_OF_BOUNDS exception will be thrown. + +=over 4 + +=item * set_pmc_keyed_int + +Set a PMC at the specified location. The rules for extracting a complex value +out of the input PMC are the same as used for set_pmc_keyed() + +=item * set_number_keyed_int + +Set the complex value at the specified location to the given real value. The +complex part of the value will be 0. + +=item * set_integer_keyed_int + +Set the complex value at the specified location to the given real value. The +complex part of the value will be 0 + +=item * set_string_keyed_int + +Convert the string to a Complex PMC, and set that at the specified location. + +=back + +=cut + +*/ + VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL value) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); const INTVAL total_size = attrs->rows * attrs->cols; if (key >= total_size) { Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, @@ -402,19 +511,39 @@ must have two elements. VTABLE_set_number_keyed_int(INTERP, SELF, key, f); } + /* TODO: set_string_keyed_int */ + /* -=item* get_string +=head3 Miscellaneous VTABLEs + +=over 4 + +=item * get_string + +Get a string representation of the matrix, suitable for printing to the console + +=item * get_attr_string + +Get a named attribute. The name can be one of "rows", "cols", or "size". + +=item * clone -Get a string representation of the matrix. +Clone the matrix + +=item * is_equal + +Determine if two matrices are equal in size and composition. + +=back =cut */ VTABLE STRING *get_string() { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); INTVAL i, j; PMC * const builder = Parrot_pmc_new(INTERP, enum_class_StringBuilder); STRING * const newline = Parrot_str_new(INTERP, "\n", 1); @@ -437,25 +566,157 @@ Get a string representation of the matrix. return VTABLE_get_string(INTERP, builder); } + VTABLE PMC * get_attr_str(STRING * idx) { + DECLATTRS_NumMatrix2D(SELF, attrs); + if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "rows"))) { + PMC * const rows = Parrot_pmc_new(INTERP, enum_class_Integer); + VTABLE_set_integer_native(INTERP, rows, attrs->rows); + return rows; + } + else if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "cols"))) { + PMC * const cols = Parrot_pmc_new(INTERP, enum_class_Integer); + VTABLE_set_integer_native(INTERP, cols, attrs->cols); + return cols; + } + else if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "size"))) { + PMC * const size = Parrot_pmc_new(INTERP, enum_class_Integer); + VTABLE_set_integer_native(INTERP, size, attrs->cols * attrs->rows); + return size; + } + return PMCNULL; + } + + VTABLE PMC * clone() { + PMC * const c = Parrot_pmc_new(INTERP, SELF->vtable->base_type); + DECLATTRS_NumMatrix2D(SELF, old_atts); + DECLATTRS_NumMatrix2D(c, new_atts); + INTVAL const newsize = old_atts->rows * old_atts->cols; + FLOATVAL * const old_s = old_atts->storage; + FLOATVAL * const new_s = ALLOCATE_STORAGE_NumMatrix2D(newsize); + memcpy(new_s, old_s, newsize * sizeof(FLOATVAL)); + memcpy(new_atts, old_atts, sizeof(Parrot_NumMatrix2D_attributes)); + new_atts->storage = new_s; + return c; + } + + VTABLE INTVAL is_equal(PMC * other) { + if (other->vtable->base_type == SELF->vtable->base_type) { + DECLATTRS_NumMatrix2D(SELF, self_attrs); + DECLATTRS_NumMatrix2D(other, other_attrs); + const INTVAL self_rows = self_attrs->rows; + const INTVAL self_cols = self_attrs->cols; + const INTVAL self_flags = self_attrs->flags; + const INTVAL other_rows = other_attrs->rows; + const INTVAL other_cols = other_attrs->cols; + const INTVAL other_flags = other_attrs->flags; + FLOATVAL * const self_s = self_attrs->storage; + FLOATVAL * const other_s = other_attrs->storage; + INTVAL i, j; + + if (self_rows != other_rows || self_cols != other_cols) + return 0; + + for (i = 0; i < self_rows; i++) { + for (j = 0; j < self_cols; j++) { + const FLOATVAL self_value = + ITEM_XY(self_s, self_flags, self_rows, self_cols, i, j); + const FLOATVAL other_value = + ITEM_XY(other_s, other_flags, other_rows, other_cols, i, j); + if (!floats_are_equal(self_value, other_value)) + return 0; + } + } + return 1; + } + return 0; + } + /* -=item* add(NumMatrix) +=head3 Serialization/Deserialization VTABLEs -Add two matrices together, element-by-element. +=over 4 -=item* add(DEFAULT) +=item * freeze + +Freeze the PMC for serialization to a string suitable for long-term storage in +a file. + +=item * thaw + +Thaw a serialized PMC + +=back + +=cut + +*/ + + VTABLE void freeze(PMC *info) { + DECLATTRS_NumMatrix2D(SELF, attrs); + INTVAL const rows = attrs->rows; + INTVAL const cols = attrs->cols; + INTVAL const flags = attrs->flags; + INTVAL i, j; + FLOATVAL * const s = attrs->storage; + VTABLE_push_integer(INTERP, info, rows); + VTABLE_push_integer(INTERP, info, cols); + VTABLE_push_integer(INTERP, info, flags); + for (i = 0; i < rows; i++) { + for (j = 0; j < cols; j++) { + const FLOATVAL f = ITEM_XY(s, flags, rows, cols, i, j); + VTABLE_push_float(INTERP, info, f); + } + } + } + + VTABLE void thaw(PMC *info) { + DECLATTRS_NumMatrix2D(SELF, attrs); + INTVAL const rows = VTABLE_shift_integer(INTERP, info); + INTVAL const cols = VTABLE_shift_integer(INTERP, info); + INTVAL const flags = VTABLE_shift_integer(INTERP, info); + INTVAL i, j; + FLOATVAL * s; + attrs->rows = 0; + attrs->cols = 0; + attrs->storage = NULL; + attrs->flags = 0; + resize_matrix(INTERP, SELF, rows - 1, cols - 1); + s = attrs->storage; + attrs->flags = flags; + for (i = 0; i < rows; i++) { + for (j = 0; j < cols; j++) { + const FLOATVAL f = VTABLE_shift_float(INTERP, info); + ITEM_XY(s, flags, rows, cols, i, j) = f; + } + } + } + +/* + +=head3 Addition VTABLEs + +=over 4 + +=item * add(NumMatrix2D) + +Add two matrices together, element-by-element + +=item * add(DEFAULT) Take the number representation of the argument PMC and add it to every element in the matrix. -=item* add_float +=item * add_float Add the float value to every element in the matrix. -=item* add_int +=item * add_int Add the integer value to every element in the matrix. +=back + =cut */ @@ -487,13 +748,19 @@ Add the integer value to every element in the matrix. /* -=item* i_add(NumMatrix2D) +=head3 In-Place Addition VTABLEs -=item* i_add(DEFAULT) +=over 4 + +=item * i_add(NumMatrix2D) + +=item * i_add(DEFAULT) -=item* i_add_int +=item * i_add_int -=item* i_add_float +=item * i_add_float + +=back =cut @@ -518,23 +785,29 @@ Add the integer value to every element in the matrix. /* -=item* subtract(NumMatrix) +=head3 Subtraction VTABLEs + +=over 4 + +=item * subtract(NumMatrix) Add two matrices together, element-by-element. -=item* subtract(DEFAULT) +=item * subtract(DEFAULT) Take the number representation of the argument PMC and add it to every element in the matrix. -=item* subtract_float +=item * subtract_float Add the float value to every element in the matrix. -=item* subtract_int +=item * subtract_int Add the integer value to every element in the matrix. +=back + =cut */ @@ -566,13 +839,19 @@ Add the integer value to every element in the matrix. /* -=item* i_subtract(NumMatrix2D) +=head3 In-Place Subtraction VTABLES + +=over 4 + +=item * i_subtract(NumMatrix2D) + +=item * i_subtract(DEFAULT) -=item* i_subtract(DEFAULT) +=item * i_subtract_int -=item* i_subtract_int +=item * i_subtract_float -=item* i_subtract_float +=back =cut @@ -597,17 +876,27 @@ Add the integer value to every element in the matrix. /* -=item* multiply(NumMatrix2D) +=head3 Multiplication VTABLEs + +=over 4 + +=item * multiply(NumMatrix2D) -=item* multiply(DEFAULT) +=item * multiply(DEFAULT) + +=item * multiply_int + +=item * multiply_float + +=back =cut */ MULTI PMC *multiply(NumMatrix2D *value, PMC *dest) { - DECLATTRS(SELF, selfattr); - DECLATTRS(value, valattr); + DECLATTRS_NumMatrix2D(SELF, selfattr); + DECLATTRS_NumMatrix2D(value, valattr); const INTVAL new_rows = selfattr->rows; const INTVAL new_cols = valattr->cols; @@ -618,7 +907,7 @@ Add the integer value to every element in the matrix. dest = Parrot_pmc_new(INTERP, VTABLE_type(INTERP, SELF)); resize_matrix(INTERP, dest, new_rows - 1, new_cols - 1); - call_dgemm(INTERP, 1.0, SELF, value, 0.0, dest); + multiply_matrices(INTERP, 1.0, SELF, value, 0.0, dest); return dest; } @@ -629,16 +918,6 @@ Add the integer value to every element in the matrix. return dest; } -/* - -=item* multiply_int - -=item* multiply_float - -=cut - -*/ - VTABLE PMC *multiply_int(INTVAL v, PMC * dest) { dest = VTABLE_clone(INTERP, SELF); multiply_scalar_float(INTERP, dest, (FLOATVAL)v); @@ -653,9 +932,19 @@ Add the integer value to every element in the matrix. /* -=item* i_multiply(NumMatrix2D) +=head3 In-Place Multiplication VTABLEs -=item* i_multiply(DEFAULT) +=over 4 + +=item * i_multiply(NumMatrix2D) + +=item * i_multiply(DEFAULT) + +=item * i_multiply_int + +=item * i_multiply_float + +=back =cut @@ -663,7 +952,7 @@ Add the integer value to every element in the matrix. MULTI void i_multiply(NumMatrix2D* value) { PMC * const temp = VTABLE_clone(INTERP, SELF); - call_dgemm(INTERP, 1.0, temp, value, 0.0, SELF); + multiply_matrices(INTERP, 1.0, temp, value, 0.0, SELF); } MULTI void i_multiply(DEFAULT* value) { @@ -671,16 +960,6 @@ Add the integer value to every element in the matrix. multiply_scalar_float(INTERP, SELF, v); } -/* - -=item i_multiply_int - -=item i_multiply_float - -=cut - -*/ - VTABLE void i_multiply_int(INTVAL v) { multiply_scalar_float(INTERP, SELF, (FLOATVAL)v); } @@ -689,154 +968,19 @@ Add the integer value to every element in the matrix. multiply_scalar_float(INTERP, SELF, v); } -/* - -=item* get_attr_str - -=cut - -*/ - - VTABLE PMC * get_attr_str(STRING * idx) { - DECLATTRS(SELF, attrs); - if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "rows"))) { - PMC * const rows = Parrot_pmc_new(INTERP, enum_class_Integer); - VTABLE_set_integer_native(INTERP, rows, attrs->rows); - return rows; - } - else if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "cols"))) { - PMC * const cols = Parrot_pmc_new(INTERP, enum_class_Integer); - VTABLE_set_integer_native(INTERP, cols, attrs->cols); - return cols; - } - else if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "size"))) { - PMC * const size = Parrot_pmc_new(INTERP, enum_class_Integer); - VTABLE_set_integer_native(INTERP, size, attrs->cols * attrs->rows); - return size; - } - return PMCNULL; - } - -/* - -=item* clone - -=cut - -*/ - - VTABLE PMC * clone() { - PMC * const c = Parrot_pmc_new(INTERP, SELF->vtable->base_type); - DECLATTRS(SELF, old_atts); - DECLATTRS(c, new_atts); - INTVAL const newsize = old_atts->rows * old_atts->cols; - FLOATVAL * const old_s = old_atts->storage; - FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize); - memcpy(new_s, old_s, newsize * sizeof(FLOATVAL)); - memcpy(new_atts, old_atts, sizeof(Parrot_NumMatrix2D_attributes)); - new_atts->storage = new_s; - return c; - } /* -=item* is_equal +=head2 METHODS -=cut - -*/ - - VTABLE INTVAL is_equal(PMC * other) { - if (other->vtable->base_type == SELF->vtable->base_type) { - DECLATTRS(SELF, self_attrs); - DECLATTRS(other, other_attrs); - const INTVAL self_rows = self_attrs->rows; - const INTVAL self_cols = self_attrs->cols; - const INTVAL self_flags = self_attrs->flags; - const INTVAL other_rows = other_attrs->rows; - const INTVAL other_cols = other_attrs->cols; - const INTVAL other_flags = other_attrs->flags; - FLOATVAL * const self_s = self_attrs->storage; - FLOATVAL * const other_s = other_attrs->storage; - INTVAL i, j; - - if (self_rows != other_rows || self_cols != other_cols) - return 0; - - for (i = 0; i < self_rows; i++) { - for (j = 0; j < self_cols; j++) { - const FLOATVAL self_value = - ITEM_XY(self_s, self_flags, self_rows, self_cols, i, j); - const FLOATVAL other_value = - ITEM_XY(other_s, other_flags, other_rows, other_cols, i, j); - if (!floats_are_equal(self_value, other_value)) - return 0; - } - } - return 1; - } - return 0; - } - -/* - -=item* freeze - -=item* thaw - -=cut - -*/ - - VTABLE void freeze(PMC *info) { - DECLATTRS(SELF, attrs); - INTVAL const rows = attrs->rows; - INTVAL const cols = attrs->cols; - INTVAL const flags = attrs->flags; - INTVAL i, j; - FLOATVAL * const s = attrs->storage; - VTABLE_push_integer(INTERP, info, rows); - VTABLE_push_integer(INTERP, info, cols); - VTABLE_push_integer(INTERP, info, flags); - for (i = 0; i < rows; i++) { - for (j = 0; j < cols; j++) { - const FLOATVAL f = ITEM_XY(s, flags, rows, cols, i, j); - VTABLE_push_float(INTERP, info, f); - } - } - } - - VTABLE void thaw(PMC *info) { - DECLATTRS(SELF, attrs); - INTVAL const rows = VTABLE_shift_integer(INTERP, info); - INTVAL const cols = VTABLE_shift_integer(INTERP, info); - INTVAL const flags = VTABLE_shift_integer(INTERP, info); - INTVAL i, j; - FLOATVAL * s; - attrs->rows = 0; - attrs->cols = 0; - attrs->storage = NULL; - attrs->flags = 0; - resize_matrix(INTERP, SELF, rows - 1, cols - 1); - s = attrs->storage; - attrs->flags = flags; - for (i = 0; i < rows; i++) { - for (j = 0; j < cols; j++) { - const FLOATVAL f = VTABLE_shift_float(INTERP, info); - ITEM_XY(s, flags, rows, cols, i, j) = f; - } - } - } - -/* - -=back +=over 4 -=head1 METHODS +=item * resize() -=over 4 +Resize the matrix to include at least the specified number of rows and columns. -=item resize() +Resizing the matrix never causes the matrix to shrink. If you need a subset of +the matrix, use get_block instead. =cut @@ -848,10 +992,11 @@ Add the integer value to every element in the matrix. /* -=item fill() +=item * fill() Fill the matrix with a single value. if sizes are provided, fill to those -sizes, growing the matrix if needed. +sizes, growing the matrix if needed. Elements outside the specified area are +unaffected. Calling fill() never causes the matrix to shrink. =cut @@ -861,7 +1006,7 @@ sizes, growing the matrix if needed. INTVAL rows_size :optional, INTVAL has_rows_size :opt_flag, INTVAL cols_size :optional, INTVAL has_cols_size :opt_flag ) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); FLOATVAL * s = attrs->storage; INTVAL const curr_rows_size = attrs->rows; INTVAL const curr_cols_size = attrs->cols; @@ -887,7 +1032,9 @@ sizes, growing the matrix if needed. /* -=item item_at() +=item * item_at() + +Return a single Float PMC from the item at the specified coordinates =cut @@ -896,7 +1043,7 @@ sizes, growing the matrix if needed. METHOD item_at(INTVAL row, INTVAL col, FLOATVAL value :optional, INTVAL has_value :opt_flag) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; if (row >= rows || col >= cols || row < 0 || col < 0) { @@ -912,16 +1059,19 @@ sizes, growing the matrix if needed. /* -=item transpose() +=item * transpose() -Transposes the matrix lazily. +Transposes the matrix lazily. This operation is O(1). Some operations, such as +mathematical operations do not work on a matrix which has been lazily +transposed, so those operations will force the matrix memory to be eagerly +transposed. =cut */ METHOD transpose() { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); INTVAL transposed = IS_TRANSPOSED(attrs->flags); SWAP_XY(attrs); @@ -936,21 +1086,22 @@ Transposes the matrix lazily. /* -=item mem_transpose() +=item * mem_transpose() -Transposes the actual data storage of the matrix. More expensive up-front -than the transpose() method. +Transposes the actual data storage of the matrix. More expensive O(n) up-front +than the transpose() method, but the resulting memory structure is more suitable +for use in certain mathematical operations. =cut */ METHOD mem_transpose() { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; const INTVAL newsize = rows * cols; - FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize); + FLOATVAL * const new_s = ALLOCATE_STORAGE_NumMatrix2D(newsize); FLOATVAL * const old_s = attrs->storage; INTVAL i, j; @@ -968,12 +1119,12 @@ than the transpose() method. /* -=item iterate_function_inplace() +=item * iterate_function_inplace() Calls a function for every element in the matrix, replacing the current value with the return value of the called function. -=item iterate_function_external() +=item * iterate_function_external() Calls a function for every element in the matrix, adding the result of each computation to a new matrix. Return the new matrix of results. @@ -983,13 +1134,13 @@ computation to a new matrix. Return the new matrix of results. */ METHOD iterate_function_inplace(PMC * func, PMC * args :slurpy) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; const INTVAL newsize = rows * cols; const INTVAL flags = attrs->flags; FLOATVAL * const old_s = attrs->storage; - FLOATVAL * const new_s = ALLOCATE_STORAGE(newsize); + FLOATVAL * const new_s = ALLOCATE_STORAGE_NumMatrix2D(newsize); INTVAL i, j; if (newsize == 0 || old_s == NULL) @@ -1009,7 +1160,7 @@ computation to a new matrix. Return the new matrix of results. } METHOD iterate_function_external(PMC * func, PMC * args :slurpy) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type); Parrot_NumMatrix2D_attributes * new_attrs; const INTVAL rows = attrs->rows; @@ -1041,12 +1192,12 @@ computation to a new matrix. Return the new matrix of results. /* -=item initialize_from_array() +=item * initialize_from_array() Initialize matrix values from a linear array, filling each row with data in order. -=item initialize_from_args() +=item * initialize_from_args() Initialize matrix values from an array of function arguments, filling each row with data in order. @@ -1067,9 +1218,14 @@ with data in order. /* -=item get_block +=item * get_block + +Get a specified sub-block of the matrix. If the bounds of the sub-block are +outside the bounds of the matrix, an OUT_OF_BOUNDS exception is thrown. -=item set_block +=item * set_block + +Set a block in the matrix, growing it if needed. =cut @@ -1077,7 +1233,7 @@ with data in order. METHOD get_block(INTVAL rows_idx, INTVAL cols_idx, INTVAL rows_size, INTVAL cols_size) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1094,7 +1250,7 @@ with data in order. PLATYPENAME ": Can not get block with negative size"); else { PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type); - DECLATTRS(new_matrix, new_attrs); + DECLATTRS_NumMatrix2D(new_matrix, new_attrs); FLOATVAL * new_s; resize_matrix(INTERP, new_matrix, rows_size - 1, cols_size - 1); new_s = new_attrs->storage; @@ -1109,8 +1265,8 @@ with data in order. } METHOD set_block(INTVAL rows_idx, INTVAL cols_idx, PMC * blck) { - DECLATTRS(SELF, self_attrs); - DECLATTRS(blck, blck_attrs); + DECLATTRS_NumMatrix2D(SELF, self_attrs); + DECLATTRS_NumMatrix2D(blck, blck_attrs); FLOATVAL * self_s = self_attrs->storage; FLOATVAL * const blck_s = blck_attrs->storage; INTVAL self_rows = self_attrs->rows; @@ -1145,12 +1301,16 @@ with data in order. /* -=item* gemm +=item * gemm Calculates the matrix equation: Z = aAB + bC +The matrices must all be NumMatrix2D, or must be convertable to it. The +matrix SELF is not used in the calculation, but the result matrix will have the +same type as SELF. + =cut */ @@ -1159,30 +1319,34 @@ Calculates the matrix equation: A = convert_to_NumMatrix2D(interp, A, 0); B = convert_to_NumMatrix2D(interp, B, 0); C = convert_to_NumMatrix2D(interp, C, 1); - call_dgemm(INTERP, alpha, A, B, beta, C); + multiply_matrices(INTERP, alpha, A, B, beta, C); RETURN(PMC* C); } /* -=item row_combine(srcidx, destidx, gain) +=item * row_combine(srcidx, destidx, gain) -add a multiple of the source row to the destination row. +add a multiple of the source row to the destination row. If either of the row +indices are outside the bounds of the matrix, an OUT_OF_BOUNDS exception is +thrown. -=item row_scale(idx, gain) +=item * row_scale(idx, gain) -Multiply all elements in the row by a gain factor. +Multiply all elements in the row by a gain factor. If the row index is outside +the bounds of the matrix and OUT_OF_BOUNDS exception is thrown. -=item row_swap(idx_a, idx_b) +=item * row_swap(idx_a, idx_b) -Swap two rows +Swap two rows. If either of the row indices are outside the bounds of the +matrix, an OUT_OF_BOUNDS exception is thrown. =cut */ METHOD row_combine(INTVAL srcidx, INTVAL destidx, FLOATVAL gain) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1198,7 +1362,7 @@ Swap two rows } METHOD row_scale(INTVAL idx, FLOATVAL gain) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1214,7 +1378,7 @@ Swap two rows METHOD row_swap(INTVAL idx_a, INTVAL idx_b) { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); FLOATVAL * const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -1233,17 +1397,20 @@ Swap two rows /* -=item convert_to_number_matrix +=item * convert_to_number_matrix -Get a NumMatrix2D from the current matrix +Get a NumMatrix2D from the current matrix. If the matrix is already a +NumMatrix2D, return a clone. -=item convert_to_complex_matrix +=item * convert_to_complex_matrix -Get a ComplexMatrix2D from the current matrix +Get a ComplexMatrix2D from the current matrix. If the matrix is already a +ComplexMatrix2D, return a clone. -=item convert_to_pmc_matrix +=item * convert_to_pmc_matrix -Get a PMCMatrix2D from the current matrix +Get a PMCMatrix2D from the current matrix. If the matrix is already a +PMCMatrix2D, return a clone. =cut @@ -1256,7 +1423,7 @@ Get a PMCMatrix2D from the current matrix } METHOD convert_to_complex_matrix() { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_ComplexMatrix2D); const INTVAL totalsize = attrs->rows * attrs->cols; PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize")); @@ -1271,7 +1438,7 @@ Get a PMCMatrix2D from the current matrix } METHOD convert_to_pmc_matrix() { - DECLATTRS(SELF, attrs); + DECLATTRS_NumMatrix2D(SELF, attrs); PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_PMCMatrix2D); const INTVAL totalsize = attrs->rows * attrs->cols; PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize")); @@ -1285,13 +1452,11 @@ Get a PMCMatrix2D from the current matrix RETURN(PMC * d); } - - /* =back -=end +=cut */ } diff --git a/src/pmc/pmcmatrix2d.pmc b/src/pmc/pmcmatrix2d.pmc index 6fc44c2..9a7e957 100644 --- a/src/pmc/pmcmatrix2d.pmc +++ b/src/pmc/pmcmatrix2d.pmc @@ -1,9 +1,18 @@ #include "pla.h" - -#define ALLOCATE_STORAGE(s) (PMC **)mem_sys_allocate_zeroed(s * sizeof (PMC *)) #define PLATYPENAME "PMCMatrix2D" -#define DECLATTRS(p, a) Parrot_PMCMatrix2D_attributes * const (a) = \ - (Parrot_PMCMatrix2D_attributes *)((p)->data) + +/* + +=head1 PMCMatrix2D + +=head2 Description + +PMCMatrix2D is a 2-dimensional matrix-like container object for the Parrot +Virtual Machine. + +=cut + +*/ INTVAL __PLA_Type_PMCMatrix2D; @@ -21,7 +30,7 @@ INTVAL __PLA_Type_PMCMatrix2D; static void resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) { - DECLATTRS(self, attrs); + DECLATTRS_PMCMatrix2D(self, attrs); /* Store the old values */ const INTVAL old_rows = attrs->rows; const INTVAL old_cols = attrs->cols; @@ -32,7 +41,7 @@ resize_matrix(PARROT_INTERP, PMC * self, INTVAL row, INTVAL col) const INTVAL new_rows = INDEX_MAX(old_rows, row + 1); const INTVAL new_cols = INDEX_MAX(old_cols, col + 1); const INTVAL newsize = new_rows * new_cols; - PMC ** new_s = ALLOCATE_STORAGE(newsize); + PMC ** new_s = ALLOCATE_STORAGE_PMCMatrix2D(newsize); INTVAL i, j; for (i = 0; i < old_rows; i++) { @@ -59,7 +68,7 @@ static void init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, INTVAL cols_size, PMC * values) { - DECLATTRS(self, attrs); + DECLATTRS_PMCMatrix2D(self, attrs); PMC ** s; INTVAL self_rows, self_cols, i, j, num = 0; const INTVAL init_elems = VTABLE_elements(interp, values); @@ -86,12 +95,12 @@ init_from_pmc_array(PARROT_INTERP, PMC * self, INTVAL rows_size, static void normalize_lazy_transpose(PARROT_INTERP, PMC * self) { - DECLATTRS(self, attrs); + DECLATTRS_PMCMatrix2D(self, attrs); if (IS_TRANSPOSED(attrs->flags)) { const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL size = rows_size * cols_size; - PMC ** const new_s = ALLOCATE_STORAGE(size); + PMC ** const new_s = ALLOCATE_STORAGE_PMCMatrix2D(size); PMC ** const old_s = attrs->storage; INTVAL i, j; @@ -118,8 +127,34 @@ pmclass PMCMatrix2D dynpmc auto_attrs provides matrix { __PLA_Type_PMCMatrix2D = entry; } +/* + +=head1 VTABLEs + +=head2 System VTABLEs + +=over 4 + +=item * init + +Initialize the new PMC + +=item * destroy + +Destroy the PMC and free all associated memory + +=item * mark + +Mark the contents of the matrix for GC + +=back + +=cut + +*/ + VTABLE void init() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); attrs->storage = NULL; attrs->rows = 0; attrs->cols = 0; @@ -129,14 +164,14 @@ pmclass PMCMatrix2D dynpmc auto_attrs provides matrix { } VTABLE void destroy() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC ** const s = attrs->storage; if (s) mem_sys_free(s); } VTABLE void mark() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC ** s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -149,10 +184,42 @@ pmclass PMCMatrix2D dynpmc auto_attrs provides matrix { } } -/* Get VTABLEs */ +/* + +=head2 Keyed Lookup VTABLEs + +In each of these cases, the specified Key PMC must have exactly two elements +to specify a location in the matrix. + +Attempting to retrieve a value outside the boundaries of the matrix will throw +an OUT_OF_BOUNDS exception. + +=over 4 + +=item * get_number_keyed + +Get the PMC at the location X, Y. + +=item * get_integer_keyed + +Get the integer from the PMC at the location X, Y. + +=item * get_string_keyed + +Get a string representation of the PMC at the location X, Y. + +=item * get_pmc_keyed + +Get a PMC at the location X, Y. + +=back + +=cut + +*/ VTABLE PMC * get_pmc_keyed(PMC * key) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); INTVAL cols, rows, cols_size, rows_size; cols_size = attrs->cols; rows_size = attrs->rows; @@ -163,8 +230,57 @@ pmclass PMCMatrix2D dynpmc auto_attrs provides matrix { return ITEM_XY(attrs->storage, attrs->flags, rows_size, cols_size, rows, cols); } + VTABLE INTVAL get_integer_keyed(PMC * key) { + PMC * const item = VTABLE_get_pmc_keyed(INTERP, SELF, key); + return VTABLE_get_integer(INTERP, item); + } + + VTABLE FLOATVAL get_number_keyed(PMC * key) { + PMC * const item = VTABLE_get_pmc_keyed(INTERP, SELF, key); + return VTABLE_get_number(INTERP, item); + } + + VTABLE STRING * get_string_keyed(PMC * key) { + PMC * const item = VTABLE_get_pmc_keyed(INTERP, SELF, key); + return VTABLE_get_string(INTERP, item); + } + +/* + +=head2 Integer-Keyed Lookup VTABLES + +These VTABLEs treat the matrix, which is a contiguous region in memory, as a +linear array of values. The matrix data is stored by rows. + +These routines are used for low-level access. Attempting to access a value +outside the bounds of the matrix will throw an OUT_OF_BOUNDS exception. + +=over 4 + +=item * get_number_keyed_int + +Get a floating point number from the PMC at the specified location. + +=item * get_integer_keyed_int + +Get an integer from the PMC at the specifed location. + +=item * get_string_keyed_int + +Get the string representation of the PMC at the specified location. + +=item * get_pmc_keyed_int + +Get the PMC at the specified location + +=back + +=cut + +*/ + VTABLE PMC * get_pmc_keyed_int(INTVAL key) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); const INTVAL total_size = attrs->cols * attrs->rows; if (key >= total_size || key < 0) { Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, @@ -173,31 +289,16 @@ pmclass PMCMatrix2D dynpmc auto_attrs provides matrix { return attrs->storage[key]; } - VTABLE INTVAL get_integer_keyed(PMC * key) { - PMC * const item = VTABLE_get_pmc_keyed(INTERP, SELF, key); - return VTABLE_get_integer(INTERP, item); - } - VTABLE INTVAL get_integer_keyed_int(INTVAL key) { PMC * const item = VTABLE_get_pmc_keyed_int(INTERP, SELF, key); return VTABLE_get_integer(INTERP, item); } - VTABLE FLOATVAL get_number_keyed(PMC * key) { - PMC * const item = VTABLE_get_pmc_keyed(INTERP, SELF, key); - return VTABLE_get_number(INTERP, item); - } - VTABLE FLOATVAL get_number_keyed_int(INTVAL key) { PMC * const item = VTABLE_get_pmc_keyed_int(INTERP, SELF, key); return VTABLE_get_number(INTERP, item); } - VTABLE STRING * get_string_keyed(PMC * key) { - PMC * const item = VTABLE_get_pmc_keyed(INTERP, SELF, key); - return VTABLE_get_string(INTERP, item); - } - VTABLE STRING * get_string_keyed_int(INTVAL key) { PMC * const item = VTABLE_get_pmc_keyed_int(INTERP, SELF, key); return VTABLE_get_string(INTERP, item); @@ -205,10 +306,46 @@ pmclass PMCMatrix2D dynpmc auto_attrs provides matrix { /* Set VTABLEs TODO: Update all these to follow HLL mappings +*/ +/* + +=head2 Keyed Setter VTABLES + +These VTABLEs insert new values into the matrix at a point specified by the +Key PMC. The Key PMC must have exactly two elements. If the matrix is not large +enough to accomodate the specified location, it will be grown with zero-padding +so that it is at least large enough to hold the specified point and all existing +data. + +=over 4 + +=item * set_number_keyed + +Create a Float PMC with the specified value, and insert it at the specified +location + +=item * set_integer_keyed + +Create an Integer PMC with the specified value, and insert it at the specified +location + +=item * set_pmc_keyed + +Set the PMC at the specified location. + +=item * set_string_keyed + +Create a String PMC with the specified value, and insert it at the specified +location + +=back + +=cut + */ VTABLE void set_pmc_keyed(PMC * key, PMC * value) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); INTVAL cols, rows, cols_size, rows_size; cols_size = attrs->cols; rows_size = attrs->rows; @@ -221,54 +358,115 @@ TODO: Update all these to follow HLL mappings ITEM_XY(attrs->storage, attrs->flags, rows_size, cols_size, rows, cols) = value; } - VTABLE void set_pmc_keyed_int(INTVAL key, PMC * value) { - DECLATTRS(SELF, attrs); - const INTVAL total_size = attrs->cols * attrs->rows; - if (key >= total_size || key < 0) { - Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, - PLATYPENAME ": indices out of bounds."); - } - attrs->storage[key] = value; - } - VTABLE void set_integer_keyed(PMC * key, INTVAL value) { PMC * const item = Parrot_pmc_new(INTERP, enum_class_Integer); VTABLE_set_integer_native(INTERP, item, value); VTABLE_set_pmc_keyed(INTERP, SELF, key, item); } - VTABLE void set_integer_keyed_int(INTVAL key, INTVAL value) { - PMC * const item = Parrot_pmc_new(INTERP, enum_class_Integer); - VTABLE_set_integer_native(INTERP, item, value); - VTABLE_set_pmc_keyed_int(INTERP, SELF, key, item); - } - VTABLE void set_number_keyed(PMC * key, FLOATVAL value) { PMC * const item = Parrot_pmc_new(INTERP, enum_class_Float); VTABLE_set_number_native(INTERP, item, value); VTABLE_set_pmc_keyed(INTERP, SELF, key, item); } - VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL value) { - PMC * const item = Parrot_pmc_new(INTERP, enum_class_Float); - VTABLE_set_number_native(INTERP, item, value); - VTABLE_set_pmc_keyed_int(INTERP, SELF, key, item); - } - VTABLE void set_string_keyed(PMC * key, STRING * value) { PMC * const item = Parrot_pmc_new(INTERP, enum_class_String); VTABLE_set_string_native(INTERP, item, value); VTABLE_set_pmc_keyed(INTERP, SELF, key, item); } +/* + +=head2 Integer-Keyed Setter VTABLEs + +These VTABLEs treat the matrix as a linear array in memory and allow fast +lookup based on the integer offset of values in the array. These are low-level +routines and are not intended for general use. + +Unlike the PMC-keyed VTABLEs, these routines will not automatically grow the +matrix if an index is provided which is outside the boundaries of the matrix. +In that case, an OUT_OF_BOUNDS exception will be thrown. + +=over 4 + +=item * set_pmc_keyed_int + +Set a PMC at the specified location. + +=item * set_number_keyed_int + +Create a Float PMC with the specified value, and insert it at the specified +location + +=item * set_integer_keyed_int + +Create an Integer PMC with the specified value, and insert it at the specified +location + +=item * set_string_keyed_int + +Create a String PMC with the specified value, and insert it at the specified +location + +=back + +=cut + +*/ + + VTABLE void set_pmc_keyed_int(INTVAL key, PMC * value) { + DECLATTRS_PMCMatrix2D(SELF, attrs); + const INTVAL total_size = attrs->cols * attrs->rows; + if (key >= total_size || key < 0) { + Parrot_ex_throw_from_c_args(INTERP, NULL, EXCEPTION_OUT_OF_BOUNDS, + PLATYPENAME ": indices out of bounds."); + } + attrs->storage[key] = value; + } + VTABLE void set_string_keyed_int(INTVAL key, STRING * value) { PMC * const item = Parrot_pmc_new(INTERP, enum_class_String); VTABLE_set_string_native(INTERP, item, value); VTABLE_set_pmc_keyed_int(INTERP, SELF, key, item); } + VTABLE void set_number_keyed_int(INTVAL key, FLOATVAL value) { + PMC * const item = Parrot_pmc_new(INTERP, enum_class_Float); + VTABLE_set_number_native(INTERP, item, value); + VTABLE_set_pmc_keyed_int(INTERP, SELF, key, item); + } + +/* + +=head2 Miscellaneous VTABLEs + +=over 4 + +=item * get_string + +Get a string representation of the matrix, suitable for printing to the console + +=item * get_attr_string + +Get a named attribute. The name can be one of "rows", "cols", or "size". + +=item * clone + +Clone the matrix + +=item * is_equal + +Determine if two matrices are equal in size and composition. + +=back + +=cut + +*/ + VTABLE STRING *get_string() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); INTVAL rows, cols; PMC * builder = Parrot_pmc_new(INTERP, enum_class_StringBuilder); STRING * const strend = Parrot_str_new(INTERP, "}", 1); @@ -300,7 +498,7 @@ TODO: Update all these to follow HLL mappings /* TODO: Update this to account for transpositions */ VTABLE PMC * get_attr_str(STRING * idx) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); if (Parrot_str_equal(INTERP, idx, CONST_STRING(INTERP, "rows"))) { PMC * const rows = Parrot_pmc_new(INTERP, enum_class_Integer); VTABLE_set_integer_native(INTERP, rows, attrs->rows); @@ -320,8 +518,8 @@ TODO: Update all these to follow HLL mappings } VTABLE INTVAL is_equal(PMC *other) { - DECLATTRS(SELF, attrs); - DECLATTRS(other, oattr); + DECLATTRS_PMCMatrix2D(SELF, attrs); + DECLATTRS_PMCMatrix2D(other, oattr); PMC ** const s = attrs->storage; PMC ** const o = oattr->storage; const INTVAL sflags = attrs->flags; @@ -354,8 +552,8 @@ TODO: Update all these to follow HLL mappings VTABLE PMC* clone() { PMC * const other = Parrot_pmc_new(INTERP, SELF->vtable->base_type); - DECLATTRS(SELF, attrs); - DECLATTRS(other, oattr); + DECLATTRS_PMCMatrix2D(SELF, attrs); + DECLATTRS_PMCMatrix2D(other, oattr); PMC ** s, ** o; INTVAL self_rows, self_cols, i, j, o_rows, o_cols; @@ -381,16 +579,27 @@ TODO: Update all these to follow HLL mappings /* -=item* freeze +=head2 Serialization/Deserialization VTABLEs + +=over 4 + +=item * freeze + +Freeze the PMC for serialization to a string suitable for long-term storage in +a file. -=item* thaw +=item * thaw + +Thaw a serialized PMC + +=back =cut */ VTABLE void freeze(PMC *info) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); INTVAL const rows = attrs->rows; INTVAL const cols = attrs->cols; INTVAL const flags = attrs->flags; @@ -408,7 +617,7 @@ TODO: Update all these to follow HLL mappings } VTABLE void thaw(PMC *info) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); INTVAL const rows = VTABLE_shift_integer(INTERP, info); INTVAL const cols = VTABLE_shift_integer(INTERP, info); INTVAL const flags = VTABLE_shift_integer(INTERP, info); @@ -429,13 +638,22 @@ TODO: Update all these to follow HLL mappings } } - METHOD initialize_from_array(INTVAL rows_size, INTVAL cols_size, PMC *values) { - init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values); - } +/* - METHOD initialize_from_args(INTVAL rows_size, INTVAL cols_size, PMC *values :slurpy) { - init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values); - } +=head1 METHODS + +=over 4 + +=item * resize() + +Resize the matrix to include at least the specified number of rows and columns. + +Resizing the matrix never causes the matrix to shrink. If you need a subset of +the matrix, use get_block instead. + +=cut + +*/ METHOD resize(INTVAL rows, INTVAL cols) { resize_matrix(INTERP, SELF, rows - 1, cols - 1); @@ -443,10 +661,11 @@ TODO: Update all these to follow HLL mappings /* -=item fill() +=item * fill() Fill the matrix with a single value. if sizes are provided, fill to those -sizes, growing the matrix if needed. +sizes, growing the matrix if needed. Elements outside the specified area are +unaffected. Calling fill() never causes the matrix to shrink. =cut @@ -456,7 +675,7 @@ sizes, growing the matrix if needed. INTVAL x_size :optional, INTVAL has_rows_size :opt_flag, INTVAL y_size :optional, INTVAL has_cols_size :opt_flag ) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC ** s = attrs->storage; INTVAL const curr_rows_size = attrs->rows; INTVAL const curr_cols_size = attrs->cols; @@ -487,7 +706,9 @@ sizes, growing the matrix if needed. /* -=item item_at() +=item * item_at() + +Return a single PMC at the specified coordinates =cut @@ -495,7 +716,7 @@ sizes, growing the matrix if needed. METHOD item_at(INTVAL row, INTVAL col, PMC * value :optional, INTVAL has_value :opt_flag) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; const INTVAL flags = attrs->flags; @@ -513,16 +734,16 @@ sizes, growing the matrix if needed. /* -=item transpose() +=item * transpose() -Transposes the matrix. +Transposes the matrix lazily. This operation is O(1). =cut */ METHOD transpose() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); INTVAL tmp = 0; INTVAL transposed = IS_TRANSPOSED(attrs->flags); @@ -538,9 +759,9 @@ Transposes the matrix. /* -=item mem_transpose() +=item * mem_transpose() -Transposes the actual data storage of the matrix. More expensive up-front +Transposes the actual data storage of the matrix. More expensive O(n) up-front than the transpose() method. =cut @@ -548,11 +769,11 @@ than the transpose() method. */ METHOD mem_transpose() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL newsize = rows_size * cols_size; - PMC ** new_s = ALLOCATE_STORAGE(newsize); + PMC ** new_s = ALLOCATE_STORAGE_PMCMatrix2D(newsize); PMC ** old_s = attrs->storage; INTVAL i, j; @@ -570,22 +791,27 @@ than the transpose() method. /* -=item iterate_function_inplace() +=item * iterate_function_inplace() -Calls a function for every element in the array, replacing the current +Calls a function for every element in the matrix, replacing the current value with the return value of the called function. +=item * iterate_function_external() + +Calls a function for every element in the matrix, adding the result of each +computation to a new matrix. Return the new matrix of results. + =cut */ METHOD iterate_function_inplace(PMC * func, PMC * args :slurpy) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL newsize = rows_size * cols_size; PMC ** old_s = attrs->storage; - PMC ** new_s = ALLOCATE_STORAGE(newsize); + PMC ** new_s = ALLOCATE_STORAGE_PMCMatrix2D(newsize); INTVAL i, j; if (newsize == 0 || old_s == NULL) @@ -605,9 +831,9 @@ value with the return value of the called function. } METHOD iterate_function_external(PMC * func, PMC * args :slurpy) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC * const new_matrix = Parrot_pmc_new(INTERP, SELF->vtable->base_type); - DECLATTRS(new_matrix, new_attrs); + DECLATTRS_PMCMatrix2D(new_matrix, new_attrs); const INTVAL rows_size = attrs->rows; const INTVAL cols_size = attrs->cols; const INTVAL newsize = rows_size * cols_size; @@ -634,9 +860,48 @@ value with the return value of the called function. RETURN(PMC * new_matrix); } +/* + +=item * initialize_from_array() + +Initialize matrix values from a linear array, filling each row with data +in order. + +=item * initialize_from_args() + +Initialize matrix values from an array of function arguments, filling each row +with data in order. + +=cut + +*/ + + METHOD initialize_from_array(INTVAL rows_size, INTVAL cols_size, PMC *values) { + init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values); + } + + METHOD initialize_from_args(INTVAL rows_size, INTVAL cols_size, PMC *values :slurpy) { + init_from_pmc_array(INTERP, SELF, rows_size, cols_size, values); + } + +/* + +=item * get_block + +Get a specified sub-block of the matrix. If the bounds of the sub-block are +outside the bounds of the matrix, an OUT_OF_BOUNDS exception is thrown. + +=item * set_block + +Set a block in the matrix, growing it if needed. + +=cut + +*/ + METHOD get_block(INTVAL rows_idx, INTVAL cols_idx, INTVAL rows_size, INTVAL cols_size) { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC ** const s = attrs->storage; const INTVAL rows = attrs->rows; const INTVAL cols = attrs->cols; @@ -668,8 +933,8 @@ value with the return value of the called function. } METHOD set_block(INTVAL rows_idx, INTVAL cols_idx, PMC * blck) { - DECLATTRS(SELF, self_attrs); - DECLATTRS(blck, blck_attrs); + DECLATTRS_PMCMatrix2D(SELF, self_attrs); + DECLATTRS_PMCMatrix2D(blck, blck_attrs); PMC ** self_s = self_attrs->storage; PMC ** const blck_s = blck_attrs->storage; INTVAL self_rows = self_attrs->rows; @@ -704,24 +969,28 @@ value with the return value of the called function. /* -=item convert_to_number_matrix +=item * convert_to_number_matrix -Get a NumMatrix2D from the current matrix +Get a NumMatrix2D from the current matrix. If the matrix is already a +NumMatrix2D, return a clone. -=item convert_to_complex_matrix +=item * convert_to_complex_matrix -Get a ComplexMatrix2D from the current matrix +Get a ComplexMatrix2D from the current matrix. If the matrix is already a +ComplexMatrix2D, return a clone. -=item convert_to_pmc_matrix +=item * convert_to_pmc_matrix -Get a PMCMatrix2D from the current matrix +Get a PMCMatrix2D from the current matrix. If the matrix is already a +PMCMatrix2D, return a clone. =cut */ + METHOD convert_to_number_matrix() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_NumMatrix2D); const INTVAL totalsize = attrs->rows * attrs->cols; PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize")); @@ -736,7 +1005,7 @@ Get a PMCMatrix2D from the current matrix } METHOD convert_to_complex_matrix() { - DECLATTRS(SELF, attrs); + DECLATTRS_PMCMatrix2D(SELF, attrs); PMC * const d = Parrot_pmc_new(INTERP, __PLA_Type_ComplexMatrix2D); const INTVAL totalsize = attrs->rows * attrs->cols; PMC * const meth = VTABLE_find_method(INTERP, d, CONST_STRING(INTERP, "resize")); @@ -756,5 +1025,12 @@ Get a PMCMatrix2D from the current matrix RETURN(PMC * d); } +/* + +=back + +=cut + +*/ }