diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-11-10 12:47:42 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-11-10 12:47:42 +0100 |
commit | 1879403562f453bd981ea45254865f51f6efc5c5 (patch) | |
tree | 6f3bb4b8881ed5c5a91964e9ae53a85a3ab53266 /Eigen | |
parent | 1333fe651d2b73df92cec8738097f893f698f468 (diff) |
mv the Skyline module to unsupported/
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/Skyline | 34 | ||||
-rw-r--r-- | Eigen/src/Core/util/Constants.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 3 | ||||
-rw-r--r-- | Eigen/src/Skyline/SkylineInplaceLU.h | 361 | ||||
-rw-r--r-- | Eigen/src/Skyline/SkylineMatrix.h | 870 | ||||
-rw-r--r-- | Eigen/src/Skyline/SkylineMatrixBase.h | 223 | ||||
-rw-r--r-- | Eigen/src/Skyline/SkylineProduct.h | 314 | ||||
-rw-r--r-- | Eigen/src/Skyline/SkylineStorage.h | 269 | ||||
-rw-r--r-- | Eigen/src/Skyline/SkylineUtil.h | 96 |
9 files changed, 1 insertions, 2171 deletions
diff --git a/Eigen/Skyline b/Eigen/Skyline deleted file mode 100644 index 5bcb87c76..000000000 --- a/Eigen/Skyline +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef EIGEN_SKYLINE_MODULE_H -#define EIGEN_SKYLINE_MODULE_H - - -#include "Eigen/Core" - -#include "Eigen/src/Core/util/DisableMSVCWarnings.h" - -#include <map> -#include <cstdlib> -#include <cstring> -#include <algorithm> - -namespace Eigen { - - /** \defgroup Skyline_Module Skyline module - * - * \nonstableyet - * - * - */ - -#include "src/Skyline/SkylineUtil.h" -#include "src/Skyline/SkylineMatrixBase.h" -#include "src/Skyline/SkylineStorage.h" -#include "src/Skyline/SkylineMatrix.h" -#include "src/Skyline/SkylineInplaceLU.h" -#include "src/Skyline/SkylineProduct.h" - -} // namespace Eigen - -#include "Eigen/src/Core/util/EnableMSVCWarnings.h" - -#endif // EIGEN_SKYLINE_MODULE_H diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 2ab427274..c9735b6e4 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -201,7 +201,7 @@ enum { AsRequested=0, EnforceAlignedAccess=2 }; enum { ConditionalJumpCost = 5 }; enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal, BothDirections }; -enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct, SkylineTimeDenseProduct }; +enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct }; enum { /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index af3ad8921..35e6dacf6 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -74,7 +74,6 @@ template<typename DecompositionType> struct ei_image_retval_base; template<typename DecompositionType> struct ei_image_retval; template<typename _Scalar, int Rows=Dynamic, int Cols=Dynamic, int Supers=Dynamic, int Subs=Dynamic, int Options=0> class BandMatrix; -template<typename _Scalar, int Rows=Dynamic, int Cols=Dynamic, int ColsUsed=Dynamic, int Options=0> class ColumnMatrix; template<typename Lhs, typename Rhs> struct ei_product_type; template<typename Lhs, typename Rhs, @@ -150,8 +149,6 @@ template<typename MatrixType,int Direction> class Homogeneous; // Sparse module: template<typename Lhs, typename Rhs, int ProductMode> class SparseProduct; -// Skyline module: -template<typename Lhs, typename Rhs, int ProductMode> class SkylineProduct; #endif // EIGEN_FORWARDDECLARATIONS_H diff --git a/Eigen/src/Skyline/SkylineInplaceLU.h b/Eigen/src/Skyline/SkylineInplaceLU.h deleted file mode 100644 index feed564c5..000000000 --- a/Eigen/src/Skyline/SkylineInplaceLU.h +++ /dev/null @@ -1,361 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SKYLINEINPLACELU_H -#define EIGEN_SKYLINEINPLACELU_H - -/** \ingroup Skyline_Module - * - * \class SkylineInplaceLU - * - * \brief Inplace LU decomposition of a skyline matrix and associated features - * - * \param MatrixType the type of the matrix of which we are computing the LU factorization - * - */ -template<typename MatrixType> -class SkylineInplaceLU { -protected: - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - -public: - - /** Creates a LU object and compute the respective factorization of \a matrix using - * flags \a flags. */ - SkylineInplaceLU(MatrixType& matrix, int flags = 0) - : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) { - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar > (); - m_lu.IsRowMajor ? computeRowMajor() : compute(); - } - - /** Sets the relative threshold value used to prune zero coefficients during the decomposition. - * - * Setting a value greater than zero speeds up computation, and yields to an imcomplete - * factorization with fewer non zero coefficients. Such approximate factors are especially - * useful to initialize an iterative solver. - * - * Note that the exact meaning of this parameter might depends on the actual - * backend. Moreover, not all backends support this feature. - * - * \sa precision() */ - void setPrecision(RealScalar v) { - m_precision = v; - } - - /** \returns the current precision. - * - * \sa setPrecision() */ - RealScalar precision() const { - return m_precision; - } - - /** Sets the flags. Possible values are: - * - CompleteFactorization - * - IncompleteFactorization - * - MemoryEfficient - * - one of the ordering methods - * - etc... - * - * \sa flags() */ - void setFlags(int f) { - m_flags = f; - } - - /** \returns the current flags */ - int flags() const { - return m_flags; - } - - void setOrderingMethod(int m) { - m_flags = m; - } - - int orderingMethod() const { - return m_flags; - } - - /** Computes/re-computes the LU factorization */ - void compute(); - void computeRowMajor(); - - /** \returns the lower triangular matrix L */ - //inline const MatrixType& matrixL() const { return m_matrixL; } - - /** \returns the upper triangular matrix U */ - //inline const MatrixType& matrixU() const { return m_matrixU; } - - template<typename BDerived, typename XDerived> - bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, - const int transposed = 0) const; - - /** \returns true if the factorization succeeded */ - inline bool succeeded(void) const { - return m_succeeded; - } - -protected: - RealScalar m_precision; - int m_flags; - mutable int m_status; - bool m_succeeded; - MatrixType& m_lu; -}; - -/** Computes / recomputes the in place LU decomposition of the SkylineInplaceLU. - * using the default algorithm. - */ -template<typename MatrixType> -//template<typename _Scalar> -void SkylineInplaceLU<MatrixType>::compute() { - const size_t rows = m_lu.rows(); - const size_t cols = m_lu.cols(); - - ei_assert(rows == cols && "We do not (yet) support rectangular LU."); - ei_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage"); - - for (unsigned int row = 0; row < rows; row++) { - const double pivot = m_lu.coeffDiag(row); - - //Lower matrix Columns update - const unsigned int& col = row; - for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) { - lIt.valueRef() /= pivot; - } - - //Upper matrix update -> contiguous memory access - typename MatrixType::InnerLowerIterator lIt(m_lu, col); - for (unsigned int rrow = row + 1; rrow < m_lu.rows(); rrow++) { - typename MatrixType::InnerUpperIterator uItPivot(m_lu, row); - typename MatrixType::InnerUpperIterator uIt(m_lu, rrow); - const double coef = lIt.value(); - - uItPivot += (rrow - row - 1); - - //update upper part -> contiguous memory access - for (++uItPivot; uIt && uItPivot;) { - uIt.valueRef() -= uItPivot.value() * coef; - - ++uIt; - ++uItPivot; - } - ++lIt; - } - - //Upper matrix update -> non contiguous memory access - typename MatrixType::InnerLowerIterator lIt3(m_lu, col); - for (unsigned int rrow = row + 1; rrow < m_lu.rows(); rrow++) { - typename MatrixType::InnerUpperIterator uItPivot(m_lu, row); - const double coef = lIt3.value(); - - //update lower part -> non contiguous memory access - for (unsigned int i = 0; i < rrow - row - 1; i++) { - m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef; - ++uItPivot; - } - ++lIt3; - } - //update diag -> contiguous - typename MatrixType::InnerLowerIterator lIt2(m_lu, col); - for (unsigned int rrow = row + 1; rrow < m_lu.rows(); rrow++) { - - typename MatrixType::InnerUpperIterator uItPivot(m_lu, row); - typename MatrixType::InnerUpperIterator uIt(m_lu, rrow); - const double coef = lIt2.value(); - - uItPivot += (rrow - row - 1); - m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef; - ++lIt2; - } - } -} - -template<typename MatrixType> -void SkylineInplaceLU<MatrixType>::computeRowMajor() { - const size_t rows = m_lu.rows(); - const size_t cols = m_lu.cols(); - - ei_assert(rows == cols && "We do not (yet) support rectangular LU."); - ei_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !"); - - for (unsigned int row = 0; row < rows; row++) { - typename MatrixType::InnerLowerIterator llIt(m_lu, row); - - - for (unsigned int col = llIt.col(); col < row; col++) { - if (m_lu.coeffExistLower(row, col)) { - const double diag = m_lu.coeffDiag(col); - - typename MatrixType::InnerLowerIterator lIt(m_lu, row); - typename MatrixType::InnerUpperIterator uIt(m_lu, col); - - - const int offset = lIt.col() - uIt.row(); - - - int stop = offset > 0 ? col - lIt.col() : col - uIt.row(); - - //#define VECTORIZE -#ifdef VECTORIZE - Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop); - Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop); - - - Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal); -#else - if (offset > 0) //Skip zero value of lIt - uIt += offset; - else //Skip zero values of uIt - lIt += -offset; - Scalar newCoeff = m_lu.coeffLower(row, col); - - for (int k = 0; k < stop; ++k) { - const Scalar tmp = newCoeff; - newCoeff = tmp - lIt.value() * uIt.value(); - ++lIt; - ++uIt; - } -#endif - - m_lu.coeffRefLower(row, col) = newCoeff / diag; - } - } - - //Upper matrix update - const int col = row; - typename MatrixType::InnerUpperIterator uuIt(m_lu, col); - for (unsigned int rrow = uuIt.row(); rrow < col; rrow++) { - - typename MatrixType::InnerLowerIterator lIt(m_lu, rrow); - typename MatrixType::InnerUpperIterator uIt(m_lu, col); - const int offset = lIt.col() - uIt.row(); - - int stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row(); - -#ifdef VECTORIZE - Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop); - Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop); - - Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal); -#else - if (offset > 0) //Skip zero value of lIt - uIt += offset; - else //Skip zero values of uIt - lIt += -offset; - Scalar newCoeff = m_lu.coeffUpper(rrow, col); - for (int k = 0; k < stop; ++k) { - const Scalar tmp = newCoeff; - newCoeff = tmp - lIt.value() * uIt.value(); - - ++lIt; - ++uIt; - } -#endif - m_lu.coeffRefUpper(rrow, col) = newCoeff; - } - - - //Diag matrix update - typename MatrixType::InnerLowerIterator lIt(m_lu, row); - typename MatrixType::InnerUpperIterator uIt(m_lu, row); - - const int offset = lIt.col() - uIt.row(); - - - int stop = offset > 0 ? lIt.size() : uIt.size(); -#ifdef VECTORIZE - Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop); - Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop); - Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal); -#else - if (offset > 0) //Skip zero value of lIt - uIt += offset; - else //Skip zero values of uIt - lIt += -offset; - Scalar newCoeff = m_lu.coeffDiag(row); - for (unsigned int k = 0; k < stop; ++k) { - const Scalar tmp = newCoeff; - newCoeff = tmp - lIt.value() * uIt.value(); - ++lIt; - ++uIt; - } -#endif - m_lu.coeffRefDiag(row) = newCoeff; - } -} - -/** Computes *x = U^-1 L^-1 b - * - * If \a transpose is set to SvTranspose or SvAdjoint, the solution - * of the transposed/adjoint system is computed instead. - * - * Not all backends implement the solution of the transposed or - * adjoint system. - */ -template<typename MatrixType> -template<typename BDerived, typename XDerived> -bool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const { - const size_t rows = m_lu.rows(); - const size_t cols = m_lu.cols(); - - - for (int row = 0; row < rows; row++) { - x->coeffRef(row) = b.coeff(row); - Scalar newVal = x->coeff(row); - typename MatrixType::InnerLowerIterator lIt(m_lu, row); - - unsigned int col = lIt.col(); - while (lIt.col() < row) { - - newVal -= x->coeff(col++) * lIt.value(); - ++lIt; - } - - x->coeffRef(row) = newVal; - } - - - for (int col = rows - 1; col > 0; col--) { - x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col); - - const Scalar x_col = x->coeff(col); - - typename MatrixType::InnerUpperIterator uIt(m_lu, col); - uIt += uIt.size()-1; - - - while (uIt) { - x->coeffRef(uIt.row()) -= x_col * uIt.value(); - //TODO : introduce --operator - uIt += -1; - } - - - } - x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0); - - return true; -} - -#endif // EIGEN_SKYLINELU_H diff --git a/Eigen/src/Skyline/SkylineMatrix.h b/Eigen/src/Skyline/SkylineMatrix.h deleted file mode 100644 index 8b3f385f6..000000000 --- a/Eigen/src/Skyline/SkylineMatrix.h +++ /dev/null @@ -1,870 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SKYLINEMATRIX_H -#define EIGEN_SKYLINEMATRIX_H - -#include "SkylineStorage.h" -#include "SkylineMatrixBase.h" - -/** \ingroup Skyline_Module - * - * \class SkylineMatrix - * - * \brief The main skyline matrix class - * - * This class implements a skyline matrix using the very uncommon storage - * scheme. - * - * \param _Scalar the scalar type, i.e. the type of the coefficients - * \param _Options Union of bit flags controlling the storage scheme. Currently the only possibility - * is RowMajor. The default is 0 which means column-major. - * - * - */ -template<typename _Scalar, int _Options> -struct ei_traits<SkylineMatrix<_Scalar, _Options> > { - typedef _Scalar Scalar; - - enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, - Flags = SkylineBit | _Options, - CoeffReadCost = NumTraits<Scalar>::ReadCost, - }; -}; - -template<typename _Scalar, int _Options> -class SkylineMatrix -: public SkylineMatrixBase<SkylineMatrix<_Scalar, _Options> > { -public: - EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(SkylineMatrix) - EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, +=) - EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, -=) - - using Base::IsRowMajor; - -protected: - - typedef SkylineMatrix<Scalar, (Flags&~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0) > TransposedSkylineMatrix; - - int m_outerSize; - int m_innerSize; - -public: - int* m_colStartIndex; - int* m_rowStartIndex; - SkylineStorage<Scalar> m_data; - -public: - - inline int rows() const { - return IsRowMajor ? m_outerSize : m_innerSize; - } - - inline int cols() const { - return IsRowMajor ? m_innerSize : m_outerSize; - } - - inline int innerSize() const { - return m_innerSize; - } - - inline int outerSize() const { - return m_outerSize; - } - - inline int upperNonZeros() const { - return m_data.upperSize(); - } - - inline int lowerNonZeros() const { - return m_data.lowerSize(); - } - - inline int upperNonZeros(int j) const { - return m_colStartIndex[j + 1] - m_colStartIndex[j]; - } - - inline int lowerNonZeros(int j) const { - return m_rowStartIndex[j + 1] - m_rowStartIndex[j]; - } - - inline const Scalar* _diagPtr() const { - return &m_data.diag(0); - } - - inline Scalar* _diagPtr() { - return &m_data.diag(0); - } - - inline const Scalar* _upperPtr() const { - return &m_data.upper(0); - } - - inline Scalar* _upperPtr() { - return &m_data.upper(0); - } - - inline const Scalar* _lowerPtr() const { - return &m_data.lower(0); - } - - inline Scalar* _lowerPtr() { - return &m_data.lower(0); - } - - inline const int* _upperProfilePtr() const { - return &m_data.upperProfile(0); - } - - inline int* _upperProfilePtr() { - return &m_data.upperProfile(0); - } - - inline const int* _lowerProfilePtr() const { - return &m_data.lowerProfile(0); - } - - inline int* _lowerProfilePtr() { - return &m_data.lowerProfile(0); - } - - inline Scalar coeff(int row, int col) const { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - - if (outer == inner) - return this->m_data.diag(outer); - - if (IsRowMajor) { - if (inner > outer) //upper matrix - { - const int minOuterIndex = inner - m_data.upperProfile(inner); - if (outer >= minOuterIndex) - return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); - else - return Scalar(0); - } - if (inner < outer) //lower matrix - { - const int minInnerIndex = outer - m_data.lowerProfile(outer); - if (inner >= minInnerIndex) - return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); - else - return Scalar(0); - } - return m_data.upper(m_colStartIndex[inner] + outer - inner); - } else { - if (outer > inner) //upper matrix - { - const int maxOuterIndex = inner + m_data.upperProfile(inner); - if (outer <= maxOuterIndex) - return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); - else - return Scalar(0); - } - if (outer < inner) //lower matrix - { - const int maxInnerIndex = outer + m_data.lowerProfile(outer); - - if (inner <= maxInnerIndex) - return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); - else - return Scalar(0); - } - } - } - - inline Scalar& coeffRef(int row, int col) { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - - if (outer == inner) - return this->m_data.diag(outer); - - if (IsRowMajor) { - if (col > row) //upper matrix - { - const int minOuterIndex = inner - m_data.upperProfile(inner); - ei_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); - } - if (col < row) //lower matrix - { - const int minInnerIndex = outer - m_data.lowerProfile(outer); - ei_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); - } - } else { - if (outer > inner) //upper matrix - { - const int maxOuterIndex = inner + m_data.upperProfile(inner); - ei_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); - } - if (outer < inner) //lower matrix - { - const int maxInnerIndex = outer + m_data.lowerProfile(outer); - ei_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); - } - } - } - - inline Scalar coeffDiag(int idx) const { - ei_assert(idx < outerSize()); - ei_assert(idx < innerSize()); - return this->m_data.diag(idx); - } - - inline Scalar coeffLower(int row, int col) const { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - ei_assert(inner != outer); - - if (IsRowMajor) { - const int minInnerIndex = outer - m_data.lowerProfile(outer); - if (inner >= minInnerIndex) - return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); - else - return Scalar(0); - - } else { - const int maxInnerIndex = outer + m_data.lowerProfile(outer); - if (inner <= maxInnerIndex) - return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); - else - return Scalar(0); - } - } - - inline Scalar coeffUpper(int row, int col) const { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - ei_assert(inner != outer); - - if (IsRowMajor) { - const int minOuterIndex = inner - m_data.upperProfile(inner); - if (outer >= minOuterIndex) - return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); - else - return Scalar(0); - } else { - const int maxOuterIndex = inner + m_data.upperProfile(inner); - if (outer <= maxOuterIndex) - return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); - else - return Scalar(0); - } - } - - inline Scalar& coeffRefDiag(int idx) { - ei_assert(idx < outerSize()); - ei_assert(idx < innerSize()); - return this->m_data.diag(idx); - } - - inline Scalar& coeffRefLower(int row, int col) { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - ei_assert(inner != outer); - - if (IsRowMajor) { - const int minInnerIndex = outer - m_data.lowerProfile(outer); - ei_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); - } else { - const int maxInnerIndex = outer + m_data.lowerProfile(outer); - ei_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); - } - } - - inline bool coeffExistLower(int row, int col) { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - ei_assert(inner != outer); - - if (IsRowMajor) { - const int minInnerIndex = outer - m_data.lowerProfile(outer); - return inner >= minInnerIndex; - } else { - const int maxInnerIndex = outer + m_data.lowerProfile(outer); - return inner <= maxInnerIndex; - } - } - - inline Scalar& coeffRefUpper(int row, int col) { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - ei_assert(inner != outer); - - if (IsRowMajor) { - const int minOuterIndex = inner - m_data.upperProfile(inner); - ei_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); - } else { - const int maxOuterIndex = inner + m_data.upperProfile(inner); - ei_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage"); - return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); - } - } - - inline bool coeffExistUpper(int row, int col) { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - ei_assert(inner != outer); - - if (IsRowMajor) { - const int minOuterIndex = inner - m_data.upperProfile(inner); - return outer >= minOuterIndex; - } else { - const int maxOuterIndex = inner + m_data.upperProfile(inner); - return outer <= maxOuterIndex; - } - } - - -protected: - -public: - class InnerUpperIterator; - class InnerLowerIterator; - - class OuterUpperIterator; - class OuterLowerIterator; - - /** Removes all non zeros */ - inline void setZero() { - m_data.clear(); - memset(m_colStartIndex, 0, (m_outerSize + 1) * sizeof (int)); - memset(m_rowStartIndex, 0, (m_outerSize + 1) * sizeof (int)); - } - - /** \returns the number of non zero coefficients */ - inline int nonZeros() const { - return m_data.diagSize() + m_data.upperSize() + m_data.lowerSize(); - } - - /** Preallocates \a reserveSize non zeros */ - inline void reserve(int reserveSize, int reserveUpperSize, int reserveLowerSize) { - m_data.reserve(reserveSize, reserveUpperSize, reserveLowerSize); - } - - /** \returns a reference to a novel non zero coefficient with coordinates \a row x \a col. - - * - * \warning This function can be extremely slow if the non zero coefficients - * are not inserted in a coherent order. - * - * After an insertion session, you should call the finalize() function. - */ - EIGEN_DONT_INLINE Scalar & insert(int row, int col) { - const int outer = IsRowMajor ? row : col; - const int inner = IsRowMajor ? col : row; - - ei_assert(outer < outerSize()); - ei_assert(inner < innerSize()); - - if (outer == inner) - return m_data.diag(col); - - if (IsRowMajor) { - if (outer < inner) //upper matrix - { - int minOuterIndex = 0; - minOuterIndex = inner - m_data.upperProfile(inner); - - if (outer < minOuterIndex) //The value does not yet exist - { - const int previousProfile = m_data.upperProfile(inner); - - m_data.upperProfile(inner) = inner - outer; - - - const int bandIncrement = m_data.upperProfile(inner) - previousProfile; - //shift data stored after this new one - const int stop = m_colStartIndex[cols()]; - const int start = m_colStartIndex[inner]; - - - for (int innerIdx = stop; innerIdx >= start; innerIdx--) { - m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx); - } - - for (int innerIdx = cols(); innerIdx > inner; innerIdx--) { - m_colStartIndex[innerIdx] += bandIncrement; - } - - //zeros new data - memset(this->_upperPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar)); - - return m_data.upper(m_colStartIndex[inner]); - } else { - return m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); - } - } - - if (outer > inner) //lower matrix - { - const int minInnerIndex = outer - m_data.lowerProfile(outer); - if (inner < minInnerIndex) //The value does not yet exist - { - const int previousProfile = m_data.lowerProfile(outer); - m_data.lowerProfile(outer) = outer - inner; - - const int bandIncrement = m_data.lowerProfile(outer) - previousProfile; - //shift data stored after this new one - const int stop = m_rowStartIndex[rows()]; - const int start = m_rowStartIndex[outer]; - - - for (int innerIdx = stop; innerIdx >= start; innerIdx--) { - m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx); - } - - for (int innerIdx = rows(); innerIdx > outer; innerIdx--) { - m_rowStartIndex[innerIdx] += bandIncrement; - } - - //zeros new data - memset(this->_lowerPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar)); - return m_data.lower(m_rowStartIndex[outer]); - } else { - return m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); - } - } - } else { - if (outer > inner) //upper matrix - { - const int maxOuterIndex = inner + m_data.upperProfile(inner); - if (outer > maxOuterIndex) //The value does not yet exist - { - const int previousProfile = m_data.upperProfile(inner); - m_data.upperProfile(inner) = outer - inner; - - const int bandIncrement = m_data.upperProfile(inner) - previousProfile; - //shift data stored after this new one - const int stop = m_rowStartIndex[rows()]; - const int start = m_rowStartIndex[inner + 1]; - - for (int innerIdx = stop; innerIdx >= start; innerIdx--) { - m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx); - } - - for (int innerIdx = inner + 1; innerIdx < outerSize() + 1; innerIdx++) { - m_rowStartIndex[innerIdx] += bandIncrement; - } - memset(this->_upperPtr() + m_rowStartIndex[inner] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar)); - return m_data.upper(m_rowStartIndex[inner] + m_data.upperProfile(inner)); - } else { - return m_data.upper(m_rowStartIndex[inner] + (outer - inner)); - } - } - - if (outer < inner) //lower matrix - { - const int maxInnerIndex = outer + m_data.lowerProfile(outer); - if (inner > maxInnerIndex) //The value does not yet exist - { - const int previousProfile = m_data.lowerProfile(outer); - m_data.lowerProfile(outer) = inner - outer; - - const int bandIncrement = m_data.lowerProfile(outer) - previousProfile; - //shift data stored after this new one - const int stop = m_colStartIndex[cols()]; - const int start = m_colStartIndex[outer + 1]; - - for (int innerIdx = stop; innerIdx >= start; innerIdx--) { - m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx); - } - - for (int innerIdx = outer + 1; innerIdx < outerSize() + 1; innerIdx++) { - m_colStartIndex[innerIdx] += bandIncrement; - } - memset(this->_lowerPtr() + m_colStartIndex[outer] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar)); - return m_data.lower(m_colStartIndex[outer] + m_data.lowerProfile(outer)); - } else { - return m_data.lower(m_colStartIndex[outer] + (inner - outer)); - } - } - } - } - - /** Must be called after inserting a set of non zero entries. - */ - inline void finalize() { - if (IsRowMajor) { - if (rows() > cols()) - m_data.resize(cols(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1); - else - m_data.resize(rows(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1); - - // ei_assert(rows() == cols() && "memory reorganisatrion only works with suare matrix"); - // - // Scalar* newArray = new Scalar[m_colStartIndex[cols()] + 1 + m_rowStartIndex[rows()] + 1]; - // unsigned int dataIdx = 0; - // for (unsigned int row = 0; row < rows(); row++) { - // - // const unsigned int nbLowerElts = m_rowStartIndex[row + 1] - m_rowStartIndex[row]; - // // std::cout << "nbLowerElts" << nbLowerElts << std::endl; - // memcpy(newArray + dataIdx, m_data.m_lower + m_rowStartIndex[row], nbLowerElts * sizeof (Scalar)); - // m_rowStartIndex[row] = dataIdx; - // dataIdx += nbLowerElts; - // - // const unsigned int nbUpperElts = m_colStartIndex[row + 1] - m_colStartIndex[row]; - // memcpy(newArray + dataIdx, m_data.m_upper + m_colStartIndex[row], nbUpperElts * sizeof (Scalar)); - // m_colStartIndex[row] = dataIdx; - // dataIdx += nbUpperElts; - // - // - // } - // //todo : don't access m_data profile directly : add an accessor from SkylineMatrix - // m_rowStartIndex[rows()] = m_rowStartIndex[rows()-1] + m_data.lowerProfile(rows()-1); - // m_colStartIndex[cols()] = m_colStartIndex[cols()-1] + m_data.upperProfile(cols()-1); - // - // delete[] m_data.m_lower; - // delete[] m_data.m_upper; - // - // m_data.m_lower = newArray; - // m_data.m_upper = newArray; - } else { - if (rows() > cols()) - m_data.resize(cols(), rows(), cols(), m_rowStartIndex[cols()] + 1, m_colStartIndex[cols()] + 1); - else - m_data.resize(rows(), rows(), cols(), m_rowStartIndex[rows()] + 1, m_colStartIndex[rows()] + 1); - } - } - - inline void squeeze() { - finalize(); - m_data.squeeze(); - } - - void prune(Scalar reference, RealScalar epsilon = precision<RealScalar > ()) { - //TODO - } - - /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero - * \sa resizeNonZeros(int), reserve(), setZero() - */ - void resize(size_t rows, size_t cols) { - const int diagSize = rows > cols ? cols : rows; - m_innerSize = IsRowMajor ? cols : rows; - - ei_assert(rows == cols && "Skyline matrix must be square matrix"); - - if (diagSize % 2) { // diagSize is odd - const int k = (diagSize - 1) / 2; - - m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols, - 2 * k * k + k + 1, - 2 * k * k + k + 1); - - } else // diagSize is even - { - const int k = diagSize / 2; - m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols, - 2 * k * k - k + 1, - 2 * k * k - k + 1); - } - - if (m_colStartIndex && m_rowStartIndex) { - delete[] m_colStartIndex; - delete[] m_rowStartIndex; - } - m_colStartIndex = new int [cols + 1]; - m_rowStartIndex = new int [rows + 1]; - m_outerSize = diagSize; - - m_data.reset(); - m_data.clear(); - - m_outerSize = diagSize; - memset(m_colStartIndex, 0, (cols + 1) * sizeof (int)); - memset(m_rowStartIndex, 0, (rows + 1) * sizeof (int)); - } - - void resizeNonZeros(int size) { - m_data.resize(size); - } - - inline SkylineMatrix() - : m_outerSize(-1), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { - resize(0, 0); - } - - inline SkylineMatrix(size_t rows, size_t cols) - : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { - resize(rows, cols); - } - - template<typename OtherDerived> - inline SkylineMatrix(const SkylineMatrixBase<OtherDerived>& other) - : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { - *this = other.derived(); - } - - inline SkylineMatrix(const SkylineMatrix & other) - : Base(), m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { - *this = other.derived(); - } - - inline void swap(SkylineMatrix & other) { - //EIGEN_DBG_SKYLINE(std::cout << "SkylineMatrix:: swap\n"); - std::swap(m_colStartIndex, other.m_colStartIndex); - std::swap(m_rowStartIndex, other.m_rowStartIndex); - std::swap(m_innerSize, other.m_innerSize); - std::swap(m_outerSize, other.m_outerSize); - m_data.swap(other.m_data); - } - - inline SkylineMatrix & operator=(const SkylineMatrix & other) { - std::cout << "SkylineMatrix& operator=(const SkylineMatrix& other)\n"; - if (other.isRValue()) { - swap(other.const_cast_derived()); - } else { - resize(other.rows(), other.cols()); - memcpy(m_colStartIndex, other.m_colStartIndex, (m_outerSize + 1) * sizeof (int)); - memcpy(m_rowStartIndex, other.m_rowStartIndex, (m_outerSize + 1) * sizeof (int)); - m_data = other.m_data; - } - return *this; - } - - template<typename OtherDerived> - inline SkylineMatrix & operator=(const SkylineMatrixBase<OtherDerived>& other) { - const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); - if (needToTranspose) { - // TODO - // return *this; - } else { - // there is no special optimization - return SkylineMatrixBase<SkylineMatrix>::operator=(other.derived()); - } - } - - friend std::ostream & operator <<(std::ostream & s, const SkylineMatrix & m) { - - EIGEN_DBG_SKYLINE( - std::cout << "upper elements : " << std::endl; - for (unsigned int i = 0; i < m.m_data.upperSize(); i++) - std::cout << m.m_data.upper(i) << "\t"; - std::cout << std::endl; - std::cout << "upper profile : " << std::endl; - for (unsigned int i = 0; i < m.m_data.upperProfileSize(); i++) - std::cout << m.m_data.upperProfile(i) << "\t"; - std::cout << std::endl; - std::cout << "lower startIdx : " << std::endl; - for (unsigned int i = 0; i < m.m_data.upperProfileSize(); i++) - std::cout << (IsRowMajor ? m.m_colStartIndex[i] : m.m_rowStartIndex[i]) << "\t"; - std::cout << std::endl; - - - std::cout << "lower elements : " << std::endl; - for (unsigned int i = 0; i < m.m_data.lowerSize(); i++) - std::cout << m.m_data.lower(i) << "\t"; - std::cout << std::endl; - std::cout << "lower profile : " << std::endl; - for (unsigned int i = 0; i < m.m_data.lowerProfileSize(); i++) - std::cout << m.m_data.lowerProfile(i) << "\t"; - std::cout << std::endl; - std::cout << "lower startIdx : " << std::endl; - for (unsigned int i = 0; i < m.m_data.lowerProfileSize(); i++) - std::cout << (IsRowMajor ? m.m_rowStartIndex[i] : m.m_colStartIndex[i]) << "\t"; - std::cout << std::endl; - ); - for (unsigned int rowIdx = 0; rowIdx < m.rows(); rowIdx++) { - for (unsigned int colIdx = 0; colIdx < m.cols(); colIdx++) { - s << m.coeff(rowIdx, colIdx) << "\t"; - } - s << std::endl; - } - return s; - } - - /** Destructor */ - inline ~SkylineMatrix() { - delete[] m_colStartIndex; - delete[] m_rowStartIndex; - } - - /** Overloaded for performance */ - Scalar sum() const; -}; - -template<typename Scalar, int _Options> -class SkylineMatrix<Scalar, _Options>::InnerUpperIterator { -public: - - InnerUpperIterator(const SkylineMatrix& mat, int outer) - : m_matrix(mat), m_outer(outer), - m_id(_Options == RowMajor ? mat.m_colStartIndex[outer] : mat.m_rowStartIndex[outer] + 1), - m_start(m_id), - m_end(_Options == RowMajor ? mat.m_colStartIndex[outer + 1] : mat.m_rowStartIndex[outer + 1] + 1) { - } - - inline InnerUpperIterator & operator++() { - m_id++; - return *this; - } - - inline InnerUpperIterator & operator+=(unsigned int shift) { - m_id += shift; - return *this; - } - - inline Scalar value() const { - return m_matrix.m_data.upper(m_id); - } - - inline Scalar* valuePtr() { - return const_cast<Scalar*> (&(m_matrix.m_data.upper(m_id))); - } - - inline Scalar& valueRef() { - return const_cast<Scalar&> (m_matrix.m_data.upper(m_id)); - } - - inline int index() const { - return IsRowMajor ? m_outer - m_matrix.m_data.upperProfile(m_outer) + (m_id - m_start) : - m_outer + (m_id - m_start) + 1; - } - - inline int row() const { - return IsRowMajor ? index() : m_outer; - } - - inline int col() const { - return IsRowMajor ? m_outer : index(); - } - - inline size_t size() const { - return m_matrix.m_data.upperProfile(m_outer); - } - - inline operator bool() const { - return (m_id < m_end) && (m_id >= m_start); - } - -protected: - const SkylineMatrix& m_matrix; - const int m_outer; - int m_id; - const int m_start; - const int m_end; -}; - -template<typename Scalar, int _Options> -class SkylineMatrix<Scalar, _Options>::InnerLowerIterator { -public: - - InnerLowerIterator(const SkylineMatrix& mat, int outer) - : m_matrix(mat), - m_outer(outer), - m_id(_Options == RowMajor ? mat.m_rowStartIndex[outer] : mat.m_colStartIndex[outer] + 1), - m_start(m_id), - m_end(_Options == RowMajor ? mat.m_rowStartIndex[outer + 1] : mat.m_colStartIndex[outer + 1] + 1) { - } - - inline InnerLowerIterator & operator++() { - m_id++; - return *this; - } - - inline InnerLowerIterator & operator+=(unsigned int shift) { - m_id += shift; - return *this; - } - - inline Scalar value() const { - return m_matrix.m_data.lower(m_id); - } - - inline Scalar* valuePtr() { - return const_cast<Scalar*> (&(m_matrix.m_data.lower(m_id))); - } - - inline Scalar& valueRef() { - return const_cast<Scalar&> (m_matrix.m_data.lower(m_id)); - } - - inline int index() const { - return IsRowMajor ? m_outer - m_matrix.m_data.lowerProfile(m_outer) + (m_id - m_start) : - m_outer + (m_id - m_start) + 1; - ; - } - - inline int row() const { - return IsRowMajor ? m_outer : index(); - } - - inline int col() const { - return IsRowMajor ? index() : m_outer; - } - - inline size_t size() const { - return m_matrix.m_data.lowerProfile(m_outer); - } - - inline operator bool() const { - return (m_id < m_end) && (m_id >= m_start); - } - -protected: - const SkylineMatrix& m_matrix; - const int m_outer; - int m_id; - const int m_start; - const int m_end; -}; - -#endif // EIGEN_SkylineMatrix_H diff --git a/Eigen/src/Skyline/SkylineMatrixBase.h b/Eigen/src/Skyline/SkylineMatrixBase.h deleted file mode 100644 index 3389a15c0..000000000 --- a/Eigen/src/Skyline/SkylineMatrixBase.h +++ /dev/null @@ -1,223 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SKYLINEMATRIXBASE_H -#define EIGEN_SKYLINEMATRIXBASE_H - -#include "SkylineUtil.h" - -/** \ingroup Skyline_Module - * - * \class SkylineMatrixBase - * - * \brief Base class of any skyline matrices or skyline expressions - * - * \param Derived - * - * - * - */ -template<typename Derived> class SkylineMatrixBase : public AnyMatrixBase<Derived> { -public: - - typedef typename ei_traits<Derived>::Scalar Scalar; - - enum { - RowsAtCompileTime = ei_traits<Derived>::RowsAtCompileTime, - /**< The number of rows at compile-time. This is just a copy of the value provided - * by the \a Derived type. If a value is not known at compile-time, - * it is set to the \a Dynamic constant. - * \sa MatrixBase::rows(), MatrixBase::cols(), ColsAtCompileTime, SizeAtCompileTime */ - - ColsAtCompileTime = ei_traits<Derived>::ColsAtCompileTime, - /**< The number of columns at compile-time. This is just a copy of the value provided - * by the \a Derived type. If a value is not known at compile-time, - * it is set to the \a Dynamic constant. - * \sa MatrixBase::rows(), MatrixBase::cols(), RowsAtCompileTime, SizeAtCompileTime */ - - - SizeAtCompileTime = (ei_size_at_compile_time<ei_traits<Derived>::RowsAtCompileTime, - ei_traits<Derived>::ColsAtCompileTime>::ret), - /**< This is equal to the number of coefficients, i.e. the number of - * rows times the number of columns, or to \a Dynamic if this is not - * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */ - - MaxRowsAtCompileTime = RowsAtCompileTime, - MaxColsAtCompileTime = ColsAtCompileTime, - - MaxSizeAtCompileTime = (ei_size_at_compile_time<MaxRowsAtCompileTime, - MaxColsAtCompileTime>::ret), - - IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1, - /**< This is set to true if either the number of rows or the number of - * columns is known at compile-time to be equal to 1. Indeed, in that case, - * we are dealing with a column-vector (if there is only one column) or with - * a row-vector (if there is only one row). */ - - Flags = ei_traits<Derived>::Flags, - /**< This stores expression \ref flags flags which may or may not be inherited by new expressions - * constructed from this one. See the \ref flags "list of flags". - */ - - CoeffReadCost = ei_traits<Derived>::CoeffReadCost, - /**< This is a rough measure of how expensive it is to read one coefficient from - * this expression. - */ - - IsRowMajor = Flags & RowMajorBit ? 1 : 0 - }; - -#ifndef EIGEN_PARSED_BY_DOXYGEN - /** This is the "real scalar" type; if the \a Scalar type is already real numbers - * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If - * \a Scalar is \a std::complex<T> then RealScalar is \a T. - * - * \sa class NumTraits - */ - typedef typename NumTraits<Scalar>::Real RealScalar; - - /** type of the equivalent square matrix */ - typedef Matrix<Scalar, EIGEN_ENUM_MAX(RowsAtCompileTime, ColsAtCompileTime), - EIGEN_ENUM_MAX(RowsAtCompileTime, ColsAtCompileTime) > SquareMatrixType; - - inline const Derived& derived() const { - return *static_cast<const Derived*> (this); - } - - inline Derived& derived() { - return *static_cast<Derived*> (this); - } - - inline Derived& const_cast_derived() const { - return *static_cast<Derived*> (const_cast<SkylineMatrixBase*> (this)); - } -#endif // not EIGEN_PARSED_BY_DOXYGEN - - /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ - inline int rows() const { - return derived().rows(); - } - - /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ - inline int cols() const { - return derived().cols(); - } - - /** \returns the number of coefficients, which is \a rows()*cols(). - * \sa rows(), cols(), SizeAtCompileTime. */ - inline int size() const { - return rows() * cols(); - } - - /** \returns the number of nonzero coefficients which is in practice the number - * of stored coefficients. */ - inline int nonZeros() const { - return derived().nonZeros(); - } - - /** \returns the size of the storage major dimension, - * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */ - int outerSize() const { - return (int(Flags) & RowMajorBit) ? this->rows() : this->cols(); - } - - /** \returns the size of the inner dimension according to the storage order, - * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */ - int innerSize() const { - return (int(Flags) & RowMajorBit) ? this->cols() : this->rows(); - } - - bool isRValue() const { - return m_isRValue; - } - - Derived& markAsRValue() { - m_isRValue = true; - return derived(); - } - - SkylineMatrixBase() : m_isRValue(false) { - /* TODO check flags */ - } - - inline Derived & operator=(const Derived& other) { - this->operator=<Derived > (other); - return derived(); - } - - template<typename OtherDerived> - inline void assignGeneric(const OtherDerived& other) { - derived().resize(other.rows(), other.cols()); - for (unsigned int row = 0; row < rows(); row++) - for (unsigned int col = 0; col < cols(); col++) { - if (other.coeff(row, col) != Scalar(0)) - derived().insert(row, col) = other.coeff(row, col); - } - derived().finalize(); - } - - template<typename OtherDerived> - inline Derived & operator=(const SkylineMatrixBase<OtherDerived>& other) { - //TODO - } - - template<typename Lhs, typename Rhs> - inline Derived & operator=(const SkylineProduct<Lhs, Rhs, SkylineTimeSkylineProduct>& product); - - friend std::ostream & operator <<(std::ostream & s, const SkylineMatrixBase& m) { - s << m.derived(); - return s; - } - - template<typename OtherDerived> - const typename SkylineProductReturnType<Derived, OtherDerived>::Type - operator*(const MatrixBase<OtherDerived> &other) const; - - /** \internal use operator= */ - template<typename DenseDerived> - void evalTo(MatrixBase<DenseDerived>& dst) const { - dst.setZero(); - for (unsigned int i = 0; i < rows(); i++) - for (unsigned int j = 0; j < rows(); j++) - dst(i, j) = derived().coeff(i, j); - } - - Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime> toDense() const { - return derived(); - } - - /** \returns the matrix or vector obtained by evaluating this expression. - * - * Notice that in the case of a plain matrix or vector (not an expression) this function just returns - * a const reference, in order to avoid a useless copy. - */ - EIGEN_STRONG_INLINE const typename ei_eval<Derived, IsSkyline>::type eval() const { - return typename ei_eval<Derived>::type(derived()); - } - -protected: - bool m_isRValue; -}; - -#endif // EIGEN_SkylineMatrixBase_H diff --git a/Eigen/src/Skyline/SkylineProduct.h b/Eigen/src/Skyline/SkylineProduct.h deleted file mode 100644 index 85ccacac8..000000000 --- a/Eigen/src/Skyline/SkylineProduct.h +++ /dev/null @@ -1,314 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SKYLINEPRODUCT_H -#define EIGEN_SKYLINEPRODUCT_H - -template<typename Lhs, typename Rhs> struct ei_skyline_product_mode { - - enum { - value = (Rhs::Flags & Lhs::Flags & SkylineBit) == SkylineBit - ? SkylineTimeSkylineProduct - : (Lhs::Flags & SkylineBit) == SkylineBit - ? SkylineTimeDenseProduct - : DenseTimeSkylineProduct - }; -}; - -template<typename Lhs, typename Rhs, int ProductMode> -struct SkylineProductReturnType { - typedef const typename ei_nested<Lhs, Rhs::RowsAtCompileTime>::type LhsNested; - typedef const typename ei_nested<Rhs, Lhs::RowsAtCompileTime>::type RhsNested; - - typedef SkylineProduct<LhsNested, RhsNested, ProductMode> Type; -}; - -template<typename LhsNested, typename RhsNested, int ProductMode> -struct ei_traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > { - // clean the nested types: - typedef typename ei_cleantype<LhsNested>::type _LhsNested; - typedef typename ei_cleantype<RhsNested>::type _RhsNested; - typedef typename _LhsNested::Scalar Scalar; - - enum { - LhsCoeffReadCost = _LhsNested::CoeffReadCost, - RhsCoeffReadCost = _RhsNested::CoeffReadCost, - LhsFlags = _LhsNested::Flags, - RhsFlags = _RhsNested::Flags, - - RowsAtCompileTime = _LhsNested::RowsAtCompileTime, - ColsAtCompileTime = _RhsNested::ColsAtCompileTime, - InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), - - MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime, - MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, - - EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), - ResultIsSkyline = ProductMode == SkylineTimeSkylineProduct, - - RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSkyline ? 0 : SkylineBit)), - - Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) - | EvalBeforeAssigningBit - | EvalBeforeNestingBit, - - CoeffReadCost = Dynamic - }; - - typedef typename ei_meta_if<ResultIsSkyline, - SkylineMatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> >, - MatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> > >::ret Base; -}; - -template<typename LhsNested, typename RhsNested, int ProductMode> -class SkylineProduct : ei_no_assignment_operator, -public ei_traits<SkylineProduct<LhsNested, RhsNested, ProductMode> >::Base { -public: - - EIGEN_GENERIC_PUBLIC_INTERFACE(SkylineProduct) - -private: - - typedef typename ei_traits<SkylineProduct>::_LhsNested _LhsNested; - typedef typename ei_traits<SkylineProduct>::_RhsNested _RhsNested; - -public: - - template<typename Lhs, typename Rhs> - EIGEN_STRONG_INLINE SkylineProduct(const Lhs& lhs, const Rhs& rhs) - : m_lhs(lhs), m_rhs(rhs) { - ei_assert(lhs.cols() == rhs.rows()); - - enum { - ProductIsValid = _LhsNested::ColsAtCompileTime == Dynamic - || _RhsNested::RowsAtCompileTime == Dynamic - || int(_LhsNested::ColsAtCompileTime) == int(_RhsNested::RowsAtCompileTime), - AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime, - SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested, _RhsNested) - }; - // note to the lost user: - // * for a dot product use: v1.dot(v2) - // * for a coeff-wise product use: v1.cwise()*v2 - EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), - INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) - EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), - INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) - EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) - } - - EIGEN_STRONG_INLINE int rows() const { - return m_lhs.rows(); - } - - EIGEN_STRONG_INLINE int cols() const { - return m_rhs.cols(); - } - - EIGEN_STRONG_INLINE const _LhsNested& lhs() const { - return m_lhs; - } - - EIGEN_STRONG_INLINE const _RhsNested& rhs() const { - return m_rhs; - } - -protected: - LhsNested m_lhs; - RhsNested m_rhs; -}; - -// dense = skyline * dense -// Note that here we force no inlining and separate the setZero() because GCC messes up otherwise - -template<typename Lhs, typename Rhs, typename Dest> -EIGEN_DONT_INLINE void ei_skyline_row_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) { - typedef typename ei_cleantype<Lhs>::type _Lhs; - typedef typename ei_cleantype<Rhs>::type _Rhs; - typedef typename ei_traits<Lhs>::Scalar Scalar; - - enum { - LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit, - LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit, - ProcessFirstHalf = LhsIsSelfAdjoint - && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0) - || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor) - || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)), - ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf) - }; - - //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix. - for (unsigned int col = 0; col < rhs.cols(); col++) { - for (unsigned int row = 0; row < lhs.rows(); row++) { - dst(row, col) = lhs.coeffDiag(row) * rhs(row, col); - } - } - //Use matrix lower triangular part - for (unsigned int row = 0; row < lhs.rows(); row++) { - typename _Lhs::InnerLowerIterator lIt(lhs, row); - const int stop = lIt.col() + lIt.size(); - for (unsigned int col = 0; col < rhs.cols(); col++) { - - unsigned int k = lIt.col(); - Scalar tmp = 0; - while (k < stop) { - tmp += - lIt.value() * - rhs(k++, col); - ++lIt; - } - dst(row, col) += tmp; - lIt += -lIt.size(); - } - - } - - //Use matrix upper triangular part - for (unsigned int lhscol = 0; lhscol < lhs.cols(); lhscol++) { - typename _Lhs::InnerUpperIterator uIt(lhs, lhscol); - const int stop = uIt.size() + uIt.row(); - for (unsigned int rhscol = 0; rhscol < rhs.cols(); rhscol++) { - - - const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol); - unsigned int k = uIt.row(); - while (k < stop) { - dst(k++, rhscol) += - uIt.value() * - rhsCoeff; - ++uIt; - } - uIt += -uIt.size(); - } - } - -} - -template<typename Lhs, typename Rhs, typename Dest> -EIGEN_DONT_INLINE void ei_skyline_col_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) { - typedef typename ei_cleantype<Lhs>::type _Lhs; - typedef typename ei_cleantype<Rhs>::type _Rhs; - typedef typename ei_traits<Lhs>::Scalar Scalar; - - enum { - LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit, - LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit, - ProcessFirstHalf = LhsIsSelfAdjoint - && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0) - || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor) - || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)), - ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf) - }; - - //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix. - for (unsigned int col = 0; col < rhs.cols(); col++) { - for (unsigned int row = 0; row < lhs.rows(); row++) { - dst(row, col) = lhs.coeffDiag(row) * rhs(row, col); - } - } - - //Use matrix upper triangular part - for (unsigned int row = 0; row < lhs.rows(); row++) { - typename _Lhs::InnerUpperIterator uIt(lhs, row); - const int stop = uIt.col() + uIt.size(); - for (unsigned int col = 0; col < rhs.cols(); col++) { - - unsigned int k = uIt.col(); - Scalar tmp = 0; - while (k < stop) { - tmp += - uIt.value() * - rhs(k++, col); - ++uIt; - } - - - dst(row, col) += tmp; - uIt += -uIt.size(); - } - } - - //Use matrix lower triangular part - for (unsigned int lhscol = 0; lhscol < lhs.cols(); lhscol++) { - typename _Lhs::InnerLowerIterator lIt(lhs, lhscol); - const int stop = lIt.size() + lIt.row(); - for (unsigned int rhscol = 0; rhscol < rhs.cols(); rhscol++) { - - const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol); - unsigned int k = lIt.row(); - while (k < stop) { - dst(k++, rhscol) += - lIt.value() * - rhsCoeff; - ++lIt; - } - lIt += -lIt.size(); - } - } - -} - -template<typename Lhs, typename Rhs, typename ResultType, - int LhsStorageOrder = ei_traits<Lhs>::Flags&RowMajorBit> - struct ei_skyline_product_selector; - -template<typename Lhs, typename Rhs, typename ResultType> -struct ei_skyline_product_selector<Lhs, Rhs, ResultType, RowMajor> { - typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) { - ei_skyline_row_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res); - } -}; - -template<typename Lhs, typename Rhs, typename ResultType> -struct ei_skyline_product_selector<Lhs, Rhs, ResultType, ColMajor> { - typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar; - - static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) { - ei_skyline_col_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res); - } -}; - -template<typename Derived> -template<typename Lhs, typename Rhs > -Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, SkylineTimeDenseProduct>& product) { - typedef typename ei_cleantype<Lhs>::type _Lhs; - ei_skyline_product_selector<typename ei_cleantype<Lhs>::type, - typename ei_cleantype<Rhs>::type, - Derived>::run(product.lhs(), product.rhs(), derived()); - - return derived(); -} - -// skyline * dense - -template<typename Derived> -template<typename OtherDerived > -EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived, OtherDerived>::Type -SkylineMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const { - - return typename SkylineProductReturnType<Derived, OtherDerived>::Type(derived(), other.derived()); -} - -#endif // EIGEN_SKYLINEPRODUCT_H diff --git a/Eigen/src/Skyline/SkylineStorage.h b/Eigen/src/Skyline/SkylineStorage.h deleted file mode 100644 index 97bdd3636..000000000 --- a/Eigen/src/Skyline/SkylineStorage.h +++ /dev/null @@ -1,269 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SKYLINE_STORAGE_H -#define EIGEN_SKYLINE_STORAGE_H - -/** Stores a skyline set of values in three structures : - * The diagonal elements - * The upper elements - * The lower elements - * - */ -template<typename Scalar> -class SkylineStorage { - typedef typename NumTraits<Scalar>::Real RealScalar; -public: - - SkylineStorage() - : m_diag(0), - m_lower(0), - m_upper(0), - m_lowerProfile(0), - m_upperProfile(0), - m_diagSize(0), - m_upperSize(0), - m_lowerSize(0), - m_upperProfileSize(0), - m_lowerProfileSize(0), - m_allocatedSize(0) { - } - - SkylineStorage(const SkylineStorage& other) - : m_diag(0), - m_lower(0), - m_upper(0), - m_lowerProfile(0), - m_upperProfile(0), - m_diagSize(0), - m_upperSize(0), - m_lowerSize(0), - m_upperProfileSize(0), - m_lowerProfileSize(0), - m_allocatedSize(0) { - *this = other; - } - - SkylineStorage & operator=(const SkylineStorage& other) { - resize(other.diagSize(), other.m_upperProfileSize, other.m_lowerProfileSize, other.upperSize(), other.lowerSize()); - memcpy(m_diag, other.m_diag, m_diagSize * sizeof (Scalar)); - memcpy(m_upper, other.m_upper, other.upperSize() * sizeof (Scalar)); - memcpy(m_lower, other.m_lower, other.lowerSize() * sizeof (Scalar)); - memcpy(m_upperProfile, other.m_upperProfile, m_upperProfileSize * sizeof (int)); - memcpy(m_lowerProfile, other.m_lowerProfile, m_lowerProfileSize * sizeof (int)); - return *this; - } - - void swap(SkylineStorage& other) { - std::swap(m_diag, other.m_diag); - std::swap(m_upper, other.m_upper); - std::swap(m_lower, other.m_lower); - std::swap(m_upperProfile, other.m_upperProfile); - std::swap(m_lowerProfile, other.m_lowerProfile); - std::swap(m_diagSize, other.m_diagSize); - std::swap(m_upperSize, other.m_upperSize); - std::swap(m_lowerSize, other.m_lowerSize); - std::swap(m_allocatedSize, other.m_allocatedSize); - } - - ~SkylineStorage() { - delete[] m_diag; - delete[] m_upper; - if (m_upper != m_lower) - delete[] m_lower; - delete[] m_upperProfile; - delete[] m_lowerProfile; - } - - void reserve(size_t size, size_t upperProfileSize, size_t lowerProfileSize, size_t upperSize, size_t lowerSize) { - int newAllocatedSize = size + upperSize + lowerSize; - if (newAllocatedSize > m_allocatedSize) - reallocate(size, upperProfileSize, lowerProfileSize, upperSize, lowerSize); - } - - void squeeze() { - if (m_allocatedSize > m_diagSize + m_upperSize + m_lowerSize) - reallocate(m_diagSize, m_upperProfileSize, m_lowerProfileSize, m_upperSize, m_lowerSize); - } - - void resize(size_t diagSize, size_t upperProfileSize, size_t lowerProfileSize, size_t upperSize, size_t lowerSize, float reserveSizeFactor = 0) { - if (m_allocatedSize < diagSize + upperSize + lowerSize) - reallocate(diagSize, upperProfileSize, lowerProfileSize, upperSize + size_t(reserveSizeFactor * upperSize), lowerSize + size_t(reserveSizeFactor * lowerSize)); - m_diagSize = diagSize; - m_upperSize = upperSize; - m_lowerSize = lowerSize; - m_upperProfileSize = upperProfileSize; - m_lowerProfileSize = lowerProfileSize; - } - - inline size_t diagSize() const { - return m_diagSize; - } - - inline size_t upperSize() const { - return m_upperSize; - } - - inline size_t lowerSize() const { - return m_lowerSize; - } - - inline size_t upperProfileSize() const { - return m_upperProfileSize; - } - - inline size_t lowerProfileSize() const { - return m_lowerProfileSize; - } - - inline size_t allocatedSize() const { - return m_allocatedSize; - } - - inline void clear() { - m_diagSize = 0; - } - - inline Scalar& diag(size_t i) { - return m_diag[i]; - } - - inline const Scalar& diag(size_t i) const { - return m_diag[i]; - } - - inline Scalar& upper(size_t i) { - return m_upper[i]; - } - - inline const Scalar& upper(size_t i) const { - return m_upper[i]; - } - - inline Scalar& lower(size_t i) { - return m_lower[i]; - } - - inline const Scalar& lower(size_t i) const { - return m_lower[i]; - } - - inline int& upperProfile(size_t i) { - return m_upperProfile[i]; - } - - inline const int& upperProfile(size_t i) const { - return m_upperProfile[i]; - } - - inline int& lowerProfile(size_t i) { - return m_lowerProfile[i]; - } - - inline const int& lowerProfile(size_t i) const { - return m_lowerProfile[i]; - } - - static SkylineStorage Map(int* upperProfile, int* lowerProfile, Scalar* diag, Scalar* upper, Scalar* lower, size_t size, size_t upperSize, size_t lowerSize) { - SkylineStorage res; - res.m_upperProfile = upperProfile; - res.m_lowerProfile = lowerProfile; - res.m_diag = diag; - res.m_upper = upper; - res.m_lower = lower; - res.m_allocatedSize = res.m_diagSize = size; - res.m_upperSize = upperSize; - res.m_lowerSize = lowerSize; - return res; - } - - inline void reset() { - memset(m_diag, 0, m_diagSize * sizeof (Scalar)); - memset(m_upper, 0, m_upperSize * sizeof (Scalar)); - memset(m_lower, 0, m_lowerSize * sizeof (Scalar)); - memset(m_upperProfile, 0, m_diagSize * sizeof (int)); - memset(m_lowerProfile, 0, m_diagSize * sizeof (int)); - } - - void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) { - //TODO - } - -protected: - - inline void reallocate(size_t diagSize, size_t upperProfileSize, size_t lowerProfileSize, size_t upperSize, size_t lowerSize) { - - Scalar* diag = new Scalar[diagSize]; - Scalar* upper = new Scalar[upperSize]; - Scalar* lower = new Scalar[lowerSize]; - int* upperProfile = new int[upperProfileSize]; - int* lowerProfile = new int[lowerProfileSize]; - - size_t copyDiagSize = std::min(diagSize, m_diagSize); - size_t copyUpperSize = std::min(upperSize, m_upperSize); - size_t copyLowerSize = std::min(lowerSize, m_lowerSize); - size_t copyUpperProfileSize = std::min(upperProfileSize, m_upperProfileSize); - size_t copyLowerProfileSize = std::min(lowerProfileSize, m_lowerProfileSize); - - // copy - memcpy(diag, m_diag, copyDiagSize * sizeof (Scalar)); - memcpy(upper, m_upper, copyUpperSize * sizeof (Scalar)); - memcpy(lower, m_lower, copyLowerSize * sizeof (Scalar)); - memcpy(upperProfile, m_upperProfile, copyUpperProfileSize * sizeof (int)); - memcpy(lowerProfile, m_lowerProfile, copyLowerProfileSize * sizeof (int)); - - - - // delete old stuff - delete[] m_diag; - delete[] m_upper; - delete[] m_lower; - delete[] m_upperProfile; - delete[] m_lowerProfile; - m_diag = diag; - m_upper = upper; - m_lower = lower; - m_upperProfile = upperProfile; - m_lowerProfile = lowerProfile; - m_allocatedSize = diagSize + upperSize + lowerSize; - m_upperSize = upperSize; - m_lowerSize = lowerSize; - } - -public: - Scalar* m_diag; - Scalar* m_upper; - Scalar* m_lower; - int* m_upperProfile; - int* m_lowerProfile; - size_t m_diagSize; - size_t m_upperSize; - size_t m_lowerSize; - size_t m_upperProfileSize; - size_t m_lowerProfileSize; - size_t m_allocatedSize; - -}; - -#endif // EIGEN_COMPRESSED_STORAGE_H diff --git a/Eigen/src/Skyline/SkylineUtil.h b/Eigen/src/Skyline/SkylineUtil.h deleted file mode 100644 index 949508994..000000000 --- a/Eigen/src/Skyline/SkylineUtil.h +++ /dev/null @@ -1,96 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Guillaume Saupin <guillaume.saupin@cea.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SKYLINEUTIL_H -#define EIGEN_SKYLINEUTIL_H - -#ifdef NDEBUG -#define EIGEN_DBG_SKYLINE(X) -#else -#define EIGEN_DBG_SKYLINE(X) X -#endif - -const unsigned int SkylineBit = 0x1200; -template<typename Lhs, typename Rhs, int ProductMode> class SkylineProduct; -enum AdditionalProductEvaluationMode {SkylineTimeSkylineProduct, DenseTimeSkylineProduct}; -enum {IsSkyline = SkylineBit}; - - -#define EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \ -template<typename OtherDerived> \ -EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::SkylineMatrixBase<OtherDerived>& other) \ -{ \ - return Base::operator Op(other.derived()); \ -} \ -EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \ -{ \ - return Base::operator Op(other); \ -} - -#define EIGEN_SKYLINE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \ -template<typename Other> \ -EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \ -{ \ - return Base::operator Op(scalar); \ -} - -#define EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ -EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \ -EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \ -EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \ -EIGEN_SKYLINE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \ -EIGEN_SKYLINE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=) - -#define _EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \ -typedef BaseClass Base; \ -typedef typename Eigen::ei_traits<Derived>::Scalar Scalar; \ -typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \ -enum { Flags = Eigen::ei_traits<Derived>::Flags, }; - -#define EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(Derived) \ -_EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::SkylineMatrixBase<Derived>) - -template<typename Derived> class SkylineMatrixBase; -template<typename _Scalar, int _Flags = 0> class SkylineMatrix; -template<typename _Scalar, int _Flags = 0> class DynamicSkylineMatrix; -template<typename _Scalar, int _Flags = 0> class SkylineVector; -template<typename _Scalar, int _Flags = 0> class MappedSkylineMatrix; - -template<typename Lhs, typename Rhs> struct ei_skyline_product_mode; -template<typename Lhs, typename Rhs, int ProductMode = ei_skyline_product_mode<Lhs,Rhs>::value> struct SkylineProductReturnType; - - -template<typename T> class ei_eval<T,IsSkyline> -{ - typedef typename ei_traits<T>::Scalar _Scalar; - enum { - _Flags = ei_traits<T>::Flags - }; - - public: - typedef SkylineMatrix<_Scalar, _Flags> type; -}; - - -#endif // EIGEN_SKYLINEUTIL_H |