diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h index 562128b69..a86dac93f 100644 --- a/Eigen/src/SparseLU/SparseLU_pivotL.h +++ b/Eigen/src/SparseLU/SparseLU_pivotL.h @@ -71,7 +71,7 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScal // Determine the largest abs numerical value for partial pivoting Index diagind = iperm_c(jcol); // diagonal index - RealScalar pivmax = 0.0; + RealScalar pivmax(-1.0); Index pivptr = nsupc; Index diag = emptyIdxLU; RealScalar rtemp; @@ -87,8 +87,9 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScal } // Test for singularity - if ( pivmax == 0.0 ) { - pivrow = lsub_ptr[pivptr]; + if ( pivmax <= RealScalar(0.0) ) { + // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero + pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr]; perm_r(pivrow) = StorageIndex(jcol); return (jcol+1); } @@ -104,13 +105,13 @@ Index SparseLUImpl::pivotL(const Index jcol, const RealScal // Diagonal element exists using std::abs; rtemp = abs(lu_col_ptr[diag]); - if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag; + if (rtemp != RealScalar(0.0) && rtemp >= thresh) pivptr = diag; } pivrow = lsub_ptr[pivptr]; } // Record pivot row - perm_r(pivrow) = StorageIndex(jcol); + perm_r(pivrow) = StorageIndex(jcol); // Interchange row subscripts if (pivptr != nsupc ) { diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 64f6ac4fe..90362ca0f 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -332,7 +332,18 @@ Index generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, De return size; } -template void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000) + +struct prune_column { + Index m_col; + prune_column(Index col) : m_col(col) {} + template + bool operator()(Index, Index col, const Scalar&) const { + return col != m_col; + } +}; + + +template void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000, bool checkDeficient = false) { typedef typename Solver::MatrixType Mat; typedef typename Mat::Scalar Scalar; @@ -364,6 +375,13 @@ template void check_sparse_square_solving(Solver& solver, int m b = DenseVector::Zero(size); check_sparse_solving(solver, A, b, dA, b); } + // regression test for Bug 792 (structurally rank deficient matrices): + if(checkDeficient && size>1) { + Index col = internal::random(0,size-1); + A.prune(prune_column(col)); + solver.compute(A); + VERIFY_IS_EQUAL(solver.info(), NumericalIssue); + } } // First, get the folder diff --git a/test/sparselu.cpp b/test/sparselu.cpp index 231c857ad..c725847d8 100644 --- a/test/sparselu.cpp +++ b/test/sparselu.cpp @@ -42,8 +42,8 @@ template void test_sparselu_T() SparseLU, NaturalOrdering > sparselu_natural; check_sparse_square_solving(sparselu_colamd); - check_sparse_square_solving(sparselu_amd, 300, 2000); - check_sparse_square_solving(sparselu_natural, 300, 2000); + check_sparse_square_solving(sparselu_amd, 300, 2000, !true); // FIXME AMD ordering fails for structurally deficient matrices! + check_sparse_square_solving(sparselu_natural, 300, 2000, true); check_sparse_square_abs_determinant(sparselu_colamd); check_sparse_square_abs_determinant(sparselu_amd); diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h b/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h index 59ae4a2d0..85ae9dd6a 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h @@ -66,7 +66,7 @@ class BaseTensorContractionMapper { const bool left = (side == Lhs); Index nocontract_val = left ? row : col; Index linidx = 0; - for (int i = array_size::value - 1; i > 0; i--) { + for (int i = static_cast(array_size::value) - 1; i > 0; i--) { const Index idx = nocontract_val / m_ij_strides[i]; linidx += idx * m_nocontract_strides[i]; nocontract_val -= idx * m_ij_strides[i]; @@ -81,17 +81,19 @@ class BaseTensorContractionMapper { } Index contract_val = left ? col : row; - for (int i = array_size::value - 1; i > 0; i--) { + for (int i = static_cast(array_size::value) - 1; i > 0; i--) { const Index idx = contract_val / m_k_strides[i]; linidx += idx * m_contract_strides[i]; contract_val -= idx * m_k_strides[i]; } - EIGEN_STATIC_ASSERT(array_size::value > 0, YOU_MADE_A_PROGRAMMING_MISTAKE); - if (side == Rhs && inner_dim_contiguous) { - eigen_assert(m_contract_strides[0] == 1); - linidx += contract_val; - } else { - linidx += contract_val * m_contract_strides[0]; + + if(array_size::value > 0) { + if (side == Rhs && inner_dim_contiguous) { + eigen_assert(m_contract_strides[0] == 1); + linidx += contract_val; + } else { + linidx += contract_val * m_contract_strides[0]; + } } return linidx; @@ -102,7 +104,7 @@ class BaseTensorContractionMapper { const bool left = (side == Lhs); Index nocontract_val[2] = {left ? row : col, left ? row + distance : col}; Index linidx[2] = {0, 0}; - for (int i = array_size::value - 1; i > 0; i--) { + for (int i = static_cast(array_size::value) - 1; i > 0; i--) { const Index idx0 = nocontract_val[0] / m_ij_strides[i]; const Index idx1 = nocontract_val[1] / m_ij_strides[i]; linidx[0] += idx0 * m_nocontract_strides[i]; @@ -122,7 +124,7 @@ class BaseTensorContractionMapper { } Index contract_val[2] = {left ? col : row, left ? col : row + distance}; - for (int i = array_size::value - 1; i > 0; i--) { + for (int i = static_cast(array_size::value) - 1; i > 0; i--) { const Index idx0 = contract_val[0] / m_k_strides[i]; const Index idx1 = contract_val[1] / m_k_strides[i]; linidx[0] += idx0 * m_contract_strides[i]; @@ -130,7 +132,7 @@ class BaseTensorContractionMapper { contract_val[0] -= idx0 * m_k_strides[i]; contract_val[1] -= idx1 * m_k_strides[i]; } - EIGEN_STATIC_ASSERT(array_size::value > 0, YOU_MADE_A_PROGRAMMING_MISTAKE); + if (side == Rhs && inner_dim_contiguous) { eigen_assert(m_contract_strides[0] == 1); linidx[0] += contract_val[0]; @@ -509,8 +511,6 @@ struct TensorContractionEvaluatorBase static_cast(TensorEvaluator::Layout)), YOU_MADE_A_PROGRAMMING_MISTAKE); - eigen_assert((internal::array_size::value > 0) && "Must contract on some indices"); - DSizes eval_left_dims; DSizes eval_right_dims; @@ -558,7 +558,9 @@ struct TensorContractionEvaluatorBase m_i_strides[0] = 1; m_j_strides[0] = 1; - m_k_strides[0] = 1; + if(ContractDims) { + m_k_strides[0] = 1; + } m_i_size = 1; m_j_size = 1; diff --git a/unsupported/test/cxx11_tensor_contraction.cpp b/unsupported/test/cxx11_tensor_contraction.cpp index f4acdc504..b0d52c6cf 100644 --- a/unsupported/test/cxx11_tensor_contraction.cpp +++ b/unsupported/test/cxx11_tensor_contraction.cpp @@ -448,6 +448,31 @@ static void test_small_blocking_factors() } } +template +static void test_tensor_product() +{ + Tensor mat1(2, 3); + Tensor mat2(4, 1); + mat1.setRandom(); + mat2.setRandom(); + + Tensor result = mat1.contract(mat2, Eigen::array{{}}); + + VERIFY_IS_EQUAL(result.dimension(0), 2); + VERIFY_IS_EQUAL(result.dimension(1), 3); + VERIFY_IS_EQUAL(result.dimension(2), 4); + VERIFY_IS_EQUAL(result.dimension(3), 1); + for (int i = 0; i < result.dimension(0); ++i) { + for (int j = 0; j < result.dimension(1); ++j) { + for (int k = 0; k < result.dimension(2); ++k) { + for (int l = 0; l < result.dimension(3); ++l) { + VERIFY_IS_APPROX(result(i, j, k, l), mat1(i, j) * mat2(k, l) ); + } + } + } + } +} + void test_cxx11_tensor_contraction() { @@ -477,4 +502,6 @@ void test_cxx11_tensor_contraction() CALL_SUBTEST(test_tensor_vector()); CALL_SUBTEST(test_small_blocking_factors()); CALL_SUBTEST(test_small_blocking_factors()); + CALL_SUBTEST(test_tensor_product()); + CALL_SUBTEST(test_tensor_product()); }