diff --git a/ext/symengine/CMakeLists.txt b/ext/symengine/CMakeLists.txt index 2fe3b65..1b2a08a 100644 --- a/ext/symengine/CMakeLists.txt +++ b/ext/symengine/CMakeLists.txt @@ -12,6 +12,7 @@ set(RUBY_WRAPPER_SRC ruby_function.c ruby_ntheory.c ruby_utils.c + ruby_matrix.c symengine_utils.c symengine.c ) diff --git a/ext/symengine/ruby_matrix.c b/ext/symengine/ruby_matrix.c new file mode 100644 index 0000000..00303a1 --- /dev/null +++ b/ext/symengine/ruby_matrix.c @@ -0,0 +1,684 @@ +#include "ruby_matrix.h" + +void cmatrix_dense_free(void *ptr) +{ + CDenseMatrix *mat_ptr = ptr; + dense_matrix_free(mat_ptr); +} + +void cmatrix_sparse_free(void *ptr) +{ + CSparseMatrix *mat_ptr = ptr; + sparse_matrix_free(mat_ptr); +} + +VALUE cmatrix_dense_alloc(VALUE klass) +{ + CDenseMatrix *mat_ptr = dense_matrix_new(); + return Data_Wrap_Struct(klass, NULL, cmatrix_dense_free, mat_ptr); +} + +VALUE cmatrix_sparse_alloc(VALUE klass) +{ + CSparseMatrix *mat_ptr = sparse_matrix_new(); + return Data_Wrap_Struct(klass, NULL, cmatrix_sparse_free, mat_ptr); +} + +VALUE cmatrix_dense_to_str(VALUE self) +{ + CDenseMatrix *this; + char *str_ptr; + VALUE result; + + Data_Get_Struct(self, CDenseMatrix, this); + + str_ptr = dense_matrix_str(this); + result = rb_str_new_cstr(str_ptr); + basic_str_free(str_ptr); + + return result; +} + +VALUE cmatrix_sparse_to_str(VALUE self) +{ + CSparseMatrix *this; + char *str_ptr; + VALUE result; + + Data_Get_Struct(self, CSparseMatrix, this); + + str_ptr = sparse_matrix_str(this); + result = rb_str_new_cstr(str_ptr); + basic_str_free(str_ptr); + + return result; +} + +VALUE cmatrix_dense_init(VALUE self, VALUE args) +{ + int argc = RARRAY_LEN(args); + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + if (argc == 1) { + // SymEngine::DenseMatrix(SymEngine::DenseMatrix) OR + // SymEngine::DenseMatrix(NMatrix) OR + // SymEngine::DenseMatrix(Array) + + VALUE operand = rb_ary_shift(args); + char *s = rb_obj_classname(operand); + + if (strcmp(s, "SymEngine::DenseMatrix") == 0) { + + // SymEngine::DenseMatrix(SymEngine::DenseMatrix) + CDenseMatrix *temp; + Data_Get_Struct(operand, CDenseMatrix, temp); + dense_matrix_set(this, temp); + + } else if (strcmp(s, "NMatrix") == 0) { + + // SymEngine::DenseMatrix(NMatrix) + rb_raise(rb_eTypeError, "TODO"); + + } else if (strcmp(s, "Array") == 0) { + + // SymEngine::DenseMatrix(Array) + int counter = 0; + + int rows = RARRAY_LEN(operand); + int cols = -1; + CVecBasic *cargs = vecbasic_new(); + basic x; + basic_new_stack(x); + int i; + for (i = 0; i < rows; i++) { + int j; + VALUE row = rb_ary_shift(operand); + + s = rb_obj_classname(row); + if (strcmp(s, "Array") != 0) { + rb_raise(rb_eTypeError, "Not a 2D Array"); + } + + if (cols == -1) { + cols = RARRAY_LEN(row); + // Checking all rows for same col length + } else if (cols != RARRAY_LEN(row)) { + rb_raise(rb_eTypeError, + "2D Array's rows contain different column counts"); + } + + for (j = 0; j < cols; j++) { + sympify(rb_ary_shift(row), x); + vecbasic_push_back(cargs, x); + counter++; + } + } + this = dense_matrix_new_vec(rows, cols, cargs); + + basic_free_stack(x); + vecbasic_free(cargs); + + } else { + rb_raise(rb_eTypeError, "Invalid Arguments. A single argument" + "NMatrix, " + "or a single Array expected."); + } + + } else { + rb_raise(rb_eTypeError, "Invalid Arguments. A single argument" + "NMatrix, " + "or a single Array expected."); + } + + return self; +} + +VALUE cmatrix_sparse_init(VALUE self, VALUE args) +{ + int argc = RARRAY_LEN(args); + CSparseMatrix *this; + Data_Get_Struct(self, CSparseMatrix, this); + + if (argc == 0) { + + // SymEngine::SparseMatrix() + sparse_matrix_init(this); + + } else if (argc == 2) { + + // SymEngine::SparseMatrix(no_rows, no_cols) + VALUE val1 = rb_ary_shift(args); + VALUE val2 = rb_ary_shift(args); + + if ((TYPE(val1) == T_FIXNUM || TYPE(val1) == T_BIGNUM) + && (TYPE(val2) == T_FIXNUM || TYPE(val2) == T_BIGNUM)) { + + unsigned long int rows = NUM2ULONG(val1); + unsigned long int cols = NUM2ULONG(val2); + sparse_matrix_rows_cols(this, rows, cols); + + } else { + + rb_raise( + rb_eTypeError, + "Invalid Arguments. No Arguments, or two Numerics expected."); + } + } else { + + rb_raise(rb_eTypeError, + "Invalid Arguments. No Arguments, or two Numerics expected."); + } + + return self; +} + +VALUE cmatrix_dense_get(VALUE self, VALUE r, VALUE c) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + basic_struct *cresult; + VALUE result; + + unsigned long cbasic_r; + cbasic_r = NUM2ULONG(r); + unsigned long cbasic_c; + cbasic_c = NUM2ULONG(c); + + cresult = basic_new_heap(); + + dense_matrix_get_basic(cresult, this, cbasic_r, cbasic_c); + result = Data_Wrap_Struct(Klass_of_Basic(cresult), NULL, basic_free_heap, + cresult); + return result; +} + +VALUE cmatrix_dense_set(VALUE self, VALUE r, VALUE c, VALUE operand) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + basic cbasic_operand; + basic_new_stack(cbasic_operand); + sympify(operand, cbasic_operand); + + unsigned long cbasic_r; + cbasic_r = NUM2ULONG(r); + unsigned long cbasic_c; + cbasic_c = NUM2ULONG(c); + + dense_matrix_set_basic(this, cbasic_r, cbasic_c, cbasic_operand); + + basic_free_stack(cbasic_operand); + + return self; +} + +VALUE cmatrix_dense_det(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + basic_struct *cresult; + VALUE result; + + cresult = basic_new_heap(); + + dense_matrix_det(cresult, this); + result = Data_Wrap_Struct(Klass_of_Basic(cresult), NULL, basic_free_heap, + cresult); + return result; +} +VALUE cmatrix_dense_inv(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + + cresult = dense_matrix_new(); + + dense_matrix_inv(cresult, this); + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + return result; +} + +VALUE cmatrix_dense_transpose(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + + cresult = dense_matrix_new(); + + dense_matrix_transpose(cresult, this); + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + return result; +} + +VALUE cmatrix_dense_submatrix(VALUE self, VALUE r1, VALUE c1, VALUE r2, + VALUE c2, VALUE r_, VALUE c_) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + + unsigned long cbasic_r1; + cbasic_r1 = NUM2ULONG(r1); + unsigned long cbasic_c1; + cbasic_c1 = NUM2ULONG(c1); + + unsigned long cbasic_r2; + cbasic_r2 = NUM2ULONG(r2); + unsigned long cbasic_c2; + cbasic_c2 = NUM2ULONG(c2); + + unsigned long cbasic_r; + cbasic_r = NUM2ULONG(r_); + unsigned long cbasic_c; + cbasic_c = NUM2ULONG(c_); + + cresult = dense_matrix_new(); + + dense_matrix_submatrix(cresult, this, cbasic_r1, cbasic_c1, cbasic_r2, + cbasic_c2, cbasic_r, cbasic_c); + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + return result; +} + +VALUE cmatrix_dense_rows(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + return ULONG2NUM(dense_matrix_rows(this)); +} + +VALUE cmatrix_dense_cols(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + return ULONG2NUM(dense_matrix_cols(this)); +} + +VALUE cmatrix_dense_add(VALUE self, VALUE operand) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + + cresult = dense_matrix_new(); + + char *s = rb_obj_classname(operand); + + if (strcmp(s, "SymEngine::DenseMatrix") == 0) { + // Matrix Addition + CDenseMatrix *coperand; + Data_Get_Struct(operand, CDenseMatrix, coperand); + + dense_matrix_add_matrix(cresult, this, coperand); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, + cresult); + + } else { + // Scalar Addition + basic_struct *coperand = basic_new_heap(); + sympify(operand, coperand); + + dense_matrix_add_scalar(cresult, this, coperand); + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, + cresult); + basic_free_heap(coperand); + } + + return result; +} + +VALUE cmatrix_dense_mul(VALUE self, VALUE operand) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + + cresult = dense_matrix_new(); + + char *s = rb_obj_classname(operand); + + if (strcmp(s, "SymEngine::DenseMatrix") == 0) { + // Matrix Multiplication + CDenseMatrix *coperand; + Data_Get_Struct(operand, CDenseMatrix, coperand); + + dense_matrix_mul_matrix(cresult, this, coperand); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, + cresult); + } else { + // Scalar Multiplication + basic_struct *coperand = basic_new_heap(); + sympify(operand, coperand); + + dense_matrix_mul_scalar(cresult, this, coperand); + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, + cresult); + basic_free_heap(coperand); + } + + return result; +} + +VALUE cmatrix_dense_eq(VALUE self, VALUE operand) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + + cresult = dense_matrix_new(); + + char *s = rb_obj_classname(operand); + + if (strcmp(s, "SymEngine::DenseMatrix") == 0) { + CDenseMatrix *coperand; + Data_Get_Struct(operand, CDenseMatrix, coperand); + + result = dense_matrix_eq(this, coperand) == 1 ? Qtrue : Qfalse; + } else { + result = Qfalse; + } + return result; +} + +VALUE cmatrix_dense_neq(VALUE self, VALUE operand) +{ + VALUE result = cmatrix_dense_eq(self, operand); + rb_p(result); + return result == Qtrue ? Qfalse : Qtrue; +} + +VALUE cmatrix_dense_LU(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult_l; + VALUE result_l; + cresult_l = dense_matrix_new(); + + CDenseMatrix *cresult_u; + VALUE result_u; + cresult_u = dense_matrix_new(); + + dense_matrix_LU(cresult_l, cresult_u, this); + + result_l + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_l); + result_u + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_u); + + VALUE result; + result = rb_ary_new(); + rb_ary_push(result, result_l); + rb_ary_push(result, result_u); + + return result; +} + +VALUE cmatrix_dense_LDL(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult_l; + VALUE result_l; + cresult_l = dense_matrix_new(); + + CDenseMatrix *cresult_d; + VALUE result_d; + cresult_d = dense_matrix_new(); + + dense_matrix_LDL(cresult_l, cresult_d, this); + + result_l + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_l); + result_d + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_d); + + VALUE result; + result = rb_ary_new(); + rb_ary_push(result, result_l); + rb_ary_push(result, result_d); + + return result; +} + +VALUE cmatrix_dense_LU_solve(VALUE self, VALUE b) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult; + VALUE result; + cresult = dense_matrix_new(); + + CDenseMatrix *coperand; + Data_Get_Struct(b, CDenseMatrix, coperand); + + dense_matrix_LU_solve(cresult, this, coperand); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + return result; +} + +VALUE cmatrix_dense_FFLU(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult_lu; + VALUE result_lu; + cresult_lu = dense_matrix_new(); + + dense_matrix_FFLU(cresult_lu, this); + + result_lu + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_lu); + + return result_lu; +} + +VALUE cmatrix_dense_FFLDU(VALUE self) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CDenseMatrix, this); + + CDenseMatrix *cresult_l; + VALUE result_l; + cresult_l = dense_matrix_new(); + + CDenseMatrix *cresult_d; + VALUE result_d; + cresult_d = dense_matrix_new(); + + CDenseMatrix *cresult_u; + VALUE result_u; + cresult_u = dense_matrix_new(); + + dense_matrix_FFLDU(cresult_l, cresult_d, cresult_u, this); + + result_l + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_l); + result_d + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_d); + result_u + = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult_u); + + VALUE result; + result = rb_ary_new(); + rb_ary_push(result, result_l); + rb_ary_push(result, result_d); + rb_ary_push(result, result_u); + + return result; +} + +VALUE cmatrix_dense_ones(VALUE self, VALUE r, VALUE c) +{ + unsigned long int r_ = NUM2ULONG(r); + unsigned long int c_ = NUM2ULONG(c); + + CDenseMatrix *cresult; + cresult = dense_matrix_new(); + VALUE result; + + dense_matrix_ones(cresult, r_, c_); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + + return result; +} + +VALUE cmatrix_dense_zeros(VALUE self, VALUE r, VALUE c) +{ + unsigned long int r_ = NUM2ULONG(r); + unsigned long int c_ = NUM2ULONG(c); + + CDenseMatrix *cresult; + cresult = dense_matrix_new(); + VALUE result; + + dense_matrix_zeros(cresult, r_, c_); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + + return result; +} + +VALUE cmatrix_dense_diag(VALUE self, VALUE args) +{ + int argc = RARRAY_LEN(args); + + if (argc != 1 && argc != 2) { + rb_raise(rb_eTypeError, "Wrong number of arguments"); + } + + CDenseMatrix *cresult; + cresult = dense_matrix_new(); + VALUE result; + + VALUE operand = rb_ary_shift(args); + char *s = rb_obj_classname(operand); + + CVecBasic *cargs; + + if (strcmp(s, "Array") == 0) { + cargs = vecbasic_new(); + basic x; + basic_new_stack(x); + int cols = RARRAY_LEN(operand); + int j; + for (j = 0; j < cols; j++) { + sympify(rb_ary_shift(operand), x); + vecbasic_push_back(cargs, x); + } + basic_free_stack(x); + } else { + rb_raise(rb_eTypeError, "Invalid Argument type"); + } + + long int k = 0; + + if (argc == 2) { // With offset + k = NUM2LONG(rb_ary_shift(args)); + } + + dense_matrix_diag(cresult, cargs, k); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + vecbasic_free(cargs); + + return result; +} + +VALUE cmatrix_dense_eye(VALUE self, VALUE args) +{ + int argc = RARRAY_LEN(args); + + if (argc != 1 && argc != 3) { + rb_raise(rb_eTypeError, "Wrong number of arguments"); + } + + CDenseMatrix *cresult; + cresult = dense_matrix_new(); + VALUE result; + + VALUE operand = rb_ary_shift(args); + unsigned long int N = NUM2ULONG(operand); + unsigned long int M = N; + int k = 0; + + if (argc == 3) { + operand = rb_ary_shift(args); + M = NUM2ULONG(operand); + operand = rb_ary_shift(args); + k = NUM2INT(operand); + } + + dense_matrix_eye(cresult, N, M, k); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, cresult); + return result; +} + +VALUE cmatrix_sparse_get(VALUE self, VALUE r, VALUE c) +{ + CDenseMatrix *this; + Data_Get_Struct(self, CSparseMatrix, this); + + basic_struct *cresult; + VALUE result; + + unsigned long cbasic_r; + cbasic_r = NUM2ULONG(r); + unsigned long cbasic_c; + cbasic_c = NUM2ULONG(c); + + cresult = basic_new_heap(); + + sparse_matrix_get_basic(cresult, this, cbasic_r, cbasic_c); + result = Data_Wrap_Struct(Klass_of_Basic(cresult), NULL, basic_free_heap, + cresult); + return result; +} + +VALUE cmatrix_sparse_set(VALUE self, VALUE r, VALUE c, VALUE operand) +{ + CSparseMatrix *this; + Data_Get_Struct(self, CSparseMatrix, this); + + basic cbasic_operand; + basic_new_stack(cbasic_operand); + sympify(operand, cbasic_operand); + + unsigned long cbasic_r; + cbasic_r = NUM2ULONG(r); + unsigned long cbasic_c; + cbasic_c = NUM2ULONG(c); + + sparse_matrix_set_basic(this, cbasic_r, cbasic_c, cbasic_operand); + + basic_free_stack(cbasic_operand); + + return self; +} diff --git a/ext/symengine/ruby_matrix.h b/ext/symengine/ruby_matrix.h new file mode 100644 index 0000000..12e777e --- /dev/null +++ b/ext/symengine/ruby_matrix.h @@ -0,0 +1,45 @@ +#ifndef RUBY_MATRIX_H_ +#define RUBY_MATRIX_H_ + +#include +#include + +#include "symengine.h" +#include "ruby_basic.h" +#include "symengine_utils.h" + +void cmatrix_dense_free(void *ptr); +VALUE cmatrix_dense_alloc(VALUE klass); +VALUE cmatrix_dense_init(VALUE self, VALUE args); +VALUE cmatrix_dense_to_str(VALUE self); +VALUE cmatrix_dense_get(VALUE self, VALUE r, VALUE c); +VALUE cmatrix_dense_set(VALUE self, VALUE r, VALUE c, VALUE operand); +VALUE cmatrix_dense_det(VALUE self); +VALUE cmatrix_dense_inv(VALUE self); +VALUE cmatrix_dense_transpose(VALUE self); +VALUE cmatrix_dense_submatrix(VALUE self, VALUE r1, VALUE c1, VALUE r2, + VALUE c2, VALUE r_, VALUE c_); +VALUE cmatrix_dense_rows(VALUE self); +VALUE cmatrix_dense_cols(VALUE self); +VALUE cmatrix_dense_add(VALUE self, VALUE operand); +VALUE cmatrix_dense_mul(VALUE self, VALUE operand); +VALUE cmatrix_dense_eq(VALUE self, VALUE operand); +VALUE cmatrix_dense_neq(VALUE self, VALUE operand); +VALUE cmatrix_dense_LU(VALUE self); +VALUE cmatrix_dense_LDL(VALUE self); +VALUE cmatrix_dense_LU_solve(VALUE self, VALUE b); +VALUE cmatrix_dense_FFLU(VALUE self); +VALUE cmatrix_dense_FFLDU(VALUE self); +VALUE cmatrix_dense_ones(VALUE self, VALUE r, VALUE c); +VALUE cmatrix_dense_zeros(VALUE self, VALUE r, VALUE c); +VALUE cmatrix_dense_diag(VALUE self, VALUE args); +VALUE cmatrix_dense_eye(VALUE self, VALUE args); + +void cmatrix_sparse_free(void *ptr); +VALUE cmatrix_sparse_alloc(VALUE klass); +VALUE cmatrix_sparse_init(VALUE self, VALUE args); +VALUE cmatrix_sparse_to_str(VALUE self); +VALUE cmatrix_sparse_get(VALUE self, VALUE r, VALUE c); +VALUE cmatrix_sparse_set(VALUE self, VALUE r, VALUE c, VALUE operand); + +#endif // RUBY_MATRIX_H_ diff --git a/ext/symengine/ruby_utils.c b/ext/symengine/ruby_utils.c index 9973eb1..7f7d2e7 100644 --- a/ext/symengine/ruby_utils.c +++ b/ext/symengine/ruby_utils.c @@ -4,14 +4,59 @@ VALUE cutils_sympify(VALUE self, VALUE operand) { VALUE result; + const char *s = rb_obj_classname(operand); - basic_struct *cbasic_operand; - cbasic_operand = basic_new_heap(); + if (strcmp(s, "Array") == 0) { - sympify(operand, cbasic_operand); - result = Data_Wrap_Struct(Klass_of_Basic(cbasic_operand), NULL, - cbasic_free_heap, cbasic_operand); + CDenseMatrix *mat_result; + int counter = 0; + int rows = RARRAY_LEN(operand); + int cols = -1; + CVecBasic *cargs = vecbasic_new(); + basic x; + basic_new_stack(x); + int i; + + for (i = 0; i < rows; i++) { + int j; + VALUE row = rb_ary_shift(operand); + s = rb_obj_classname(row); + if (strcmp(s, "Array") != 0) { + rb_raise(rb_eTypeError, "2D Array is required"); + } + + if (cols == -1) { + cols = RARRAY_LEN(row); + // Checking all rows for same col length + } else if (cols != RARRAY_LEN(row)) { + rb_raise(rb_eTypeError, + "2D Array's rows contain different column counts"); + } + + for (j = 0; j < cols; j++) { + sympify(rb_ary_shift(row), x); + vecbasic_push_back(cargs, x); + counter++; + } + } + + mat_result = dense_matrix_new_vec(rows, cols, cargs); + + basic_free_stack(x); + vecbasic_free(cargs); + + result = Data_Wrap_Struct(c_dense_matrix, NULL, dense_matrix_free, + mat_result); + } else { + + basic_struct *cbasic_operand; + cbasic_operand = basic_new_heap(); + + sympify(operand, cbasic_operand); + result = Data_Wrap_Struct(Klass_of_Basic(cbasic_operand), NULL, + cbasic_free_heap, cbasic_operand); + } return result; } diff --git a/ext/symengine/ruby_utils.h b/ext/symengine/ruby_utils.h index 65edbf2..ec02e10 100644 --- a/ext/symengine/ruby_utils.h +++ b/ext/symengine/ruby_utils.h @@ -1,6 +1,7 @@ #ifndef RUBY_UTILS_H_ #define RUBY_UTILS_H_ +#include "symengine.h" #include "symengine_utils.h" // Returns the Ruby Value after going through sympify diff --git a/ext/symengine/symengine.c b/ext/symengine/symengine.c index 11af6a9..7c47c0c 100644 --- a/ext/symengine/symengine.c +++ b/ext/symengine/symengine.c @@ -12,6 +12,7 @@ #include "ruby_function.h" #include "ruby_ntheory.h" #include "ruby_utils.h" +#include "ruby_matrix.h" #include "symengine_utils.h" #include "symengine.h" @@ -241,5 +242,51 @@ void Init_symengine() rb_define_module_function(m_symengine, "lucas", cntheory_lucas, 1); rb_define_module_function(m_symengine, "binomial", cntheory_binomial, 2); + // MatrixBase Class + c_matrix_base + = rb_define_class_under(m_symengine, "MatrixBase", rb_cObject); + + // DenseMatrix and SparseMatrix Classes + c_dense_matrix + = rb_define_class_under(m_symengine, "DenseMatrix", c_matrix_base); + c_sparse_matrix + = rb_define_class_under(m_symengine, "SparseMatrix", c_matrix_base); + + // DenseMatrix Methods + rb_define_alloc_func(c_dense_matrix, cmatrix_dense_alloc); + rb_define_method(c_dense_matrix, "initialize", cmatrix_dense_init, -2); + rb_define_method(c_dense_matrix, "to_s", cmatrix_dense_to_str, 0); + rb_define_method(c_dense_matrix, "get", cmatrix_dense_get, 2); + rb_define_method(c_dense_matrix, "set", cmatrix_dense_set, 3); + rb_define_method(c_dense_matrix, "det", cmatrix_dense_det, 0); + rb_define_method(c_dense_matrix, "inv", cmatrix_dense_inv, 0); + rb_define_method(c_dense_matrix, "transpose", cmatrix_dense_transpose, 0); + rb_define_method(c_dense_matrix, "submatrix", cmatrix_dense_submatrix, 6); + rb_define_method(c_dense_matrix, "rows", cmatrix_dense_rows, 0); + rb_define_method(c_dense_matrix, "cols", cmatrix_dense_cols, 0); + rb_define_method(c_dense_matrix, "+", cmatrix_dense_add, 1); + rb_define_method(c_dense_matrix, "*", cmatrix_dense_mul, 1); + rb_define_method(c_dense_matrix, "==", cmatrix_dense_eq, 1); + rb_define_method(c_dense_matrix, "eql?", cmatrix_dense_eq, 1); + rb_define_method(c_dense_matrix, "!=", cmatrix_dense_neq, 1); + rb_define_method(c_dense_matrix, "LU", cmatrix_dense_LU, 0); + rb_define_method(c_dense_matrix, "LDL", cmatrix_dense_LDL, 0); + rb_define_method(c_dense_matrix, "LU_solve", cmatrix_dense_LU_solve, 1); + rb_define_method(c_dense_matrix, "FFLU", cmatrix_dense_FFLU, 0); + rb_define_method(c_dense_matrix, "FFLDU", cmatrix_dense_FFLDU, 0); + + // Numpy-like methods for DenseMatrix as module level methods + rb_define_module_function(m_symengine, "ones", cmatrix_dense_ones, 2); + rb_define_module_function(m_symengine, "zeros", cmatrix_dense_zeros, 2); + rb_define_module_function(m_symengine, "diag", cmatrix_dense_diag, -2); + rb_define_module_function(m_symengine, "eye", cmatrix_dense_eye, -2); + + // SparseMatrix Methods + rb_define_alloc_func(c_sparse_matrix, cmatrix_sparse_alloc); + rb_define_method(c_sparse_matrix, "initialize", cmatrix_sparse_init, -2); + rb_define_method(c_sparse_matrix, "to_s", cmatrix_sparse_to_str, 0); + rb_define_method(c_sparse_matrix, "get", cmatrix_sparse_get, 2); + rb_define_method(c_sparse_matrix, "set", cmatrix_sparse_set, 3); + symengine_print_stack_on_segfault(); } diff --git a/ext/symengine/symengine.h b/ext/symengine/symengine.h index 3f7b233..a144663 100644 --- a/ext/symengine/symengine.h +++ b/ext/symengine/symengine.h @@ -60,5 +60,8 @@ VALUE c_atanh; VALUE c_acsch; VALUE c_asech; VALUE c_acoth; +VALUE c_matrix_base; +VALUE c_dense_matrix; +VALUE c_sparse_matrix; #endif // SYMENGINE_H_ diff --git a/lib/symengine.rb b/lib/symengine.rb index 35ec80d..39a00cb 100644 --- a/lib/symengine.rb +++ b/lib/symengine.rb @@ -77,3 +77,5 @@ def lambdify_code(exp, syms) require 'symengine/complex' require 'symengine/complex_double' require 'symengine/undef_function' +require 'symengine/matrix_base' +require 'symengine/dense_matrix' diff --git a/lib/symengine/dense_matrix.rb b/lib/symengine/dense_matrix.rb new file mode 100644 index 0000000..e836df4 --- /dev/null +++ b/lib/symengine/dense_matrix.rb @@ -0,0 +1,7 @@ +module SymEngine + class DenseMatrix + def inspect + "#<#{self.class}(#{rows}x#{cols})>" + end + end +end diff --git a/lib/symengine/matrix_base.rb b/lib/symengine/matrix_base.rb new file mode 100644 index 0000000..4b871a3 --- /dev/null +++ b/lib/symengine/matrix_base.rb @@ -0,0 +1,7 @@ +module SymEngine + class MatrixBase + def inspect + "#<#{self.class}()>" + end + end +end diff --git a/spec/matrix_spec.rb b/spec/matrix_spec.rb new file mode 100644 index 0000000..494bb16 --- /dev/null +++ b/spec/matrix_spec.rb @@ -0,0 +1,66 @@ +describe SymEngine::DenseMatrix do + describe 'convert' do + subject { SymEngine([[1, 2],[3, 4]]) } + + it { is_expected.to be_a SymEngine::DenseMatrix } + its(:to_s) { is_expected.to eq "[1, 2]\n[3, 4]\n" } + end + + describe 'methods without arguments' do + subject { SymEngine([[4, 3],[3, 2]]) } + + its(:rows) { is_expected.to eq (2) } + its(:cols) { is_expected.to eq (2) } + + its(:inv) { is_expected.to be_a SymEngine::DenseMatrix } + its(:inv) { is_expected.to eq SymEngine([[-2, 3],[3, -4]]) } + + its(:FFLU) { is_expected.to be_a SymEngine::DenseMatrix } + its(:FFLU) { is_expected.to eq SymEngine([[4, 3],[3, -1]]) } + + its(:LU) { is_expected.to eq [SymEngine([[1, 0],[SymEngine(3)/SymEngine(4), 1]]), + SymEngine([[4, 3],[0, SymEngine(-1)/SymEngine(4)]])] } + + its(:LDL) { is_expected.to eq [SymEngine([[1, 0],[SymEngine(3)/SymEngine(4), 1]]), + SymEngine([[4, 0],[0, SymEngine(-1)/SymEngine(4)]])] } + + its(:FFLDU) { is_expected.to eq [SymEngine([[4, 0],[3, 1]]), SymEngine([[4, 0],[0, 4]]), + SymEngine([[4, 3],[0, -1]])] } + + its(:det) { is_expected.to eq (-1)} + + end + + describe 'methods with arguments' do + + let(:mat1) { SymEngine([[1, 2],[3, 4]]) } + let(:mat2) { SymEngine([[4, 3],[2, 1]]) } + let(:a) { SymEngine(4) } + let(:matA) { SymEngine([[3, 1, 2],[5, 7, 5], [1, 2, 3]]) } + let(:b) { SymEngine([[4],[4],[4]]) } + + it 'adds and multiplies' do + expect(mat1 + mat2).to eq(SymEngine([[5, 5],[5, 5]])) + expect(mat1 + a).to eq(SymEngine([[5, 6],[7, 8]])) + expect(mat1 * mat2).to eq(SymEngine([[8, 5],[20, 13]])) + expect(mat1 * a).to eq(SymEngine([[4, 8],[12, 16]])) + end + + it 'performs transpose and submatrix' do + expect(mat1.transpose).to eq(SymEngine([[1, 3],[2, 4]])) + expect(mat1.submatrix(0, 0, 1, 0, 1, 1)).to eq(SymEngine([[1],[3]])) + end + + it 'solves Ax = b' do + expect(matA.LU_solve(b)).to eq(SymEngine([[SymEngine(12)/SymEngine(29)], + [SymEngine(-32)/SymEngine(29)], + [SymEngine(56)/SymEngine(29)] + ])) + end + + it 'gets and sets elements' do + expect(mat1.set(0, 0, a)).to eq(SymEngine([[4, 2],[3, 4]])) + expect(mat1.get(1, 0)).to eq(SymEngine(3)) + end + end +end