| Re: [eigen] Skyline matrix |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
You'll find enclosed the code as a patch.
guillaume saupin wrote:
Hello !
I've finally (after some performances tests) implemented a basic
skyline matrix using the following storage :
- an array stores the diagonal elements
- an array stores the upper elements in a column oriented way
- an array stores the lower elements in a row oriented way
I've implemented a LU decomposition for this format. I use the
following algorithm (with minor modification to handle the sparseness
of the matrix):
for each row
-- compute the lower part
for each col < row
LU(row, col) = (LU(row, col) - sum(k = 0 ; k < col) [LU(row,
k) * LU(k, col)] ) / LU(row, row) -- elements are access contiguously
/ Can be easily vectorize
-- compute the upper part
for each rrow < row
LU(rrow, row) -= sum(k = 0 ; k < rrow) [LU(rrow, k) * LU(k, row)]
-- elements are access contiguously / Can be easily vectorize
-- compute the diag part
LU(row, row) -= sum(k = 0 ; k < row) [LU(rrow, k) * LU(k, row)] --
elements are access contiguously / Can be easily vectorize
This algorithm is well suited for our storage format. After a few
tests, it seems even more efficient to store the the upper and lower
elements in the same array, with the elements from one row followed by
the elements of the cols of same index, and so on ...
I've also implemented a product for the Skyline * Dense operation.
Do you think that it would be possible to insert this code into eigen ?
guillaume saupin wrote:
Gael Guennebaud a écrit :
On Tue, Oct 27, 2009 at 10:40 AM, guillaume saupin
<guillaume.saupin@xxxxxx <mailto:guillaume.saupin@xxxxxx>> wrote:
Gael Guennebaud a écrit :
this is not an easy question because we are still unsure how
to manage efficiently and with minimal code all these kind of
"special" matrices.
Of course an easy solution would be to make it inherit
SparseMatrixBase and let the generic algorithms based on
sparse iterators do the job. However, this approach won't be
optimal because the iterators won't take advantage of the
specificity of a SkylineMatrix which, e.g., should allow
vectorization.
Actually skyline matrices are very similar to banded matrices,
and in particular it seems to me that they are similar enough
to share the same algorithms. So for instance these two kind
of matrices could inherit the same matrix base and the
algorithms could be implemented in a generic way via the
concept of "range inner vector" which is a column (or a row)
vector with an associated start index... We also need
efficient basic operators dealing with such vector. E.g. the
addition of two such vector can be efficiently done by first
copying the non overlapping part (start and end), and then sum
the overlapping part using Eigen's core module feature.
The way we implemented skyline matrix is not really similar to
what is usually done for banded matrix, as we store the diagonal,
the upper, and the lower elements in separate arrays. This way, we
can achieve quite efficient LU decomposition (and solving) due to
coherent memory accesses.
ok I see, then it is not really what people call skyline storage
(http://www.netlib.org/linalg/html_templates/node96.html), but it is
more like an extension of the tridiagonal storage.
I think that our storage is actually quite similar to the one
proposed there for non symetric skyline (referred as NSK format by
Saad) matrix. As in these format, we store the upper element in a
column-oriented way, and the lower element in a row-oriented way.
The only difference is that we store the diagonal elements in a
separate array instead of keeping an index array to these elements.
This allows contiguous access to diagonal elements.
However, I don't see how your storage can yield to a more efficient
implementation of the LU decomposition than the standard skyline
storage. For instance here is a basic algorithm without pivoting for
the standard skyline storage:
for(int k = 0; k+1 < rows; ++k)
{
int rrows = rows-k-1;
int rsize = size-k-1;
// line 1
lu.col(k).end(rrows) /= lu.coeff(k,k);
// line 2
lu.corner(BottomRight,rrows,rsize).noalias() -=
lu.col(k).end(rrows) * lu.row(k).end(rsize);
}
where lu is the working matrix, lu.col(k) is assumed to return a
"range vector" as I described in my previous email.
Here line 1 would be trivially optimized (i.e., vectorized) since
lu.col(k).end(rrows) is just a small dense vector.
Line 2 is an outer product which again is trivially/automatically
vectorized sicne it is impelmented as a sequence of:
"col_range_vector_i -= scalar * col_range_vector".
Here the locality is pretty good because the vector
"lu.col(k).end(rrows)" which is reused multiple times is
sequentially stored in memory.
But perhaps there exists a special algorithm which perfectly fit
your storage ? Is it possible to see your code somewhere ? Finally,
if it is really more efficient then that would make sense to have it
in Eigen.
The benefit of using the NSK, i.e. of storing the lower matrix in a
column-oriented way, and the upper matrix in a row-oriented way is
that in line 1 we access contigous elements (the column k of the
lower matrix), and in line 2 we also access contiguous element (the
row elements of the upper matrix). I don't know if it's efficient for
vectorization, but without, it is quite efficient, as the locality is
very good.
gael.
So I don't think that our algorithms share common features with
band matrix algorithms. They seem closer to what can be done with
diagonal matrices.
Anyway, I've used the SparseMatrix code as a base for my
SkylineMatrix code. Currently it only supports very basic
operations. I'll implement the LU decomposition soon. I will also
try to use vectorization.
guillaume
These are just some initial thoughts and the discussion is
very open!
gael
On Mon, Oct 19, 2009 at 12:55 PM, guillaume saupin
<guillaume.saupin@xxxxxx <mailto:guillaume.saupin@xxxxxx>
<mailto:guillaume.saupin@xxxxxx
<mailto:guillaume.saupin@xxxxxx>>> wrote:
Hello,
We are planning to use your library in our projects, but we
need a
skyline matrix. Therefore I'd like to implement one for
eigen, but
I don't now where to start.
Is there a specific class that can be a good starting
point /
skeleton to start with ? The SparseMatrixBase might be a
good choice.
Should this SkylineMatrix inherit from SparseMatrixBase, or
be a
separate class ?
Thanks,
guillaume
# HG changeset patch
# User Gael Guennebaud <g.gael@xxxxxxx>
# Date 1257321890 -3600
# Node ID 3b71d57d8dda9623369038340afbebfdf613c104
# Parent 748e7700b89743270d6764609f8d577027bce1ae
an attempt to fix a compilation issue with MSVC
diff -r 748e7700b897 -r 3b71d57d8dda Eigen/src/Sparse/CholmodSupport.h
--- a/Eigen/src/Sparse/CholmodSupport.h Mon Nov 02 10:46:40 2009 +0100
+++ b/Eigen/src/Sparse/CholmodSupport.h Wed Nov 04 09:04:50 2009 +0100
@@ -126,6 +126,7 @@
typedef SparseLLT<MatrixType> Base;
typedef typename Base::Scalar Scalar;
typedef typename Base::RealScalar RealScalar;
+ typedef typename Base::CholMatrixType CholMatrixType;
using Base::MatrixLIsDirty;
using Base::SupernodalFactorIsDirty;
using Base::m_flags;
@@ -154,7 +155,7 @@
cholmod_finish(&m_cholmod);
}
- inline const typename Base::CholMatrixType& matrixL(void) const;
+ inline const CholMatrixType& matrixL() const;
template<typename Derived>
bool solveInPlace(MatrixBase<Derived> &b) const;
@@ -198,7 +199,7 @@
}
template<typename MatrixType>
-inline const typename SparseLLT<MatrixType>::CholMatrixType&
+inline const typename SparseLLT<MatrixType,Cholmod>::CholMatrixType&
SparseLLT<MatrixType,Cholmod>::matrixL() const
{
if (m_status & MatrixLIsDirty)
# HG changeset patch
# User kayhman@xxxxxxxxxxxxxxxxxxxx
# Date 1257344292 -3600
# Node ID 506ddadd3453ad6cc1dbce99d5eb5942e06ae07e
# Parent 3b71d57d8dda9623369038340afbebfdf613c104
Added basic SkylineMatrix.
diff -r 3b71d57d8dda -r 506ddadd3453 Eigen/Skyline
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/Skyline Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,38 @@
+#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 Sparse_Module Sparse module
+ *
+ * \nonstableyet
+ *
+ * See the \ref TutorialSparse "Sparse tutorial"
+ *
+ * \code
+ * #include <Eigen/QR>
+ * \endcode
+ */
+
+#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 -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Core/MatrixBase.h
--- a/Eigen/src/Core/MatrixBase.h Wed Nov 04 09:04:50 2009 +0100
+++ b/Eigen/src/Core/MatrixBase.h Wed Nov 04 15:18:12 2009 +0100
@@ -771,6 +771,11 @@
template<typename Derived1, typename Derived2>
Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
+// dense = skyline * dense
+ template<typename Derived1, typename Derived2>
+ Derived& lazyAssign(const SkylineProduct<Derived1,Derived2,SkylineTimeDenseProduct>& product);
+
+
////////// Householder module ///////////
void makeHouseholderInPlace(Scalar *tau, RealScalar *beta);
diff -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Core/util/Constants.h
--- a/Eigen/src/Core/util/Constants.h Wed Nov 04 09:04:50 2009 +0100
+++ b/Eigen/src/Core/util/Constants.h Wed Nov 04 15:18:12 2009 +0100
@@ -201,7 +201,7 @@
enum { ConditionalJumpCost = 5 };
enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
enum DirectionType { Vertical, Horizontal, BothDirections };
-enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct };
+enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct, SkylineTimeDenseProduct };
enum {
/** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment
diff -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Core/util/ForwardDeclarations.h
--- a/Eigen/src/Core/util/ForwardDeclarations.h Wed Nov 04 09:04:50 2009 +0100
+++ b/Eigen/src/Core/util/ForwardDeclarations.h Wed Nov 04 15:18:12 2009 +0100
@@ -142,5 +142,8 @@
// 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 -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Skyline/SkylineInplaceLU.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/src/Skyline/SkylineInplaceLU.h Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,372 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@xxxxxxx>
+//
+// 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
+
+//enum {
+// SvNoTrans = 0,
+// SvTranspose = 1,
+// SvAdjoint = 2
+//};
+
+/** \ingroup Skyline_Module
+ *
+ * \class SkylineInplaceLU
+ *
+ * \brief LU decomposition of a sparse matrix and associated features
+ *
+ * \param MatrixType the type of the matrix of which we are computing the LU factorization
+ *
+ * \sa class LU, class SkylineLLT
+ */
+template<typename MatrixType>
+class SkylineInplaceLU {
+protected:
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+ typedef SkylineMatrix<Scalar, LowerTriangular> LUMatrixType;
+
+ enum {
+ MatrixLUIsDirty = 0x10000
+ };
+
+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 -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Skyline/SkylineMatrix.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/src/Skyline/SkylineMatrix.h Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,895 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@xxxxxxx>
+//
+// 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, -=)
+ // FIXME: why are these operator already alvailable ???
+ // EIGEN_SKYLINE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SkylineMatrix, *=)
+ // EIGEN_SKYLINE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SkylineMatrix, /=)
+
+ // typedef MappedSkylineMatrix<Scalar,Flags> Map;
+ 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);
+ }
+
+ //--- low level purely coherent filling ---
+
+ inline Scalar & insertBack(int outer, int inner) {
+ // ei_assert(size_t(m_outerIndex[outer + 1]) == m_data.size() && "wrong sorted insertion");
+ // ei_assert((m_outerIndex[outer + 1] - m_outerIndex[outer] == 0 || m_data.index(m_data.size() - 1) < inner) && "wrong sorted insertion");
+ // int id = m_outerIndex[outer + 1];
+ // ++m_outerIndex[outer + 1];
+ // m_data.append(0, inner);
+ // return m_data.value(id);
+ //TODO
+ return 0;
+ }
+
+ inline void startVec(int outer) {
+ // ei_assert(m_outerIndex[outer] == int(m_data.size()) && "you must call startVec on each inner vec");
+ // ei_assert(m_outerIndex[outer + 1] == 0 && "you must call startVec on each inner vec");
+ // m_outerIndex[outer + 1] = m_outerIndex[outer];
+ }
+
+ //---
+
+ /** \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 -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Skyline/SkylineMatrixBase.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/src/Skyline/SkylineMatrixBase.h Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,371 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@xxxxxxx>
+//
+// 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 sparse matrices or sparse 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
+ };
+
+ /** \internal the return type of MatrixBase::conjugate() */
+ typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+ const SkylineCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Derived>,
+ const Derived&
+ >::ret ConjugateReturnType;
+ typedef SkylineCwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType;
+ /** \internal the return type of MatrixBase::imag() */
+ typedef SkylineCwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
+ /** \internal the return type of MatrixBase::adjoint() */
+ typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+ SkylineCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, SkylineNestByValue<Eigen::SkylineTranspose<Derived> > >,
+ SkylineTranspose<Derived>
+ >::ret AdjointReturnType;
+
+#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 true if either the number of rows or the number of columns is equal to 1.
+ * In other words, this function returns
+ * \code rows()==1 || cols()==1 \endcode
+ * \sa rows(), cols(), IsVectorAtCompileTime. */
+ inline bool isVector() const {
+ return rows() == 1 || cols() == 1;
+ }
+
+ /** \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;
+ }
+
+// const SkylineCwiseUnaryOp<ei_scalar_opposite_op<typename ei_traits<Derived>::Scalar>, Derived> operator-() const;
+//
+// template<typename OtherDerived>
+// const SkylineCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
+// operator+(const SkylineMatrixBase<OtherDerived> &other) const;
+//
+// template<typename OtherDerived>
+// const SkylineCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
+// operator-(const SkylineMatrixBase<OtherDerived> &other) const;
+
+// template<typename OtherDerived>
+// Derived & operator+=(const SkylineMatrixBase<OtherDerived>& other);
+// template<typename OtherDerived>
+// Derived & operator-=(const SkylineMatrixBase<OtherDerived>& other);
+
+
+// Derived & operator*=(const Scalar& other);
+// Derived & operator/=(const Scalar& other);
+//
+// const SkylineCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived>
+// operator*(const Scalar& scalar) const;
+// const SkylineCwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived>
+// operator/(const Scalar& scalar) const;
+
+// inline friend const SkylineCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived>
+// operator*(const Scalar& scalar, const SkylineMatrixBase& matrix) {
+// return matrix*scalar;
+// }
+
+
+// // skyline * skyline
+// template<typename OtherDerived>
+// const typename SkylineProductReturnType<Derived, OtherDerived>::Type
+// operator*(const SkylineMatrixBase<OtherDerived> &other) const;
+//
+// // skyline * diagonal
+// template<typename OtherDerived>
+// const SkylineDiagonalProduct<Derived, OtherDerived>
+// operator*(const DiagonalBase<OtherDerived> &other) const;
+//
+// // diagonal * skyline
+// template<typename OtherDerived> friend
+// const SkylineDiagonalProduct<OtherDerived, Derived>
+// operator*(const DiagonalBase<OtherDerived> &lhs, const SkylineMatrixBase& rhs) {
+// return SkylineDiagonalProduct<OtherDerived, Derived > (lhs.derived(), rhs.derived());
+// }
+
+// // dense * sparse (return a dense object)
+//
+// template<typename OtherDerived> friend
+// const typename SkylineProductReturnType<OtherDerived, Derived>::Type
+// operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs) {
+// return typename SkylineProductReturnType<OtherDerived, Derived>::Type(lhs.derived(), rhs);
+// }
+
+ template<typename OtherDerived>
+ const typename SkylineProductReturnType<Derived, OtherDerived>::Type
+ operator*(const MatrixBase<OtherDerived> &other) const;
+
+// template<typename OtherDerived>
+// Derived & operator*=(const SkylineMatrixBase<OtherDerived>& other);
+
+// template<int Mode>
+// inline const SkylineTriangular<Derived, Mode> triangular() const;
+
+// template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const;
+// template<typename OtherDerived> Scalar dot(const SkylineMatrixBase<OtherDerived>& other) const;
+// RealScalar squaredNorm() const;
+// RealScalar norm() const;
+
+ SkylineTranspose<Derived> transpose() {
+ return derived();
+ }
+
+ const SkylineTranspose<Derived> transpose() const {
+ return derived();
+ }
+ // void transposeInPlace();
+
+ const AdjointReturnType adjoint() const {
+ return transpose().nestByValue();
+ }
+
+ /** \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();
+ }
+//
+// template<typename OtherDerived>
+// bool isApprox(const SkylineMatrixBase<OtherDerived>& other,
+// RealScalar prec = precision<Scalar>()) const {
+// return toDense().isApprox(other.toDense(), prec);
+// }
+//
+// template<typename OtherDerived>
+// bool isApprox(const MatrixBase<OtherDerived>& other,
+// RealScalar prec = precision<Scalar>()) const {
+// return toDense().isApprox(other, prec);
+// }
+
+ template<typename NewType>
+ const SkylineCwiseUnaryOp<ei_scalar_cast_op<typename ei_traits<Derived>::Scalar, NewType>, Derived> cast() const;
+
+ /** \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());
+ }
+
+ // template<typename OtherDerived>
+ // void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other);
+
+ template<unsigned int Added>
+ const SkylineFlagged<Derived, Added, 0 > marked() const;
+ // const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const;
+
+ /** \returns number of elements to skip to pass from one row (resp. column) to another
+ * for a row-major (resp. column-major) matrix.
+ * Combined with coeffRef() and the \ref flags flags, it allows a direct access to the data
+ * of the underlying matrix.
+ */
+ // inline int stride(void) const { return derived().stride(); }
+
+ inline const SkylineNestByValue<Derived> nestByValue() const;
+
+
+ ConjugateReturnType conjugate() const;
+ const RealReturnType real() const;
+ const ImagReturnType imag() const;
+
+ template<typename CustomUnaryOp>
+ const SkylineCwiseUnaryOp<CustomUnaryOp, Derived> unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const;
+
+
+ Scalar sum() const;
+
+ const SkylineCwise<Derived> cwise() const;
+ SkylineCwise<Derived> cwise();
+
+
+protected:
+
+ bool m_isRValue;
+};
+
+#endif // EIGEN_SkylineMatrixBase_H
diff -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Skyline/SkylineProduct.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/src/Skyline/SkylineProduct.h Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,454 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@xxxxxxx>
+//
+// 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,
+
+ // LhsIsRowMajor = (LhsFlags & RowMajorBit)==RowMajorBit,
+ // RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit,
+
+ 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;
+};
+
+// perform a pseudo in-place skyline * skyline product assuming all matrices are col major
+
+template<typename Lhs, typename Rhs, typename ResultType>
+static void ei_skyline_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) {
+
+ //TODO
+
+
+ // typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
+ //
+ // // make sure to call innerSize/outerSize since we fake the storage order.
+ // int rows = lhs.innerSize();
+ // int cols = rhs.outerSize();
+ // //int size = lhs.outerSize();
+ // ei_assert(lhs.outerSize() == rhs.innerSize());
+ //
+ // // allocate a temporary buffer
+ // AmbiVector<Scalar> tempVector(rows);
+ //
+ // // estimate the number of non zero entries
+ // float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
+ // float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
+ // float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
+ //
+ // res.resize(rows, cols);
+ // res.reserve(int(ratioRes*rows*cols));
+ // for (int j=0; j<cols; ++j)
+ // {
+ // // let's do a more accurate determination of the nnz ratio for the current column j of res
+ // //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
+ // // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
+ // float ratioColRes = ratioRes;
+ // tempVector.init(ratioColRes);
+ // tempVector.setZero();
+ // for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
+ // {
+ // // FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
+ // tempVector.restart();
+ // Scalar x = rhsIt.value();
+ // for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
+ // {
+ // tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
+ // }
+ // }
+ // res.startVec(j);
+ // for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
+ // res.insertBack(j,it.index()) = it.value();
+ // }
+ // res.finalize();
+}
+
+
+
+//
+//template<typename Lhs, typename Rhs, typename ResultType>
+//struct ei_skyline_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
+//{
+// static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
+// {
+// // we need a col-major matrix to hold the result
+// typedef SkylineMatrix<typename ResultType::Scalar> SkylineTemporaryType;
+// SkylineTemporaryType _res(res.rows(), res.cols());
+// ei_skyline_product_impl<Lhs,Rhs,SkylineTemporaryType>(lhs, rhs, _res);
+// res = _res;
+// }
+//};
+//
+//template<typename Lhs, typename Rhs, typename ResultType>
+//struct ei_skyline_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
+//{
+// static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
+// {
+// // let's transpose the product to get a column x column product
+// typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
+// ei_skyline_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
+// res.swap(_res);
+// }
+//};
+//
+//template<typename Lhs, typename Rhs, typename ResultType>
+//struct ei_skyline_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
+//{
+// static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
+// {
+// // let's transpose the product to get a column x column product
+// typedef SkylineMatrix<typename ResultType::Scalar> SkylineTemporaryType;
+// SkylineTemporaryType _res(res.cols(), res.rows());
+// ei_skyline_product_impl<Rhs,Lhs,SkylineTemporaryType>(rhs, lhs, _res);
+// res = _res.transpose();
+// }
+//};
+
+//// skyline = skyline * skyline
+//template<typename Derived>
+//template<typename Lhs, typename Rhs>
+//inline Derived& SkylineMatrixBase<Derived>::operator=(const SkylineProduct<Lhs,Rhs,SkylineTimeSkylineProduct>& product)
+//{
+// ei_skyline_product_selector<
+// typename ei_cleantype<Lhs>::type,
+// typename ei_cleantype<Rhs>::type,
+// Derived>::run(product.lhs(),product.rhs(),derived());
+// return derived();
+//}
+
+// 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)
+ };
+
+ //ei_assert(_Lhs::IsRowMajor && "Not yet implemented for col major");
+
+ //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;
+ // if (_Lhs::IsRowMajor)
+ // ei_skyline_time_dense_product(product.lhs(), product.rhs(), derived());
+ // else
+ // ei_skyline_col_major_time_dense_product(product.lhs(), product.rhs(), derived());
+
+ ei_skyline_product_selector<typename ei_cleantype<Lhs>::type,
+ typename ei_cleantype<Rhs>::type,
+ Derived>::run(product.lhs(), product.rhs(), derived());
+
+ return derived();
+}
+
+// dense = dense * skyline
+
+//template<typename Derived>
+//template<typename Lhs, typename Rhs >
+//Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, DenseTimeSkylineProduct>& product) {
+// //TODO
+//
+// // typedef typename ei_cleantype<Rhs>::type _Rhs;
+// // typedef typename _Rhs::InnerIterator RhsInnerIterator;
+// // enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit };
+// // derived().setZero();
+// // for (int j=0; j<product.rhs().outerSize(); ++j)
+// // for (RhsInnerIterator i(product.rhs(),j); i; ++i)
+// // derived().col(RhsIsRowMajor ? i.index() : j) += i.value() * product.lhs().col(RhsIsRowMajor ? j : i.index());
+// // return derived();
+//}
+
+//// skyline * skyline
+//template<typename Derived>
+//template<typename OtherDerived>
+//EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived,OtherDerived>::Type
+//SkylineMatrixBase<Derived>::operator*(const SkylineMatrixBase<OtherDerived> &other) const
+//{
+// return typename SkylineProductReturnType<Derived,OtherDerived>::Type(derived(), other.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 -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Skyline/SkylineStorage.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/src/Skyline/SkylineStorage.h Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,266 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@xxxxxxx>
+//
+// 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 sparse set of values as a list of values and a list of indices.
+ *
+ */
+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 -r 3b71d57d8dda -r 506ddadd3453 Eigen/src/Skyline/SkylineUtil.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Eigen/src/Skyline/SkylineUtil.h Wed Nov 04 15:18:12 2009 +0100
@@ -0,0 +1,141 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@xxxxxxx>
+//
+// 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};
+//template<typename T, int Skyliness = ei_traits<T>::Flags&SkylineBit> class ei_eval;
+
+
+#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, };
+//typedef typename Eigen::ei_nested<Derived>::type Nested;
+//enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
+// ColsAtCompileTime = Eigen::ei_traits<Derived>::ColsAtCompileTime, \
+/ Flags = Eigen::ei_traits<Derived>::Flags, \
+// CoeffReadCost = Eigen::ei_traits<Derived>::CoeffReadCost, \
+// SizeAtCompileTime = Base::SizeAtCompileTime, \
+// IsVectorAtCompileTime = Base::IsVectorAtCompileTime };
+
+#define EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(Derived) \
+_EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::SkylineMatrixBase<Derived>)
+
+// solver flags
+//enum {
+// CompleteFactorization = 0x0000, // the default
+// IncompleteFactorization = 0x0001,
+// MemoryEfficient = 0x0002,
+//
+// // For LLT Cholesky:
+// SupernodalMultifrontal = 0x0010,
+// SupernodalLeftLooking = 0x0020,
+//
+// // Ordering methods:
+// NaturalOrdering = 0x0100, // the default
+// MinimumDegree_AT_PLUS_A = 0x0200,
+// MinimumDegree_ATA = 0x0300,
+// ColApproxMinimumDegree = 0x0400,
+// Metis = 0x0500,
+// Scotch = 0x0600,
+// Chaco = 0x0700,
+// OrderingMask = 0x0f00
+//};
+
+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 MatrixType> class SkylineNestByValue;
+template<typename MatrixType> class SkylineTranspose;
+template<typename MatrixType, int Size> class SkylineInnerVectorSet;
+template<typename Derived> class SkylineCwise;
+template<typename UnaryOp, typename MatrixType> class SkylineCwiseUnaryOp;
+template<typename BinaryOp, typename Lhs, typename Rhs> class SkylineCwiseBinaryOp;
+template<typename ExpressionType,
+ unsigned int Added, unsigned int Removed> class SkylineFlagged;
+template<typename ExpressionType, int Mode> class SkylineTriangular;
+template<typename Lhs, typename Rhs> class SkylineDiagonalProduct;
+
+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;
+
+//const int CoherentAccessPattern = 0x1;
+//const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern;
+//const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern;
+//const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern;
+
+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;
+};
+
+template<typename T> struct ei_must_nest_by_value<SkylineNestByValue<T> > { enum { ret = true }; };
+
+#endif // EIGEN_SKYLINEUTIL_H