src/Skyline/SkylineMatrix.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SKYLINEMATRIX_H
0011 #define EIGEN_SKYLINEMATRIX_H
0012
0013 #include "SkylineStorage.h"
0014 #include "SkylineMatrixBase.h"
0015
0016 namespace Eigen {
0017
0018 /** \ingroup Skyline_Module
0019  *
0020  * \class SkylineMatrix
0021  *
0022  * \brief The main skyline matrix class
0023  *
0024  * This class implements a skyline matrix using the very uncommon storage
0025  * scheme.
0026  *
0027  * \param _Scalar the scalar type, i.e. the type of the coefficients
0028  * \param _Options Union of bit flags controlling the storage scheme. Currently the only possibility
0029  *                 is RowMajor. The default is 0 which means column-major.
0030  *
0031  *
0032  */
0033 namespace internal {
0034 template<typename _Scalar, int _Options>
0035 struct traits<SkylineMatrix<_Scalar, _Options> > {
0036     typedef _Scalar Scalar;
0037     typedef Sparse StorageKind;
0038
0039     enum {
0040         RowsAtCompileTime = Dynamic,
0041         ColsAtCompileTime = Dynamic,
0042         MaxRowsAtCompileTime = Dynamic,
0043         MaxColsAtCompileTime = Dynamic,
0044         Flags = SkylineBit | _Options,
0045         CoeffReadCost = NumTraits<Scalar>::ReadCost,
0046     };
0047 };
0048 }
0049
0050 template<typename _Scalar, int _Options>
0051 class SkylineMatrix
0052 : public SkylineMatrixBase<SkylineMatrix<_Scalar, _Options> > {
0053 public:
0054     EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(SkylineMatrix)
0055     EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, +=)
0056     EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, -=)
0057
0058     using Base::IsRowMajor;
0059
0060 protected:
0061
0062     typedef SkylineMatrix<Scalar, (Flags&~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0) > TransposedSkylineMatrix;
0063
0064     Index m_outerSize;
0065     Index m_innerSize;
0066
0067 public:
0068     Index* m_colStartIndex;
0069     Index* m_rowStartIndex;
0070     SkylineStorage<Scalar> m_data;
0071
0072 public:
0073
0074     inline Index rows() const {
0075         return IsRowMajor ? m_outerSize : m_innerSize;
0076     }
0077
0078     inline Index cols() const {
0079         return IsRowMajor ? m_innerSize : m_outerSize;
0080     }
0081
0082     inline Index innerSize() const {
0083         return m_innerSize;
0084     }
0085
0086     inline Index outerSize() const {
0087         return m_outerSize;
0088     }
0089
0090     inline Index upperNonZeros() const {
0091         return m_data.upperSize();
0092     }
0093
0094     inline Index lowerNonZeros() const {
0095         return m_data.lowerSize();
0096     }
0097
0098     inline Index upperNonZeros(Index j) const {
0099         return m_colStartIndex[j + 1] - m_colStartIndex[j];
0100     }
0101
0102     inline Index lowerNonZeros(Index j) const {
0103         return m_rowStartIndex[j + 1] - m_rowStartIndex[j];
0104     }
0105
0106     inline const Scalar* _diagPtr() const {
0107         return &m_data.diag(0);
0108     }
0109
0110     inline Scalar* _diagPtr() {
0111         return &m_data.diag(0);
0112     }
0113
0114     inline const Scalar* _upperPtr() const {
0115         return &m_data.upper(0);
0116     }
0117
0118     inline Scalar* _upperPtr() {
0119         return &m_data.upper(0);
0120     }
0121
0122     inline const Scalar* _lowerPtr() const {
0123         return &m_data.lower(0);
0124     }
0125
0126     inline Scalar* _lowerPtr() {
0127         return &m_data.lower(0);
0128     }
0129
0130     inline const Index* _upperProfilePtr() const {
0131         return &m_data.upperProfile(0);
0132     }
0133
0134     inline Index* _upperProfilePtr() {
0135         return &m_data.upperProfile(0);
0136     }
0137
0138     inline const Index* _lowerProfilePtr() const {
0139         return &m_data.lowerProfile(0);
0140     }
0141
0142     inline Index* _lowerProfilePtr() {
0143         return &m_data.lowerProfile(0);
0144     }
0145
0146     inline Scalar coeff(Index row, Index col) const {
0147         const Index outer = IsRowMajor ? row : col;
0148         const Index inner = IsRowMajor ? col : row;
0149
0150         eigen_assert(outer < outerSize());
0151         eigen_assert(inner < innerSize());
0152
0153         if (outer == inner)
0154             return this->m_data.diag(outer);
0155
0156         if (IsRowMajor) {
0157             if (inner > outer) //upper matrix
0158             {
0159                 const Index minOuterIndex = inner - m_data.upperProfile(inner);
0160                 if (outer >= minOuterIndex)
0161                     return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
0162                 else
0163                     return Scalar(0);
0164             }
0165             if (inner < outer) //lower matrix
0166             {
0167                 const Index minInnerIndex = outer - m_data.lowerProfile(outer);
0168                 if (inner >= minInnerIndex)
0169                     return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
0170                 else
0171                     return Scalar(0);
0172             }
0173             return m_data.upper(m_colStartIndex[inner] + outer - inner);
0174         } else {
0175             if (outer > inner) //upper matrix
0176             {
0177                 const Index maxOuterIndex = inner + m_data.upperProfile(inner);
0178                 if (outer <= maxOuterIndex)
0179                     return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
0180                 else
0181                     return Scalar(0);
0182             }
0183             if (outer < inner) //lower matrix
0184             {
0185                 const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
0186
0187                 if (inner <= maxInnerIndex)
0188                     return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
0189                 else
0190                     return Scalar(0);
0191             }
0192         }
0193     }
0194
0195     inline Scalar& coeffRef(Index row, Index col) {
0196         const Index outer = IsRowMajor ? row : col;
0197         const Index inner = IsRowMajor ? col : row;
0198
0199         eigen_assert(outer < outerSize());
0200         eigen_assert(inner < innerSize());
0201
0202         if (outer == inner)
0203             return this->m_data.diag(outer);
0204
0205         if (IsRowMajor) {
0206             if (col > row) //upper matrix
0207             {
0208                 const Index minOuterIndex = inner - m_data.upperProfile(inner);
0209                 eigen_assert(outer >= minOuterIndex && "You tried to access a coeff that does not exist in the storage");
0210                 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
0211             }
0212             if (col < row) //lower matrix
0213             {
0214                 const Index minInnerIndex = outer - m_data.lowerProfile(outer);
0215                 eigen_assert(inner >= minInnerIndex && "You tried to access a coeff that does not exist in the storage");
0216                 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
0217             }
0218         } else {
0219             if (outer > inner) //upper matrix
0220             {
0221                 const Index maxOuterIndex = inner + m_data.upperProfile(inner);
0222                 eigen_assert(outer <= maxOuterIndex && "You tried to access a coeff that does not exist in the storage");
0223                 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
0224             }
0225             if (outer < inner) //lower matrix
0226             {
0227                 const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
0228                 eigen_assert(inner <= maxInnerIndex && "You tried to access a coeff that does not exist in the storage");
0229                 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
0230             }
0231         }
0232     }
0233
0234     inline Scalar coeffDiag(Index idx) const {
0235         eigen_assert(idx < outerSize());
0236         eigen_assert(idx < innerSize());
0237         return this->m_data.diag(idx);
0238     }
0239
0240     inline Scalar coeffLower(Index row, Index col) const {
0241         const Index outer = IsRowMajor ? row : col;
0242         const Index inner = IsRowMajor ? col : row;
0243
0244         eigen_assert(outer < outerSize());
0245         eigen_assert(inner < innerSize());
0246         eigen_assert(inner != outer);
0247
0248         if (IsRowMajor) {
0249             const Index minInnerIndex = outer - m_data.lowerProfile(outer);
0250             if (inner >= minInnerIndex)
0251                 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
0252             else
0253                 return Scalar(0);
0254
0255         } else {
0256             const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
0257             if (inner <= maxInnerIndex)
0258                 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
0259             else
0260                 return Scalar(0);
0261         }
0262     }
0263
0264     inline Scalar coeffUpper(Index row, Index col) const {
0265         const Index outer = IsRowMajor ? row : col;
0266         const Index inner = IsRowMajor ? col : row;
0267
0268         eigen_assert(outer < outerSize());
0269         eigen_assert(inner < innerSize());
0270         eigen_assert(inner != outer);
0271
0272         if (IsRowMajor) {
0273             const Index minOuterIndex = inner - m_data.upperProfile(inner);
0274             if (outer >= minOuterIndex)
0275                 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
0276             else
0277                 return Scalar(0);
0278         } else {
0279             const Index maxOuterIndex = inner + m_data.upperProfile(inner);
0280             if (outer <= maxOuterIndex)
0281                 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
0282             else
0283                 return Scalar(0);
0284         }
0285     }
0286
0287     inline Scalar& coeffRefDiag(Index idx) {
0288         eigen_assert(idx < outerSize());
0289         eigen_assert(idx < innerSize());
0290         return this->m_data.diag(idx);
0291     }
0292
0293     inline Scalar& coeffRefLower(Index row, Index col) {
0294         const Index outer = IsRowMajor ? row : col;
0295         const Index inner = IsRowMajor ? col : row;
0296
0297         eigen_assert(outer < outerSize());
0298         eigen_assert(inner < innerSize());
0299         eigen_assert(inner != outer);
0300
0301         if (IsRowMajor) {
0302             const Index minInnerIndex = outer - m_data.lowerProfile(outer);
0303             eigen_assert(inner >= minInnerIndex && "You tried to access a coeff that does not exist in the storage");
0304             return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
0305         } else {
0306             const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
0307             eigen_assert(inner <= maxInnerIndex && "You tried to access a coeff that does not exist in the storage");
0308             return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
0309         }
0310     }
0311
0312     inline bool coeffExistLower(Index row, Index col) {
0313         const Index outer = IsRowMajor ? row : col;
0314         const Index inner = IsRowMajor ? col : row;
0315
0316         eigen_assert(outer < outerSize());
0317         eigen_assert(inner < innerSize());
0318         eigen_assert(inner != outer);
0319
0320         if (IsRowMajor) {
0321             const Index minInnerIndex = outer - m_data.lowerProfile(outer);
0322             return inner >= minInnerIndex;
0323         } else {
0324             const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
0325             return inner <= maxInnerIndex;
0326         }
0327     }
0328
0329     inline Scalar& coeffRefUpper(Index row, Index col) {
0330         const Index outer = IsRowMajor ? row : col;
0331         const Index inner = IsRowMajor ? col : row;
0332
0333         eigen_assert(outer < outerSize());
0334         eigen_assert(inner < innerSize());
0335         eigen_assert(inner != outer);
0336
0337         if (IsRowMajor) {
0338             const Index minOuterIndex = inner - m_data.upperProfile(inner);
0339             eigen_assert(outer >= minOuterIndex && "You tried to access a coeff that does not exist in the storage");
0340             return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
0341         } else {
0342             const Index maxOuterIndex = inner + m_data.upperProfile(inner);
0343             eigen_assert(outer <= maxOuterIndex && "You tried to access a coeff that does not exist in the storage");
0344             return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
0345         }
0346     }
0347
0348     inline bool coeffExistUpper(Index row, Index col) {
0349         const Index outer = IsRowMajor ? row : col;
0350         const Index inner = IsRowMajor ? col : row;
0351
0352         eigen_assert(outer < outerSize());
0353         eigen_assert(inner < innerSize());
0354         eigen_assert(inner != outer);
0355
0356         if (IsRowMajor) {
0357             const Index minOuterIndex = inner - m_data.upperProfile(inner);
0358             return outer >= minOuterIndex;
0359         } else {
0360             const Index maxOuterIndex = inner + m_data.upperProfile(inner);
0361             return outer <= maxOuterIndex;
0362         }
0363     }
0364
0365
0366 protected:
0367
0368 public:
0369     class InnerUpperIterator;
0370     class InnerLowerIterator;
0371
0372     class OuterUpperIterator;
0373     class OuterLowerIterator;
0374
0375     /** Removes all non zeros */
0376     inline void setZero() {
0377         m_data.clear();
0378         memset(m_colStartIndex, 0, (m_outerSize + 1) * sizeof (Index));
0379         memset(m_rowStartIndex, 0, (m_outerSize + 1) * sizeof (Index));
0380     }
0381
0382     /** \returns the number of non zero coefficients */
0383     inline Index nonZeros() const {
0384         return m_data.diagSize() + m_data.upperSize() + m_data.lowerSize();
0385     }
0386
0387     /** Preallocates \a reserveSize non zeros */
0388     inline void reserve(Index reserveSize, Index reserveUpperSize, Index reserveLowerSize) {
0389         m_data.reserve(reserveSize, reserveUpperSize, reserveLowerSize);
0390     }
0391
0392     /** \returns a reference to a novel non zero coefficient with coordinates \a row x \a col.
0393
0394      *
0395      * \warning This function can be extremely slow if the non zero coefficients
0396      * are not inserted in a coherent order.
0397      *
0398      * After an insertion session, you should call the finalize() function.
0399      */
0400     EIGEN_DONT_INLINE Scalar & insert(Index row, Index col) {
0401         const Index outer = IsRowMajor ? row : col;
0402         const Index inner = IsRowMajor ? col : row;
0403
0404         eigen_assert(outer < outerSize());
0405         eigen_assert(inner < innerSize());
0406
0407         if (outer == inner)
0408             return m_data.diag(col);
0409
0410         if (IsRowMajor) {
0411             if (outer < inner) //upper matrix
0412             {
0413                 Index minOuterIndex = 0;
0414                 minOuterIndex = inner - m_data.upperProfile(inner);
0415
0416                 if (outer < minOuterIndex) //The value does not yet exist
0417                 {
0418                     const Index previousProfile = m_data.upperProfile(inner);
0419
0420                     m_data.upperProfile(inner) = inner - outer;
0421
0422
0423                     const Index bandIncrement = m_data.upperProfile(inner) - previousProfile;
0424                     //shift data stored after this new one
0425                     const Index stop = m_colStartIndex[cols()];
0426                     const Index start = m_colStartIndex[inner];
0427
0428
0429                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
0430                         m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx);
0431                     }
0432
0433                     for (Index innerIdx = cols(); innerIdx > inner; innerIdx--) {
0434                         m_colStartIndex[innerIdx] += bandIncrement;
0435                     }
0436
0437                     //zeros new data
0438                     memset(this->_upperPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar));
0439
0440                     return m_data.upper(m_colStartIndex[inner]);
0441                 } else {
0442                     return m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
0443                 }
0444             }
0445
0446             if (outer > inner) //lower matrix
0447             {
0448                 const Index minInnerIndex = outer - m_data.lowerProfile(outer);
0449                 if (inner < minInnerIndex) //The value does not yet exist
0450                 {
0451                     const Index previousProfile = m_data.lowerProfile(outer);
0452                     m_data.lowerProfile(outer) = outer - inner;
0453
0454                     const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile;
0455                     //shift data stored after this new one
0456                     const Index stop = m_rowStartIndex[rows()];
0457                     const Index start = m_rowStartIndex[outer];
0458
0459
0460                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
0461                         m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx);
0462                     }
0463
0464                     for (Index innerIdx = rows(); innerIdx > outer; innerIdx--) {
0465                         m_rowStartIndex[innerIdx] += bandIncrement;
0466                     }
0467
0468                     //zeros new data
0469                     memset(this->_lowerPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar));
0470                     return m_data.lower(m_rowStartIndex[outer]);
0471                 } else {
0472                     return m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
0473                 }
0474             }
0475         } else {
0476             if (outer > inner) //upper matrix
0477             {
0478                 const Index maxOuterIndex = inner + m_data.upperProfile(inner);
0479                 if (outer > maxOuterIndex) //The value does not yet exist
0480                 {
0481                     const Index previousProfile = m_data.upperProfile(inner);
0482                     m_data.upperProfile(inner) = outer - inner;
0483
0484                     const Index bandIncrement = m_data.upperProfile(inner) - previousProfile;
0485                     //shift data stored after this new one
0486                     const Index stop = m_rowStartIndex[rows()];
0487                     const Index start = m_rowStartIndex[inner + 1];
0488
0489                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
0490                         m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx);
0491                     }
0492
0493                     for (Index innerIdx = inner + 1; innerIdx < outerSize() + 1; innerIdx++) {
0494                         m_rowStartIndex[innerIdx] += bandIncrement;
0495                     }
0496                     memset(this->_upperPtr() + m_rowStartIndex[inner] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar));
0497                     return m_data.upper(m_rowStartIndex[inner] + m_data.upperProfile(inner));
0498                 } else {
0499                     return m_data.upper(m_rowStartIndex[inner] + (outer - inner));
0500                 }
0501             }
0502
0503             if (outer < inner) //lower matrix
0504             {
0505                 const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
0506                 if (inner > maxInnerIndex) //The value does not yet exist
0507                 {
0508                     const Index previousProfile = m_data.lowerProfile(outer);
0509                     m_data.lowerProfile(outer) = inner - outer;
0510
0511                     const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile;
0512                     //shift data stored after this new one
0513                     const Index stop = m_colStartIndex[cols()];
0514                     const Index start = m_colStartIndex[outer + 1];
0515
0516                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
0517                         m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx);
0518                     }
0519
0520                     for (Index innerIdx = outer + 1; innerIdx < outerSize() + 1; innerIdx++) {
0521                         m_colStartIndex[innerIdx] += bandIncrement;
0522                     }
0523                     memset(this->_lowerPtr() + m_colStartIndex[outer] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar));
0524                     return m_data.lower(m_colStartIndex[outer] + m_data.lowerProfile(outer));
0525                 } else {
0526                     return m_data.lower(m_colStartIndex[outer] + (inner - outer));
0527                 }
0528             }
0529         }
0530     }
0531
0532     /** Must be called after inserting a set of non zero entries.
0533      */
0534     inline void finalize() {
0535         if (IsRowMajor) {
0536             if (rows() > cols())
0537                 m_data.resize(cols(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1);
0538             else
0539                 m_data.resize(rows(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1);
0540
0541             //            eigen_assert(rows() == cols() && "memory reorganisatrion only works with suare matrix");
0542             //
0543             //            Scalar* newArray = new Scalar[m_colStartIndex[cols()] + 1 + m_rowStartIndex[rows()] + 1];
0544             //            Index dataIdx = 0;
0545             //            for (Index row = 0; row < rows(); row++) {
0546             //
0547             //                const Index nbLowerElts = m_rowStartIndex[row + 1] - m_rowStartIndex[row];
0548             //                //                std::cout << "nbLowerElts" << nbLowerElts << std::endl;
0549             //                memcpy(newArray + dataIdx, m_data.m_lower + m_rowStartIndex[row], nbLowerElts * sizeof (Scalar));
0550             //                m_rowStartIndex[row] = dataIdx;
0551             //                dataIdx += nbLowerElts;
0552             //
0553             //                const Index nbUpperElts = m_colStartIndex[row + 1] - m_colStartIndex[row];
0554             //                memcpy(newArray + dataIdx, m_data.m_upper + m_colStartIndex[row], nbUpperElts * sizeof (Scalar));
0555             //                m_colStartIndex[row] = dataIdx;
0556             //                dataIdx += nbUpperElts;
0557             //
0558             //
0559             //            }
0560             //            //todo : don't access m_data profile directly : add an accessor from SkylineMatrix
0561             //            m_rowStartIndex[rows()] = m_rowStartIndex[rows()-1] + m_data.lowerProfile(rows()-1);
0562             //            m_colStartIndex[cols()] = m_colStartIndex[cols()-1] + m_data.upperProfile(cols()-1);
0563             //
0564             //            delete[] m_data.m_lower;
0565             //            delete[] m_data.m_upper;
0566             //
0567             //            m_data.m_lower = newArray;
0568             //            m_data.m_upper = newArray;
0569         } else {
0570             if (rows() > cols())
0571                 m_data.resize(cols(), rows(), cols(), m_rowStartIndex[cols()] + 1, m_colStartIndex[cols()] + 1);
0572             else
0573                 m_data.resize(rows(), rows(), cols(), m_rowStartIndex[rows()] + 1, m_colStartIndex[rows()] + 1);
0574         }
0575     }
0576
0577     inline void squeeze() {
0578         finalize();
0579         m_data.squeeze();
0580     }
0581
0582     void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar > ()) {
0583         //TODO
0584     }
0585
0586     /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero
0587      * \sa resizeNonZeros(Index), reserve(), setZero()
0588      */
0589     void resize(size_t rows, size_t cols) {
0590         const Index diagSize = rows > cols ? cols : rows;
0591         m_innerSize = IsRowMajor ? cols : rows;
0592
0593         eigen_assert(rows == cols && "Skyline matrix must be square matrix");
0594
0595         if (diagSize % 2) { // diagSize is odd
0596             const Index k = (diagSize - 1) / 2;
0597
0598             m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols,
0599                     2 * k * k + k + 1,
0600                     2 * k * k + k + 1);
0601
0602         } else // diagSize is even
0603         {
0604             const Index k = diagSize / 2;
0605             m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols,
0606                     2 * k * k - k + 1,
0607                     2 * k * k - k + 1);
0608         }
0609
0610         if (m_colStartIndex && m_rowStartIndex) {
0611             delete[] m_colStartIndex;
0612             delete[] m_rowStartIndex;
0613         }
0614         m_colStartIndex = new Index [cols + 1];
0615         m_rowStartIndex = new Index [rows + 1];
0616         m_outerSize = diagSize;
0617
0618         m_data.reset();
0619         m_data.clear();
0620
0621         m_outerSize = diagSize;
0622         memset(m_colStartIndex, 0, (cols + 1) * sizeof (Index));
0623         memset(m_rowStartIndex, 0, (rows + 1) * sizeof (Index));
0624     }
0625
0626     void resizeNonZeros(Index size) {
0627         m_data.resize(size);
0628     }
0629
0630     inline SkylineMatrix()
0631     : m_outerSize(-1), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
0632         resize(0, 0);
0633     }
0634
0635     inline SkylineMatrix(size_t rows, size_t cols)
0636     : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
0637         resize(rows, cols);
0638     }
0639
0640     template<typename OtherDerived>
0641     inline SkylineMatrix(const SkylineMatrixBase<OtherDerived>& other)
0642     : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
0643         *this = other.derived();
0644     }
0645
0646     inline SkylineMatrix(const SkylineMatrix & other)
0647     : Base(), m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
0648         *this = other.derived();
0649     }
0650
0651     inline void swap(SkylineMatrix & other) {
0652         //EIGEN_DBG_SKYLINE(std::cout << "SkylineMatrix:: swap\n");
0653         std::swap(m_colStartIndex, other.m_colStartIndex);
0654         std::swap(m_rowStartIndex, other.m_rowStartIndex);
0655         std::swap(m_innerSize, other.m_innerSize);
0656         std::swap(m_outerSize, other.m_outerSize);
0657         m_data.swap(other.m_data);
0658     }
0659
0660     inline SkylineMatrix & operator=(const SkylineMatrix & other) {
0661         std::cout << "SkylineMatrix& operator=(const SkylineMatrix& other)\n";
0662         if (other.isRValue()) {
0663             swap(other.const_cast_derived());
0664         } else {
0665             resize(other.rows(), other.cols());
0666             memcpy(m_colStartIndex, other.m_colStartIndex, (m_outerSize + 1) * sizeof (Index));
0667             memcpy(m_rowStartIndex, other.m_rowStartIndex, (m_outerSize + 1) * sizeof (Index));
0668             m_data = other.m_data;
0669         }
0670         return *this;
0671     }
0672
0673     template<typename OtherDerived>
0674             inline SkylineMatrix & operator=(const SkylineMatrixBase<OtherDerived>& other) {
0675         const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
0676         if (needToTranspose) {
0677             //         TODO
0678             //            return *this;
0679         } else {
0680             // there is no special optimization
0681             return SkylineMatrixBase<SkylineMatrix>::operator=(other.derived());
0682         }
0683     }
0684
0685     friend std::ostream & operator <<(std::ostream & s, const SkylineMatrix & m) {
0686
0687         EIGEN_DBG_SKYLINE(
0688         std::cout << "upper elements : " << std::endl;
0689         for (Index i = 0; i < m.m_data.upperSize(); i++)
0690             std::cout << m.m_data.upper(i) << "\t";
0691         std::cout << std::endl;
0692         std::cout << "upper profile : " << std::endl;
0693         for (Index i = 0; i < m.m_data.upperProfileSize(); i++)
0694             std::cout << m.m_data.upperProfile(i) << "\t";
0695         std::cout << std::endl;
0696         std::cout << "lower startIdx : " << std::endl;
0697         for (Index i = 0; i < m.m_data.upperProfileSize(); i++)
0698             std::cout << (IsRowMajor ? m.m_colStartIndex[i] : m.m_rowStartIndex[i]) << "\t";
0699         std::cout << std::endl;
0700
0701
0702         std::cout << "lower elements : " << std::endl;
0703         for (Index i = 0; i < m.m_data.lowerSize(); i++)
0704             std::cout << m.m_data.lower(i) << "\t";
0705         std::cout << std::endl;
0706         std::cout << "lower profile : " << std::endl;
0707         for (Index i = 0; i < m.m_data.lowerProfileSize(); i++)
0708             std::cout << m.m_data.lowerProfile(i) << "\t";
0709         std::cout << std::endl;
0710         std::cout << "lower startIdx : " << std::endl;
0711         for (Index i = 0; i < m.m_data.lowerProfileSize(); i++)
0712             std::cout << (IsRowMajor ? m.m_rowStartIndex[i] : m.m_colStartIndex[i]) << "\t";
0713         std::cout << std::endl;
0714         );
0715         for (Index rowIdx = 0; rowIdx < m.rows(); rowIdx++) {
0716             for (Index colIdx = 0; colIdx < m.cols(); colIdx++) {
0717                 s << m.coeff(rowIdx, colIdx) << "\t";
0718             }
0719             s << std::endl;
0720         }
0721         return s;
0722     }
0723
0724     /** Destructor */
0725     inline ~SkylineMatrix() {
0726         delete[] m_colStartIndex;
0727         delete[] m_rowStartIndex;
0728     }
0729
0730     /** Overloaded for performance */
0731     Scalar sum() const;
0732 };
0733
0734 template<typename Scalar, int _Options>
0735 class SkylineMatrix<Scalar, _Options>::InnerUpperIterator {
0736 public:
0737
0738     InnerUpperIterator(const SkylineMatrix& mat, Index outer)
0739     : m_matrix(mat), m_outer(outer),
0740     m_id(_Options == RowMajor ? mat.m_colStartIndex[outer] : mat.m_rowStartIndex[outer] + 1),
0741     m_start(m_id),
0742     m_end(_Options == RowMajor ? mat.m_colStartIndex[outer + 1] : mat.m_rowStartIndex[outer + 1] + 1) {
0743     }
0744
0745     inline InnerUpperIterator & operator++() {
0746         m_id++;
0747         return *this;
0748     }
0749
0750     inline InnerUpperIterator & operator+=(Index shift) {
0751         m_id += shift;
0752         return *this;
0753     }
0754
0755     inline Scalar value() const {
0756         return m_matrix.m_data.upper(m_id);
0757     }
0758
0759     inline Scalar* valuePtr() {
0760         return const_cast<Scalar*> (&(m_matrix.m_data.upper(m_id)));
0761     }
0762
0763     inline Scalar& valueRef() {
0764         return const_cast<Scalar&> (m_matrix.m_data.upper(m_id));
0765     }
0766
0767     inline Index index() const {
0768         return IsRowMajor ? m_outer - m_matrix.m_data.upperProfile(m_outer) + (m_id - m_start) :
0769                 m_outer + (m_id - m_start) + 1;
0770     }
0771
0772     inline Index row() const {
0773         return IsRowMajor ? index() : m_outer;
0774     }
0775
0776     inline Index col() const {
0777         return IsRowMajor ? m_outer : index();
0778     }
0779
0780     inline size_t size() const {
0781         return m_matrix.m_data.upperProfile(m_outer);
0782     }
0783
0784     inline operator bool() const {
0785         return (m_id < m_end) && (m_id >= m_start);
0786     }
0787
0788 protected:
0789     const SkylineMatrix& m_matrix;
0790     const Index m_outer;
0791     Index m_id;
0792     const Index m_start;
0793     const Index m_end;
0794 };
0795
0796 template<typename Scalar, int _Options>
0797 class SkylineMatrix<Scalar, _Options>::InnerLowerIterator {
0798 public:
0799
0800     InnerLowerIterator(const SkylineMatrix& mat, Index outer)
0801     : m_matrix(mat),
0802     m_outer(outer),
0803     m_id(_Options == RowMajor ? mat.m_rowStartIndex[outer] : mat.m_colStartIndex[outer] + 1),
0804     m_start(m_id),
0805     m_end(_Options == RowMajor ? mat.m_rowStartIndex[outer + 1] : mat.m_colStartIndex[outer + 1] + 1) {
0806     }
0807
0808     inline InnerLowerIterator & operator++() {
0809         m_id++;
0810         return *this;
0811     }
0812
0813     inline InnerLowerIterator & operator+=(Index shift) {
0814         m_id += shift;
0815         return *this;
0816     }
0817
0818     inline Scalar value() const {
0819         return m_matrix.m_data.lower(m_id);
0820     }
0821
0822     inline Scalar* valuePtr() {
0823         return const_cast<Scalar*> (&(m_matrix.m_data.lower(m_id)));
0824     }
0825
0826     inline Scalar& valueRef() {
0827         return const_cast<Scalar&> (m_matrix.m_data.lower(m_id));
0828     }
0829
0830     inline Index index() const {
0831         return IsRowMajor ? m_outer - m_matrix.m_data.lowerProfile(m_outer) + (m_id - m_start) :
0832                 m_outer + (m_id - m_start) + 1;
0833         ;
0834     }
0835
0836     inline Index row() const {
0837         return IsRowMajor ? m_outer : index();
0838     }
0839
0840     inline Index col() const {
0841         return IsRowMajor ? index() : m_outer;
0842     }
0843
0844     inline size_t size() const {
0845         return m_matrix.m_data.lowerProfile(m_outer);
0846     }
0847
0848     inline operator bool() const {
0849         return (m_id < m_end) && (m_id >= m_start);
0850     }
0851
0852 protected:
0853     const SkylineMatrix& m_matrix;
0854     const Index m_outer;
0855     Index m_id;
0856     const Index m_start;
0857     const Index m_end;
0858 };
0859
0860 } // end namespace Eigen
0861
0862 #endif // EIGEN_SKYLINEMATRIX_H