From c43b9f4a22f184dbbcdd1618bd5f277a2bac110d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Dec 2014 12:58:16 +0100 Subject: [PATCH] Use fastSparseProduct(), do not use hijacked Eigen header. --- opm/autodiff/AutoDiffBlock.hpp | 5 +- .../ConservativeSparseSparseProduct.h | 332 ------------------ 2 files changed, 3 insertions(+), 334 deletions(-) delete mode 100644 opm/autodiff/ConservativeSparseSparseProduct.h diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index de452b903..4326837f6 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -22,9 +22,9 @@ #include -#include #include #include +#include #include @@ -441,7 +441,8 @@ namespace Opm std::vector::M> jac(num_blocks); assert(lhs.cols() == rhs.value().rows()); for (int block = 0; block < num_blocks; ++block) { - jac[block] = lhs*rhs.derivative()[block]; + // jac[block] = lhs*rhs.derivative()[block]; + fastSparseProduct(lhs, rhs.derivative()[block], jac[block]); } typename AutoDiffBlock::V val = lhs*rhs.value().matrix(); return AutoDiffBlock::function(val, jac); diff --git a/opm/autodiff/ConservativeSparseSparseProduct.h b/opm/autodiff/ConservativeSparseSparseProduct.h deleted file mode 100644 index dd59b2038..000000000 --- a/opm/autodiff/ConservativeSparseSparseProduct.h +++ /dev/null @@ -1,332 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2011 Gael Guennebaud -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#ifndef EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H -#define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H - -#warning "Using overloaded Eigen::ConservativeSparseSparseProduct.h" - -#include -#include -#include -#include -#include - -#include - -namespace Eigen { - -// forward declaration of SparseMatrix -template -class SparseMatrix; - - -namespace internal { - -template < unsigned int depth > -struct QuickSort -{ - template - static inline void sort(T begin, T end) - { - if (begin != end) - { - T middle = std::partition (begin, end, - std::bind2nd(std::less::value_type>(), *begin) - ); - QuickSort< depth-1 >::sort(begin, middle); - - // std::sort (max(begin + 1, middle), end); - T new_middle = begin; - QuickSort< depth-1 >::sort(++new_middle, end); - } - } -}; - -template <> -struct QuickSort< 0 > -{ - template - static inline void sort(T begin, T end) - { - // fall back to standard insertion sort - std::sort( begin, end ); - } -}; - - -template -static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) -{ - // if one of the matrices does not contain non zero elements - // the result will only contain an empty matrix - if( lhs.nonZeros() == 0 || rhs.nonZeros() == 0 ) - return ; - - typedef typename remove_all::type::Scalar Scalar; - typedef typename remove_all::type::Index Index; - - // make sure to call innerSize/outerSize since we fake the storage order. - Index rows = lhs.innerSize(); - Index cols = rhs.outerSize(); - eigen_assert(lhs.outerSize() == rhs.innerSize()); - - std::vector mask(rows,false); - Matrix values(rows); - Matrix indices(rows); - - // estimate the number of non zero entries - // given a rhs column containing Y non zeros, we assume that the respective Y columns - // of the lhs differs in average of one non zeros, thus the number of non zeros for - // the product of a rhs column with the lhs is X+Y where X is the average number of non zero - // per column of the lhs. - // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs) - Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros(); - - res.setZero(); - res.reserve(Index(estimated_nnz_prod)); - - //const Scalar epsilon = std::numeric_limits< Scalar >::epsilon(); - const Scalar epsilon = 1e-15 ; - - // we compute each column of the result, one after the other - for (Index j=0; j epsilon ) - { - const Index i = lhsIt.index(); - if(!mask[i]) - { - mask[i] = true; - values[i] = val; - indices[nnz] = i; - ++nnz; - } - else - values[i] += val; - } - } - } - - if( nnz > 1 ) - { - // sort indices for sorted insertion to avoid later copying - QuickSort< 1 >::sort( indices.data(), indices.data()+nnz ); - } - - res.startVec(j); - // ordered insertion - // still using insertBackByOuterInnerUnordered since we know what we are doing - for(Index k=0; k use a quick sort - // otherwise => loop through the entire vector - // In order to avoid to perform an expensive log2 when the - // result is clearly very sparse we use a linear bound up to 200. - //if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); - for(Index k=0; k::Flags&RowMajorBit) ? RowMajor : ColMajor, - int RhsStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor, - int ResStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor> -struct conservative_sparse_sparse_product_selector; - -template -struct conservative_sparse_sparse_product_selector -{ - typedef typename remove_all::type LhsCleaned; - typedef typename LhsCleaned::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - //typedef SparseMatrix RowMajorMatrix; - typedef SparseMatrix ColMajorMatrix; - //ColMajorMatrix resCol(lhs.rows(),rhs.cols()); - res = ColMajorMatrix(lhs.rows(),rhs.cols()); - internal::conservative_sparse_sparse_product_impl(lhs, rhs, res); - //internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); - // sort the non zeros: - //RowMajorMatrix resRow(resCol); - //res = resRow; - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; - //RowMajorMatrix rhsRow = rhs; - //RowMajorMatrix resRow(lhs.rows(), rhs.cols()); - ColMajorMatrix lhsCol = lhs; - res = ResultType( lhs.rows(), rhs.cols() ); - internal::conservative_sparse_sparse_product_impl( lhsCol, rhs, res ); - //internal::conservative_sparse_sparse_product_impl(rhsRow, lhs, resRow); - //res = resRow; - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix rhsCol = rhs; - res = ResultType( lhs.rows(), rhs.cols() ); - internal::conservative_sparse_sparse_product_impl( lhs, rhsCol, res); - /* - typedef SparseMatrix RowMajorMatrix; - RowMajorMatrix lhsRow = lhs; - RowMajorMatrix resRow(lhs.rows(), rhs.cols()); - internal::conservative_sparse_sparse_product_impl(rhs, lhsRow, resRow); - res = resRow; - */ - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorMatrix; - RowMajorMatrix resRow(lhs.rows(), rhs.cols()); - internal::conservative_sparse_sparse_product_impl(rhs, lhs, resRow); - res = resRow; - } -}; - - -template -struct conservative_sparse_sparse_product_selector -{ - typedef typename traits::type>::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix resCol(lhs.rows(), rhs.cols()); - internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); - res = resCol; - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorMatrix; - RowMajorMatrix rhsRow = rhs; - res = ResultType( lhs.rows(), rhs.cols() ); - internal::conservative_sparse_sparse_product_impl(rhsRow, lhs, res); - /* - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix lhsCol = lhs; - ColMajorMatrix resCol(lhs.rows(), rhs.cols()); - internal::conservative_sparse_sparse_product_impl(lhsCol, rhs, resCol); - res = resCol; - */ - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorMatrix; - RowMajorMatrix lhsRow = lhs; - res = RowMajorMatrix( lhs.rows(), rhs.cols() ); - internal::conservative_sparse_sparse_product_impl(rhs, lhsRow, res); - /* - typedef SparseMatrix ColMajorMatrix; - ColMajorMatrix rhsCol = rhs; - ColMajorMatrix resCol(lhs.rows(), rhs.cols()); - internal::conservative_sparse_sparse_product_impl(lhs, rhsCol, resCol); - res = resCol; - */ - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorMatrix; - //typedef SparseMatrix ColMajorMatrix; - res = RowMajorMatrix( lhs.rows(),rhs.cols() ); - //RowMajorMatrix resRow(lhs.rows(),rhs.cols()); - internal::conservative_sparse_sparse_product_impl(rhs, lhs, res); - // sort the non zeros: - //ColMajorMatrix resCol(resRow); - //res = resCol; - } -}; - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H