diff --git a/ext/nmatrix/math.cpp b/ext/nmatrix/math.cpp index a69c8f52..880c1cd4 100644 --- a/ext/nmatrix/math.cpp +++ b/ext/nmatrix/math.cpp @@ -188,7 +188,7 @@ extern "C" { // Math Functions // //////////////////// -namespace nm { +namespace nm { namespace math { /* @@ -232,8 +232,8 @@ namespace nm { DType* a = reinterpret_cast(storage->a); IType col_pos = storage->shape[0] + 1; if (M == 2) { - if (ija[2] - ija[0] == 2) { - *result = a[0] * a[1] - a[col_pos] * a[col_pos+1]; + if (ija[2] - ija[0] == 2) { + *result = a[0] * a[1] - a[col_pos] * a[col_pos+1]; } else { *result = a[0] * a[1]; } } else if (M == 3) { @@ -258,7 +258,7 @@ namespace nm { rb_raise(rb_eArgError, "some value in IJA is incorrect!"); } } - *result = + *result = m[0][0] * m[1][1] * m[2][2] + m[0][1] * m[1][2] * m[2][0] + m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] - m[0][2] * m[1][1] * m[2][0]; @@ -270,13 +270,13 @@ namespace nm { } /* - * Solve a system of linear equations using forward-substution followed by + * Solve a system of linear equations using forward-substution followed by * back substution from the LU factorization of the matrix of co-efficients. * Replaces x_elements with the result. Works only with non-integer, non-object * data types. * * args - r -> The number of rows of the matrix. - * lu_elements -> Elements of the LU decomposition of the co-efficients + * lu_elements -> Elements of the LU decomposition of the co-efficients * matrix, as a contiguos array. * b_elements -> Elements of the the right hand sides, as a contiguous array. * x_elements -> The array that will contain the results of the computation. @@ -291,7 +291,7 @@ namespace nm { const DType* b = reinterpret_cast(b_elements); DType* x = reinterpret_cast(x_elements); - for (int i = 0; i < r; ++i) { x[i] = b[i]; } + for (int i = 0; i < r; ++i) { x[i] = b[i]; } for (int i = 0; i < r; ++i) { // forward substitution loop ip = pivot[i]; sum = x[ip]; @@ -335,18 +335,18 @@ namespace nm { for (int row = k + 1; row < M; ++row) { typename MagnitudeDType::type big; big = magnitude( matrix[M*row + k] ); // element below the temp pivot - + if ( big > akk ) { interchange = row; - akk = big; + akk = big; } - } + } if (interchange != k) { // check if rows need flipping DType temp; for (int col = 0; col < M; ++col) { - NM_SWAP(matrix[interchange*M + col], matrix[k*M + col], temp); + NM_SWAP(matrix[interchange*M + col], matrix[k*M + col], temp); } } @@ -360,7 +360,7 @@ namespace nm { DType pivot = matrix[k * (M + 1)]; matrix[k * (M + 1)] = (DType)(1); // set diagonal as 1 for in-place inversion - for (int col = 0; col < M; ++col) { + for (int col = 0; col < M; ++col) { // divide each element in the kth row with the pivot matrix[k*M + col] = matrix[k*M + col] / pivot; } @@ -369,7 +369,7 @@ namespace nm { if (kk == k) continue; DType dum = matrix[k + M*kk]; - matrix[k + M*kk] = (DType)(0); // prepare for inplace inversion + matrix[k + M*kk] = (DType)(0); // prepare for inplace inversion for (int col = 0; col < M; ++col) { matrix[M*kk + col] = matrix[M*kk + col] - matrix[M*k + col] * dum; } @@ -384,7 +384,7 @@ namespace nm { for (int row = 0; row < M; ++row) { NM_SWAP(matrix[row * M + row_index[k]], matrix[row * M + col_index[k]], - temp); + temp); } } } @@ -410,14 +410,14 @@ namespace nm { DType sum_of_squares, *p_row, *psubdiag, *p_a, scale, innerproduct; int i, k, col; - // For each column use a Householder transformation to zero all entries + // For each column use a Householder transformation to zero all entries // below the subdiagonal. - for (psubdiag = a + nrows, col = 0; col < nrows - 2; psubdiag += nrows + 1, + for (psubdiag = a + nrows, col = 0; col < nrows - 2; psubdiag += nrows + 1, col++) { // Calculate the signed square root of the sum of squares of the // elements below the diagonal. - for (p_a = psubdiag, sum_of_squares = 0.0, i = col + 1; i < nrows; + for (p_a = psubdiag, sum_of_squares = 0.0, i = col + 1; i < nrows; p_a += nrows, i++) { sum_of_squares += *p_a * *p_a; } @@ -447,7 +447,7 @@ namespace nm { *p_a -= u[k] * innerproduct; } } - + // Postmultiply QA by Q for (p_row = a, i = 0; i < nrows; p_row += nrows, i++) { for (innerproduct = 0.0, k = col + 1; k < nrows; k++) { @@ -465,7 +465,7 @@ namespace nm { } void raise_not_invertible_error() { - rb_raise(nm_eNotInvertibleError, + rb_raise(nm_eNotInvertibleError, "matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)"); } @@ -473,7 +473,7 @@ namespace nm { * Calculate the exact inverse for a dense matrix (A [elements]) of size 2 or 3. Places the result in B_elements. */ template - void inverse_exact_from_dense(const int M, const void* A_elements, + void inverse_exact_from_dense(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb) { const DType* A = reinterpret_cast(A_elements); @@ -485,7 +485,7 @@ namespace nm { B[0] = A[lda+1] / det; B[1] = -A[1] / det; B[ldb] = -A[lda] / det; - B[ldb+1] = -A[0] / det; + B[ldb+1] = A[0] / det; } else if (M == 3) { // Calculate the exact determinant. @@ -510,7 +510,7 @@ namespace nm { } template - void inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, + void inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, YALE_STORAGE* inverse, const int ldb) { // inverse is a clone of storage @@ -524,18 +524,18 @@ namespace nm { if (M == 2) { IType ndnz = ija[2] - ija[0]; - if (ndnz == 2) { - det = a[0] * a[1] - a[col_pos] * a[col_pos+1]; + if (ndnz == 2) { + det = a[0] * a[1] - a[col_pos] * a[col_pos+1]; } else { det = a[0] * a[1]; } if (det == 0) { raise_not_invertible_error(); } b[0] = a[1] / det; - b[1] = -a[0] / det; - if (ndnz == 2) { + b[1] = a[0] / det; + if (ndnz == 2) { b[col_pos] = -a[col_pos] / det; b[col_pos+1] = -a[col_pos+1] / det; } - else if (ndnz == 1) { + else if (ndnz == 1) { b[col_pos] = -a[col_pos] / det; } @@ -561,7 +561,7 @@ namespace nm { rb_raise(rb_eArgError, "some value in IJA is incorrect!"); } } - det = + det = A[0] * A[lda+1] * A[2*lda+2] + A[1] * A[lda+2] * A[2*lda] + A[2] * A[lda] * A[2*lda+1] - A[0] * A[lda+2] * A[2*lda+1] - A[1] * A[lda] * A[2*lda+2] - A[2] * A[lda+1] * A[2*lda]; if (det == 0) { raise_not_invertible_error(); } @@ -1267,9 +1267,9 @@ static VALUE nm_clapack_laswp(VALUE self, VALUE n, VALUE a, VALUE lda, VALUE k1, /* * C accessor for calculating an exact determinant. Dense matrix version. */ -void nm_math_det_exact_from_dense(const int M, const void* elements, const int lda, +void nm_math_det_exact_from_dense(const int M, const void* elements, const int lda, nm::dtype_t dtype, void* result) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_dense, void, const int M, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_dense, void, const int M, const void* A_elements, const int lda, void* result_arg); ttable[dtype](M, elements, lda, result); @@ -1278,27 +1278,27 @@ void nm_math_det_exact_from_dense(const int M, const void* elements, const int l /* * C accessor for calculating an exact determinant. Yale matrix version. */ -void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, +void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, nm::dtype_t dtype, void* result) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_yale, void, const int M, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_yale, void, const int M, const YALE_STORAGE* storage, const int lda, void* result_arg); ttable[dtype](M, storage, lda, result); } /* - * C accessor for solving a system of linear equations. + * C accessor for solving a system of linear equations. */ void nm_math_solve(VALUE lu, VALUE b, VALUE x, VALUE ipiv) { int* pivot = new int[RARRAY_LEN(ipiv)]; for (int i = 0; i < RARRAY_LEN(ipiv); ++i) { - pivot[i] = FIX2INT(rb_ary_entry(ipiv, i)); + pivot[i] = FIX2INT(rb_ary_entry(ipiv, i)); } NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::solve, void, const int, const void*, const void*, void*, const int*); - ttable[NM_DTYPE(x)](NM_SHAPE0(b), NM_STORAGE_DENSE(lu)->elements, + ttable[NM_DTYPE(x)](NM_SHAPE0(b), NM_STORAGE_DENSE(lu)->elements, NM_STORAGE_DENSE(b)->elements, NM_STORAGE_DENSE(x)->elements, pivot); } @@ -1313,7 +1313,7 @@ void nm_math_hessenberg(VALUE a) { NULL, NULL, // does not support Complex NULL // no support for Ruby Object }; - + ttable[NM_DTYPE(a)](NM_SHAPE0(a), NM_STORAGE_DENSE(a)->elements); } /* @@ -1328,10 +1328,10 @@ void nm_math_inverse(const int M, void* a_elements, nm::dtype_t dtype) { /* * C accessor for calculating an exact inverse. Dense matrix version. */ -void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, +void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb, nm::dtype_t dtype) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_dense, void, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_dense, void, const int, const void*, const int, void*, const int); ttable[dtype](M, A_elements, lda, B_elements, ldb); @@ -1340,10 +1340,10 @@ void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, /* * C accessor for calculating an exact inverse. Yale matrix version. */ -void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, +void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, YALE_STORAGE* inverse, const int ldb, nm::dtype_t dtype) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_yale, void, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_yale, void, const int, const YALE_STORAGE*, const int, YALE_STORAGE*, const int); ttable[dtype](M, storage, lda, inverse, ldb); diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index d2f62710..1883da9d 100644 --- a/lib/nmatrix/math.rb +++ b/lib/nmatrix/math.rb @@ -112,6 +112,58 @@ def invert end alias :inverse :invert + # call-seq: + # exact_inverse! -> NMatrix + # + # Calulates inverse_exact of a matrix of size 2 or 3. + # Only works on dense matrices. + # + # * *Raises* : + # - +DataTypeError+ -> cannot invert an integer matrix in-place. + # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 # + def exact_inverse! + raise(ShapeError, "Cannot invert non-square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1] + raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype? + #No internal implementation of getri, so use this other function + n = self.shape[0] + if n>3 + raise(NotImplementedError, "Cannot find exact inverse of matrix of size greater than 3") + else + clond=self.clone + __inverse_exact__(clond, n, n) + end + end + + # + # call-seq: + # exact_inverse -> NMatrix + # + # Make a copy of the matrix, then invert using exact_inverse + # + # * *Returns* : + # - A dense NMatrix. Will be the same type as the input NMatrix, + # except if the input is an integral dtype, in which case it will be a + # :float64 NMatrix. + # + # * *Raises* : + # - +StorageTypeError+ -> only implemented on dense matrices. + # - +ShapeError+ -> matrix must be square. + # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 + # + def exact_inverse + #write this in terms of exact_inverse! so plugins will only have to overwrite + #exact_inverse! and not exact_inverse + if self.integer_dtype? + cloned = self.cast(dtype: :float64) + cloned.exact_inverse! + else + cloned = self.clone + cloned.exact_inverse! + end + end + alias :invert_exactly :exact_inverse + + # # call-seq: @@ -1458,4 +1510,4 @@ def dtype_for_floor_or_ceil self.__dense_map__ { |l| l.send(op,rhs) } end end -end \ No newline at end of file +end diff --git a/spec/00_nmatrix_spec.rb b/spec/00_nmatrix_spec.rb index 9ce6602b..09da20c4 100644 --- a/spec/00_nmatrix_spec.rb +++ b/spec/00_nmatrix_spec.rb @@ -73,7 +73,7 @@ e = NMatrix.new(2, [3,1,2,1], stype: :dense, dtype: :int64) inversed = e.method(:__inverse_exact__).call(e.clone, 2, 2) - f = NMatrix.new(2, [1,-1,-2,-3], stype: :dense, dtype: :int64) + f = NMatrix.new(2, [1,-1,-2,3], stype: :dense, dtype: :int64) expect(inversed).to eq(f) end @@ -91,7 +91,7 @@ e = NMatrix.new(2, [3,1,2,1], stype: :yale, dtype: :int64) inversed = e.method(:__inverse_exact__).call(e.clone, 2, 2) - f = NMatrix.new(2, [1,-1,-2,-3], stype: :yale, dtype: :int64) + f = NMatrix.new(2, [1,-1,-2,3], stype: :yale, dtype: :int64) expect(inversed).to eq(f) end @@ -104,7 +104,7 @@ c = NMatrix.new(2, [3,1,2,1], stype: :list, dtype: :int64) inversed = c.method(:__inverse_exact__).call(c.clone, 2, 2) - d = NMatrix.new(2, [1,-1,-2,-3], stype: :list, dtype: :int64) + d = NMatrix.new(2, [1,-1,-2,3], stype: :list, dtype: :int64) expect(inversed).to eq(d) end diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 14fd8e83..4276de6f 100644 --- a/spec/math_spec.rb +++ b/spec/math_spec.rb @@ -462,9 +462,9 @@ it "should correctly invert a matrix in place (bang)" do pending("not yet implemented for :object dtype") if dtype == :object - a = NMatrix.new(:dense, 5, [1, 8,-9, 7, 5, - 0, 1, 0, 4, 4, - 0, 0, 1, 2, 5, + a = NMatrix.new(:dense, 5, [1, 8,-9, 7, 5, + 0, 1, 0, 4, 4, + 0, 0, 1, 2, 5, 0, 0, 0, 1,-5, 0, 0, 0, 0, 1 ], dtype) b = NMatrix.new(:dense, 5, [1,-8, 9, 7, 17, @@ -495,6 +495,22 @@ expect(a.invert).to be_within(err).of(b) end + + it "should correctly find exact inverse" do + pending("not yet implemented for NMatrix-JRuby") if jruby? + a = NMatrix.new(:dense, 3, [1,2,3,0,1,4,5,6,0], dtype) + b = NMatrix.new(:dense, 3, [-24,18,5,20,-15,-4,-5,4,1], dtype) + + expect(a.exact_inverse).to be_within(err).of(b) + end + + it "should correctly find exact inverse" do + pending("not yet implemented for NMatrix-JRuby") if jruby? + a = NMatrix.new(:dense, 2, [1,3,3,8], dtype) + b = NMatrix.new(:dense, 2, [-8,3,3,-1], dtype) + + expect(a.exact_inverse).to be_within(err).of(b) + end end end @@ -1148,7 +1164,7 @@ 3, 4, 5, 6, 7, 8], stype: stype, dtype: dtype) end - + it "scales the matrix by a given factor and return the result" do pending("not yet implemented for :object dtype") if dtype == :object if integer_dtype? dtype @@ -1160,7 +1176,7 @@ 12, 14, 16], stype: stype, dtype: dtype)) end end - + it "scales the matrix in place by a given factor" do pending("not yet implemented for :object dtype") if dtype == :object if dtype == :int8