00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED
00013 #define OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED
00014
00015 #include <Eigen/Sparse>
00016
00017 #include <algorithm>
00018 #include <iterator>
00019 #include <functional>
00020 #include <limits>
00021 #include <vector>
00022
00023 #include <Eigen/Core>
00024
00025 namespace Opm {
00026
00027 template < unsigned int depth >
00028 struct QuickSort
00029 {
00030 template <typename T>
00031 static inline void sort(T begin, T end)
00032 {
00033 if (begin != end)
00034 {
00035 T middle = std::partition (begin, end,
00036 std::bind2nd(std::less<typename std::iterator_traits<T>::value_type>(), *begin)
00037 );
00038 QuickSort< depth-1 >::sort(begin, middle);
00039
00040
00041 T new_middle = begin;
00042 QuickSort< depth-1 >::sort(++new_middle, end);
00043 }
00044 }
00045 };
00046
00047 template <>
00048 struct QuickSort< 0 >
00049 {
00050 template <typename T>
00051 static inline void sort(T begin, T end)
00052 {
00053
00054 std::sort( begin, end );
00055 }
00056 };
00057
00058
00059 template<typename Lhs, typename Rhs, typename ResultType>
00060 void fastSparseProduct(const Lhs& lhs, const Rhs& rhs, ResultType& res)
00061 {
00062
00063 res = ResultType(lhs.rows(), rhs.cols());
00064
00065
00066
00067 if( lhs.nonZeros() == 0 || rhs.nonZeros() == 0 )
00068 return;
00069
00070 typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
00071 typedef typename Eigen::internal::remove_all<Lhs>::type::Index Index;
00072
00073
00074 Index rows = lhs.innerSize();
00075 Index cols = rhs.outerSize();
00076 eigen_assert(lhs.outerSize() == rhs.innerSize());
00077
00078 std::vector<bool> mask(rows,false);
00079 Eigen::Matrix<Scalar,Eigen::Dynamic,1> values(rows);
00080 Eigen::Matrix<Index, Eigen::Dynamic,1> indices(rows);
00081
00082
00083
00084
00085
00086
00087
00088 Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros();
00089
00090 res.setZero();
00091 res.reserve(Index(estimated_nnz_prod));
00092
00093
00094 const Scalar epsilon = 0.0;
00095
00096
00097 for (Index j=0; j<cols; ++j)
00098 {
00099 Index nnz = 0;
00100 for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
00101 {
00102 const Scalar y = rhsIt.value();
00103 for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
00104 {
00105 const Scalar val = lhsIt.value() * y;
00106 if( std::abs( val ) > epsilon )
00107 {
00108 const Index i = lhsIt.index();
00109 if(!mask[i])
00110 {
00111 mask[i] = true;
00112 values[i] = val;
00113 indices[nnz] = i;
00114 ++nnz;
00115 }
00116 else
00117 values[i] += val;
00118 }
00119 }
00120 }
00121
00122 if( nnz > 1 )
00123 {
00124
00125 QuickSort< 1 >::sort( indices.data(), indices.data()+nnz );
00126 }
00127
00128 res.startVec(j);
00129
00130
00131 for(Index k=0; k<nnz; ++k)
00132 {
00133 const Index i = indices[k];
00134 res.insertBackByOuterInnerUnordered(j,i) = values[i];
00135 mask[i] = false;
00136 }
00137
00138 }
00139 res.finalize();
00140 }
00141
00142
00143
00144
00145 inline void fastDiagSparseProduct(const std::vector<double>& lhs,
00146 const Eigen::SparseMatrix<double>& rhs,
00147 Eigen::SparseMatrix<double>& res)
00148 {
00149 res = rhs;
00150
00151
00152 int n = res.cols();
00153 for (int col = 0; col < n; ++col) {
00154 typedef Eigen::SparseMatrix<double>::InnerIterator It;
00155 for (It it(res, col); it; ++it) {
00156 it.valueRef() *= lhs[it.row()];
00157 }
00158 }
00159 }
00160
00161
00162
00163
00164 inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
00165 const std::vector<double>& rhs,
00166 Eigen::SparseMatrix<double>& res)
00167 {
00168 res = lhs;
00169
00170
00171 int n = res.cols();
00172 for (int col = 0; col < n; ++col) {
00173 typedef Eigen::SparseMatrix<double>::InnerIterator It;
00174 for (It it(res, col); it; ++it) {
00175 it.valueRef() *= rhs[col];
00176 }
00177 }
00178 }
00179
00180 template<typename Lhs, typename Rhs>
00181 inline bool
00182 equalSparsityPattern(const Lhs& lhs, const Rhs& rhs)
00183 {
00184
00185 bool equal = (Lhs::IsRowMajor == Rhs::IsRowMajor) && (lhs.nonZeros() == rhs.nonZeros());
00186
00187
00188 if( equal )
00189 {
00190 typedef std::size_t Index;
00191 const Index outerSize = lhs.outerSize();
00192 const Index rhsOuterSize = rhs.outerSize();
00193 if( outerSize != rhsOuterSize )
00194 {
00195 return false;
00196 }
00197
00198
00199 const auto rhsOuter = rhs.outerIndexPtr();
00200 const auto lhsOuter = lhs.outerIndexPtr();
00201 for(Index i=0; i<=outerSize; ++i )
00202 {
00203 if( lhsOuter[ i ] != rhsOuter[ i ] ) {
00204 return false ;
00205 }
00206 }
00207
00208
00209 const auto rhsInner = rhs.innerIndexPtr();
00210 const auto lhsInner = lhs.innerIndexPtr();
00211
00212 const Index nnz = lhs.nonZeros();
00213 for( Index i=0; i<nnz; ++i)
00214 {
00215 if( lhsInner[ i ] != rhsInner[ i ] ) {
00216 return false;
00217 }
00218 }
00219 }
00220
00221 return equal;
00222 }
00223
00224
00225
00226 template<typename Lhs, typename Rhs>
00227 inline void
00228 fastSparseAdd(Lhs& lhs, const Rhs& rhs)
00229 {
00230 if( equalSparsityPattern( lhs, rhs ) )
00231 {
00232 typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
00233 typedef std::size_t Index;
00234
00235 const Index nnz = lhs.nonZeros();
00236
00237
00238 const Scalar* rhsV = rhs.valuePtr();
00239 Scalar* lhsV = lhs.valuePtr();
00240
00241 for(Index i=0; i<nnz; ++i )
00242 {
00243 lhsV[ i ] += rhsV[ i ];
00244 }
00245 }
00246 else
00247 {
00248
00249 lhs += rhs;
00250 }
00251 }
00252
00253
00254
00255 template<typename Lhs, typename Rhs>
00256 inline void
00257 fastSparseSubstract(Lhs& lhs, const Rhs& rhs)
00258 {
00259 if( equalSparsityPattern( lhs, rhs ) )
00260 {
00261 typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
00262 typedef std::size_t Index;
00263
00264 const Index nnz = lhs.nonZeros();
00265
00266
00267 const Scalar* rhsV = rhs.valuePtr();
00268 Scalar* lhsV = lhs.valuePtr();
00269
00270 for(Index i=0; i<nnz; ++i )
00271 {
00272 lhsV[ i ] -= rhsV[ i ];
00273 }
00274 }
00275 else
00276 {
00277
00278 lhs -= rhs;
00279 }
00280 }
00281
00282 }
00283
00284 #endif // OPM_FASTSPARSEPRODUCT_HEADER_INCLUDED