SkylineMatrix.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SKYLINEMATRIX_H
00026 #define EIGEN_SKYLINEMATRIX_H
00027 
00028 #include "SkylineStorage.h"
00029 #include "SkylineMatrixBase.h"
00030 
00031 namespace Eigen { 
00032 
00048 namespace internal {
00049 template<typename _Scalar, int _Options>
00050 struct traits<SkylineMatrix<_Scalar, _Options> > {
00051     typedef _Scalar Scalar;
00052     typedef Sparse StorageKind;
00053 
00054     enum {
00055         RowsAtCompileTime = Dynamic,
00056         ColsAtCompileTime = Dynamic,
00057         MaxRowsAtCompileTime = Dynamic,
00058         MaxColsAtCompileTime = Dynamic,
00059         Flags = SkylineBit | _Options,
00060         CoeffReadCost = NumTraits<Scalar>::ReadCost,
00061     };
00062 };
00063 }
00064 
00065 template<typename _Scalar, int _Options>
00066 class SkylineMatrix
00067 : public SkylineMatrixBase<SkylineMatrix<_Scalar, _Options> > {
00068 public:
00069     EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(SkylineMatrix)
00070     EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, +=)
00071     EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, -=)
00072 
00073     using Base::IsRowMajor;
00074 
00075 protected:
00076 
00077     typedef SkylineMatrix<Scalar, (Flags&~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0) > TransposedSkylineMatrix;
00078 
00079     Index m_outerSize;
00080     Index m_innerSize;
00081 
00082 public:
00083     Index* m_colStartIndex;
00084     Index* m_rowStartIndex;
00085     SkylineStorage<Scalar> m_data;
00086 
00087 public:
00088 
00089     inline Index rows() const {
00090         return IsRowMajor ? m_outerSize : m_innerSize;
00091     }
00092 
00093     inline Index cols() const {
00094         return IsRowMajor ? m_innerSize : m_outerSize;
00095     }
00096 
00097     inline Index innerSize() const {
00098         return m_innerSize;
00099     }
00100 
00101     inline Index outerSize() const {
00102         return m_outerSize;
00103     }
00104 
00105     inline Index upperNonZeros() const {
00106         return m_data.upperSize();
00107     }
00108 
00109     inline Index lowerNonZeros() const {
00110         return m_data.lowerSize();
00111     }
00112 
00113     inline Index upperNonZeros(Index j) const {
00114         return m_colStartIndex[j + 1] - m_colStartIndex[j];
00115     }
00116 
00117     inline Index lowerNonZeros(Index j) const {
00118         return m_rowStartIndex[j + 1] - m_rowStartIndex[j];
00119     }
00120 
00121     inline const Scalar* _diagPtr() const {
00122         return &m_data.diag(0);
00123     }
00124 
00125     inline Scalar* _diagPtr() {
00126         return &m_data.diag(0);
00127     }
00128 
00129     inline const Scalar* _upperPtr() const {
00130         return &m_data.upper(0);
00131     }
00132 
00133     inline Scalar* _upperPtr() {
00134         return &m_data.upper(0);
00135     }
00136 
00137     inline const Scalar* _lowerPtr() const {
00138         return &m_data.lower(0);
00139     }
00140 
00141     inline Scalar* _lowerPtr() {
00142         return &m_data.lower(0);
00143     }
00144 
00145     inline const Index* _upperProfilePtr() const {
00146         return &m_data.upperProfile(0);
00147     }
00148 
00149     inline Index* _upperProfilePtr() {
00150         return &m_data.upperProfile(0);
00151     }
00152 
00153     inline const Index* _lowerProfilePtr() const {
00154         return &m_data.lowerProfile(0);
00155     }
00156 
00157     inline Index* _lowerProfilePtr() {
00158         return &m_data.lowerProfile(0);
00159     }
00160 
00161     inline Scalar coeff(Index row, Index col) const {
00162         const Index outer = IsRowMajor ? row : col;
00163         const Index inner = IsRowMajor ? col : row;
00164 
00165         eigen_assert(outer < outerSize());
00166         eigen_assert(inner < innerSize());
00167 
00168         if (outer == inner)
00169             return this->m_data.diag(outer);
00170 
00171         if (IsRowMajor) {
00172             if (inner > outer) //upper matrix
00173             {
00174                 const Index minOuterIndex = inner - m_data.upperProfile(inner);
00175                 if (outer >= minOuterIndex)
00176                     return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
00177                 else
00178                     return Scalar(0);
00179             }
00180             if (inner < outer) //lower matrix
00181             {
00182                 const Index minInnerIndex = outer - m_data.lowerProfile(outer);
00183                 if (inner >= minInnerIndex)
00184                     return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
00185                 else
00186                     return Scalar(0);
00187             }
00188             return m_data.upper(m_colStartIndex[inner] + outer - inner);
00189         } else {
00190             if (outer > inner) //upper matrix
00191             {
00192                 const Index maxOuterIndex = inner + m_data.upperProfile(inner);
00193                 if (outer <= maxOuterIndex)
00194                     return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
00195                 else
00196                     return Scalar(0);
00197             }
00198             if (outer < inner) //lower matrix
00199             {
00200                 const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
00201 
00202                 if (inner <= maxInnerIndex)
00203                     return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
00204                 else
00205                     return Scalar(0);
00206             }
00207         }
00208     }
00209 
00210     inline Scalar& coeffRef(Index row, Index col) {
00211         const Index outer = IsRowMajor ? row : col;
00212         const Index inner = IsRowMajor ? col : row;
00213 
00214         eigen_assert(outer < outerSize());
00215         eigen_assert(inner < innerSize());
00216 
00217         if (outer == inner)
00218             return this->m_data.diag(outer);
00219 
00220         if (IsRowMajor) {
00221             if (col > row) //upper matrix
00222             {
00223                 const Index minOuterIndex = inner - m_data.upperProfile(inner);
00224                 eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage");
00225                 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
00226             }
00227             if (col < row) //lower matrix
00228             {
00229                 const Index minInnerIndex = outer - m_data.lowerProfile(outer);
00230                 eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage");
00231                 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
00232             }
00233         } else {
00234             if (outer > inner) //upper matrix
00235             {
00236                 const Index maxOuterIndex = inner + m_data.upperProfile(inner);
00237                 eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage");
00238                 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
00239             }
00240             if (outer < inner) //lower matrix
00241             {
00242                 const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
00243                 eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage");
00244                 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
00245             }
00246         }
00247     }
00248 
00249     inline Scalar coeffDiag(Index idx) const {
00250         eigen_assert(idx < outerSize());
00251         eigen_assert(idx < innerSize());
00252         return this->m_data.diag(idx);
00253     }
00254 
00255     inline Scalar coeffLower(Index row, Index col) const {
00256         const Index outer = IsRowMajor ? row : col;
00257         const Index inner = IsRowMajor ? col : row;
00258 
00259         eigen_assert(outer < outerSize());
00260         eigen_assert(inner < innerSize());
00261         eigen_assert(inner != outer);
00262 
00263         if (IsRowMajor) {
00264             const Index minInnerIndex = outer - m_data.lowerProfile(outer);
00265             if (inner >= minInnerIndex)
00266                 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
00267             else
00268                 return Scalar(0);
00269 
00270         } else {
00271             const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
00272             if (inner <= maxInnerIndex)
00273                 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
00274             else
00275                 return Scalar(0);
00276         }
00277     }
00278 
00279     inline Scalar coeffUpper(Index row, Index col) const {
00280         const Index outer = IsRowMajor ? row : col;
00281         const Index inner = IsRowMajor ? col : row;
00282 
00283         eigen_assert(outer < outerSize());
00284         eigen_assert(inner < innerSize());
00285         eigen_assert(inner != outer);
00286 
00287         if (IsRowMajor) {
00288             const Index minOuterIndex = inner - m_data.upperProfile(inner);
00289             if (outer >= minOuterIndex)
00290                 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
00291             else
00292                 return Scalar(0);
00293         } else {
00294             const Index maxOuterIndex = inner + m_data.upperProfile(inner);
00295             if (outer <= maxOuterIndex)
00296                 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
00297             else
00298                 return Scalar(0);
00299         }
00300     }
00301 
00302     inline Scalar& coeffRefDiag(Index idx) {
00303         eigen_assert(idx < outerSize());
00304         eigen_assert(idx < innerSize());
00305         return this->m_data.diag(idx);
00306     }
00307 
00308     inline Scalar& coeffRefLower(Index row, Index col) {
00309         const Index outer = IsRowMajor ? row : col;
00310         const Index inner = IsRowMajor ? col : row;
00311 
00312         eigen_assert(outer < outerSize());
00313         eigen_assert(inner < innerSize());
00314         eigen_assert(inner != outer);
00315 
00316         if (IsRowMajor) {
00317             const Index minInnerIndex = outer - m_data.lowerProfile(outer);
00318             eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage");
00319             return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
00320         } else {
00321             const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
00322             eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage");
00323             return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
00324         }
00325     }
00326 
00327     inline bool coeffExistLower(Index row, Index col) {
00328         const Index outer = IsRowMajor ? row : col;
00329         const Index inner = IsRowMajor ? col : row;
00330 
00331         eigen_assert(outer < outerSize());
00332         eigen_assert(inner < innerSize());
00333         eigen_assert(inner != outer);
00334 
00335         if (IsRowMajor) {
00336             const Index minInnerIndex = outer - m_data.lowerProfile(outer);
00337             return inner >= minInnerIndex;
00338         } else {
00339             const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
00340             return inner <= maxInnerIndex;
00341         }
00342     }
00343 
00344     inline Scalar& coeffRefUpper(Index row, Index col) {
00345         const Index outer = IsRowMajor ? row : col;
00346         const Index inner = IsRowMajor ? col : row;
00347 
00348         eigen_assert(outer < outerSize());
00349         eigen_assert(inner < innerSize());
00350         eigen_assert(inner != outer);
00351 
00352         if (IsRowMajor) {
00353             const Index minOuterIndex = inner - m_data.upperProfile(inner);
00354             eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage");
00355             return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
00356         } else {
00357             const Index maxOuterIndex = inner + m_data.upperProfile(inner);
00358             eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage");
00359             return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
00360         }
00361     }
00362 
00363     inline bool coeffExistUpper(Index row, Index col) {
00364         const Index outer = IsRowMajor ? row : col;
00365         const Index inner = IsRowMajor ? col : row;
00366 
00367         eigen_assert(outer < outerSize());
00368         eigen_assert(inner < innerSize());
00369         eigen_assert(inner != outer);
00370 
00371         if (IsRowMajor) {
00372             const Index minOuterIndex = inner - m_data.upperProfile(inner);
00373             return outer >= minOuterIndex;
00374         } else {
00375             const Index maxOuterIndex = inner + m_data.upperProfile(inner);
00376             return outer <= maxOuterIndex;
00377         }
00378     }
00379 
00380 
00381 protected:
00382 
00383 public:
00384     class InnerUpperIterator;
00385     class InnerLowerIterator;
00386 
00387     class OuterUpperIterator;
00388     class OuterLowerIterator;
00389 
00391     inline void setZero() {
00392         m_data.clear();
00393         memset(m_colStartIndex, 0, (m_outerSize + 1) * sizeof (Index));
00394         memset(m_rowStartIndex, 0, (m_outerSize + 1) * sizeof (Index));
00395     }
00396 
00398     inline Index nonZeros() const {
00399         return m_data.diagSize() + m_data.upperSize() + m_data.lowerSize();
00400     }
00401 
00403     inline void reserve(Index reserveSize, Index reserveUpperSize, Index reserveLowerSize) {
00404         m_data.reserve(reserveSize, reserveUpperSize, reserveLowerSize);
00405     }
00406 
00415     EIGEN_DONT_INLINE Scalar & insert(Index row, Index col) {
00416         const Index outer = IsRowMajor ? row : col;
00417         const Index inner = IsRowMajor ? col : row;
00418 
00419         eigen_assert(outer < outerSize());
00420         eigen_assert(inner < innerSize());
00421 
00422         if (outer == inner)
00423             return m_data.diag(col);
00424 
00425         if (IsRowMajor) {
00426             if (outer < inner) //upper matrix
00427             {
00428                 Index minOuterIndex = 0;
00429                 minOuterIndex = inner - m_data.upperProfile(inner);
00430 
00431                 if (outer < minOuterIndex) //The value does not yet exist
00432                 {
00433                     const Index previousProfile = m_data.upperProfile(inner);
00434 
00435                     m_data.upperProfile(inner) = inner - outer;
00436 
00437 
00438                     const Index bandIncrement = m_data.upperProfile(inner) - previousProfile;
00439                     //shift data stored after this new one
00440                     const Index stop = m_colStartIndex[cols()];
00441                     const Index start = m_colStartIndex[inner];
00442 
00443 
00444                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
00445                         m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx);
00446                     }
00447 
00448                     for (Index innerIdx = cols(); innerIdx > inner; innerIdx--) {
00449                         m_colStartIndex[innerIdx] += bandIncrement;
00450                     }
00451 
00452                     //zeros new data
00453                     memset(this->_upperPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar));
00454 
00455                     return m_data.upper(m_colStartIndex[inner]);
00456                 } else {
00457                     return m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
00458                 }
00459             }
00460 
00461             if (outer > inner) //lower matrix
00462             {
00463                 const Index minInnerIndex = outer - m_data.lowerProfile(outer);
00464                 if (inner < minInnerIndex) //The value does not yet exist
00465                 {
00466                     const Index previousProfile = m_data.lowerProfile(outer);
00467                     m_data.lowerProfile(outer) = outer - inner;
00468 
00469                     const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile;
00470                     //shift data stored after this new one
00471                     const Index stop = m_rowStartIndex[rows()];
00472                     const Index start = m_rowStartIndex[outer];
00473 
00474 
00475                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
00476                         m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx);
00477                     }
00478 
00479                     for (Index innerIdx = rows(); innerIdx > outer; innerIdx--) {
00480                         m_rowStartIndex[innerIdx] += bandIncrement;
00481                     }
00482 
00483                     //zeros new data
00484                     memset(this->_lowerPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar));
00485                     return m_data.lower(m_rowStartIndex[outer]);
00486                 } else {
00487                     return m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
00488                 }
00489             }
00490         } else {
00491             if (outer > inner) //upper matrix
00492             {
00493                 const Index maxOuterIndex = inner + m_data.upperProfile(inner);
00494                 if (outer > maxOuterIndex) //The value does not yet exist
00495                 {
00496                     const Index previousProfile = m_data.upperProfile(inner);
00497                     m_data.upperProfile(inner) = outer - inner;
00498 
00499                     const Index bandIncrement = m_data.upperProfile(inner) - previousProfile;
00500                     //shift data stored after this new one
00501                     const Index stop = m_rowStartIndex[rows()];
00502                     const Index start = m_rowStartIndex[inner + 1];
00503 
00504                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
00505                         m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx);
00506                     }
00507 
00508                     for (Index innerIdx = inner + 1; innerIdx < outerSize() + 1; innerIdx++) {
00509                         m_rowStartIndex[innerIdx] += bandIncrement;
00510                     }
00511                     memset(this->_upperPtr() + m_rowStartIndex[inner] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar));
00512                     return m_data.upper(m_rowStartIndex[inner] + m_data.upperProfile(inner));
00513                 } else {
00514                     return m_data.upper(m_rowStartIndex[inner] + (outer - inner));
00515                 }
00516             }
00517 
00518             if (outer < inner) //lower matrix
00519             {
00520                 const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
00521                 if (inner > maxInnerIndex) //The value does not yet exist
00522                 {
00523                     const Index previousProfile = m_data.lowerProfile(outer);
00524                     m_data.lowerProfile(outer) = inner - outer;
00525 
00526                     const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile;
00527                     //shift data stored after this new one
00528                     const Index stop = m_colStartIndex[cols()];
00529                     const Index start = m_colStartIndex[outer + 1];
00530 
00531                     for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
00532                         m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx);
00533                     }
00534 
00535                     for (Index innerIdx = outer + 1; innerIdx < outerSize() + 1; innerIdx++) {
00536                         m_colStartIndex[innerIdx] += bandIncrement;
00537                     }
00538                     memset(this->_lowerPtr() + m_colStartIndex[outer] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar));
00539                     return m_data.lower(m_colStartIndex[outer] + m_data.lowerProfile(outer));
00540                 } else {
00541                     return m_data.lower(m_colStartIndex[outer] + (inner - outer));
00542                 }
00543             }
00544         }
00545     }
00546 
00549     inline void finalize() {
00550         if (IsRowMajor) {
00551             if (rows() > cols())
00552                 m_data.resize(cols(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1);
00553             else
00554                 m_data.resize(rows(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1);
00555 
00556             //            eigen_assert(rows() == cols() && "memory reorganisatrion only works with suare matrix");
00557             //
00558             //            Scalar* newArray = new Scalar[m_colStartIndex[cols()] + 1 + m_rowStartIndex[rows()] + 1];
00559             //            Index dataIdx = 0;
00560             //            for (Index row = 0; row < rows(); row++) {
00561             //
00562             //                const Index nbLowerElts = m_rowStartIndex[row + 1] - m_rowStartIndex[row];
00563             //                //                std::cout << "nbLowerElts" << nbLowerElts << std::endl;
00564             //                memcpy(newArray + dataIdx, m_data.m_lower + m_rowStartIndex[row], nbLowerElts * sizeof (Scalar));
00565             //                m_rowStartIndex[row] = dataIdx;
00566             //                dataIdx += nbLowerElts;
00567             //
00568             //                const Index nbUpperElts = m_colStartIndex[row + 1] - m_colStartIndex[row];
00569             //                memcpy(newArray + dataIdx, m_data.m_upper + m_colStartIndex[row], nbUpperElts * sizeof (Scalar));
00570             //                m_colStartIndex[row] = dataIdx;
00571             //                dataIdx += nbUpperElts;
00572             //
00573             //
00574             //            }
00575             //            //todo : don't access m_data profile directly : add an accessor from SkylineMatrix
00576             //            m_rowStartIndex[rows()] = m_rowStartIndex[rows()-1] + m_data.lowerProfile(rows()-1);
00577             //            m_colStartIndex[cols()] = m_colStartIndex[cols()-1] + m_data.upperProfile(cols()-1);
00578             //
00579             //            delete[] m_data.m_lower;
00580             //            delete[] m_data.m_upper;
00581             //
00582             //            m_data.m_lower = newArray;
00583             //            m_data.m_upper = newArray;
00584         } else {
00585             if (rows() > cols())
00586                 m_data.resize(cols(), rows(), cols(), m_rowStartIndex[cols()] + 1, m_colStartIndex[cols()] + 1);
00587             else
00588                 m_data.resize(rows(), rows(), cols(), m_rowStartIndex[rows()] + 1, m_colStartIndex[rows()] + 1);
00589         }
00590     }
00591 
00592     inline void squeeze() {
00593         finalize();
00594         m_data.squeeze();
00595     }
00596 
00597     void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar > ()) {
00598         //TODO
00599     }
00600 
00604     void resize(size_t rows, size_t cols) {
00605         const Index diagSize = rows > cols ? cols : rows;
00606         m_innerSize = IsRowMajor ? cols : rows;
00607 
00608         eigen_assert(rows == cols && "Skyline matrix must be square matrix");
00609 
00610         if (diagSize % 2) { // diagSize is odd
00611             const Index k = (diagSize - 1) / 2;
00612 
00613             m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols,
00614                     2 * k * k + k + 1,
00615                     2 * k * k + k + 1);
00616 
00617         } else // diagSize is even
00618         {
00619             const Index k = diagSize / 2;
00620             m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols,
00621                     2 * k * k - k + 1,
00622                     2 * k * k - k + 1);
00623         }
00624 
00625         if (m_colStartIndex && m_rowStartIndex) {
00626             delete[] m_colStartIndex;
00627             delete[] m_rowStartIndex;
00628         }
00629         m_colStartIndex = new Index [cols + 1];
00630         m_rowStartIndex = new Index [rows + 1];
00631         m_outerSize = diagSize;
00632 
00633         m_data.reset();
00634         m_data.clear();
00635 
00636         m_outerSize = diagSize;
00637         memset(m_colStartIndex, 0, (cols + 1) * sizeof (Index));
00638         memset(m_rowStartIndex, 0, (rows + 1) * sizeof (Index));
00639     }
00640 
00641     void resizeNonZeros(Index size) {
00642         m_data.resize(size);
00643     }
00644 
00645     inline SkylineMatrix()
00646     : m_outerSize(-1), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
00647         resize(0, 0);
00648     }
00649 
00650     inline SkylineMatrix(size_t rows, size_t cols)
00651     : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
00652         resize(rows, cols);
00653     }
00654 
00655     template<typename OtherDerived>
00656     inline SkylineMatrix(const SkylineMatrixBase<OtherDerived>& other)
00657     : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
00658         *this = other.derived();
00659     }
00660 
00661     inline SkylineMatrix(const SkylineMatrix & other)
00662     : Base(), m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
00663         *this = other.derived();
00664     }
00665 
00666     inline void swap(SkylineMatrix & other) {
00667         //EIGEN_DBG_SKYLINE(std::cout << "SkylineMatrix:: swap\n");
00668         std::swap(m_colStartIndex, other.m_colStartIndex);
00669         std::swap(m_rowStartIndex, other.m_rowStartIndex);
00670         std::swap(m_innerSize, other.m_innerSize);
00671         std::swap(m_outerSize, other.m_outerSize);
00672         m_data.swap(other.m_data);
00673     }
00674 
00675     inline SkylineMatrix & operator=(const SkylineMatrix & other) {
00676         std::cout << "SkylineMatrix& operator=(const SkylineMatrix& other)\n";
00677         if (other.isRValue()) {
00678             swap(other.const_cast_derived());
00679         } else {
00680             resize(other.rows(), other.cols());
00681             memcpy(m_colStartIndex, other.m_colStartIndex, (m_outerSize + 1) * sizeof (Index));
00682             memcpy(m_rowStartIndex, other.m_rowStartIndex, (m_outerSize + 1) * sizeof (Index));
00683             m_data = other.m_data;
00684         }
00685         return *this;
00686     }
00687 
00688     template<typename OtherDerived>
00689             inline SkylineMatrix & operator=(const SkylineMatrixBase<OtherDerived>& other) {
00690         const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
00691         if (needToTranspose) {
00692             //         TODO
00693             //            return *this;
00694         } else {
00695             // there is no special optimization
00696             return SkylineMatrixBase<SkylineMatrix>::operator=(other.derived());
00697         }
00698     }
00699 
00700     friend std::ostream & operator <<(std::ostream & s, const SkylineMatrix & m) {
00701 
00702         EIGEN_DBG_SKYLINE(
00703         std::cout << "upper elements : " << std::endl;
00704         for (Index i = 0; i < m.m_data.upperSize(); i++)
00705             std::cout << m.m_data.upper(i) << "\t";
00706         std::cout << std::endl;
00707         std::cout << "upper profile : " << std::endl;
00708         for (Index i = 0; i < m.m_data.upperProfileSize(); i++)
00709             std::cout << m.m_data.upperProfile(i) << "\t";
00710         std::cout << std::endl;
00711         std::cout << "lower startIdx : " << std::endl;
00712         for (Index i = 0; i < m.m_data.upperProfileSize(); i++)
00713             std::cout << (IsRowMajor ? m.m_colStartIndex[i] : m.m_rowStartIndex[i]) << "\t";
00714         std::cout << std::endl;
00715 
00716 
00717         std::cout << "lower elements : " << std::endl;
00718         for (Index i = 0; i < m.m_data.lowerSize(); i++)
00719             std::cout << m.m_data.lower(i) << "\t";
00720         std::cout << std::endl;
00721         std::cout << "lower profile : " << std::endl;
00722         for (Index i = 0; i < m.m_data.lowerProfileSize(); i++)
00723             std::cout << m.m_data.lowerProfile(i) << "\t";
00724         std::cout << std::endl;
00725         std::cout << "lower startIdx : " << std::endl;
00726         for (Index i = 0; i < m.m_data.lowerProfileSize(); i++)
00727             std::cout << (IsRowMajor ? m.m_rowStartIndex[i] : m.m_colStartIndex[i]) << "\t";
00728         std::cout << std::endl;
00729         );
00730         for (Index rowIdx = 0; rowIdx < m.rows(); rowIdx++) {
00731             for (Index colIdx = 0; colIdx < m.cols(); colIdx++) {
00732                 s << m.coeff(rowIdx, colIdx) << "\t";
00733             }
00734             s << std::endl;
00735         }
00736         return s;
00737     }
00738 
00740     inline ~SkylineMatrix() {
00741         delete[] m_colStartIndex;
00742         delete[] m_rowStartIndex;
00743     }
00744 
00746     Scalar sum() const;
00747 };
00748 
00749 template<typename Scalar, int _Options>
00750 class SkylineMatrix<Scalar, _Options>::InnerUpperIterator {
00751 public:
00752 
00753     InnerUpperIterator(const SkylineMatrix& mat, Index outer)
00754     : m_matrix(mat), m_outer(outer),
00755     m_id(_Options == RowMajor ? mat.m_colStartIndex[outer] : mat.m_rowStartIndex[outer] + 1),
00756     m_start(m_id),
00757     m_end(_Options == RowMajor ? mat.m_colStartIndex[outer + 1] : mat.m_rowStartIndex[outer + 1] + 1) {
00758     }
00759 
00760     inline InnerUpperIterator & operator++() {
00761         m_id++;
00762         return *this;
00763     }
00764 
00765     inline InnerUpperIterator & operator+=(Index shift) {
00766         m_id += shift;
00767         return *this;
00768     }
00769 
00770     inline Scalar value() const {
00771         return m_matrix.m_data.upper(m_id);
00772     }
00773 
00774     inline Scalar* valuePtr() {
00775         return const_cast<Scalar*> (&(m_matrix.m_data.upper(m_id)));
00776     }
00777 
00778     inline Scalar& valueRef() {
00779         return const_cast<Scalar&> (m_matrix.m_data.upper(m_id));
00780     }
00781 
00782     inline Index index() const {
00783         return IsRowMajor ? m_outer - m_matrix.m_data.upperProfile(m_outer) + (m_id - m_start) :
00784                 m_outer + (m_id - m_start) + 1;
00785     }
00786 
00787     inline Index row() const {
00788         return IsRowMajor ? index() : m_outer;
00789     }
00790 
00791     inline Index col() const {
00792         return IsRowMajor ? m_outer : index();
00793     }
00794 
00795     inline size_t size() const {
00796         return m_matrix.m_data.upperProfile(m_outer);
00797     }
00798 
00799     inline operator bool() const {
00800         return (m_id < m_end) && (m_id >= m_start);
00801     }
00802 
00803 protected:
00804     const SkylineMatrix& m_matrix;
00805     const Index m_outer;
00806     Index m_id;
00807     const Index m_start;
00808     const Index m_end;
00809 };
00810 
00811 template<typename Scalar, int _Options>
00812 class SkylineMatrix<Scalar, _Options>::InnerLowerIterator {
00813 public:
00814 
00815     InnerLowerIterator(const SkylineMatrix& mat, Index outer)
00816     : m_matrix(mat),
00817     m_outer(outer),
00818     m_id(_Options == RowMajor ? mat.m_rowStartIndex[outer] : mat.m_colStartIndex[outer] + 1),
00819     m_start(m_id),
00820     m_end(_Options == RowMajor ? mat.m_rowStartIndex[outer + 1] : mat.m_colStartIndex[outer + 1] + 1) {
00821     }
00822 
00823     inline InnerLowerIterator & operator++() {
00824         m_id++;
00825         return *this;
00826     }
00827 
00828     inline InnerLowerIterator & operator+=(Index shift) {
00829         m_id += shift;
00830         return *this;
00831     }
00832 
00833     inline Scalar value() const {
00834         return m_matrix.m_data.lower(m_id);
00835     }
00836 
00837     inline Scalar* valuePtr() {
00838         return const_cast<Scalar*> (&(m_matrix.m_data.lower(m_id)));
00839     }
00840 
00841     inline Scalar& valueRef() {
00842         return const_cast<Scalar&> (m_matrix.m_data.lower(m_id));
00843     }
00844 
00845     inline Index index() const {
00846         return IsRowMajor ? m_outer - m_matrix.m_data.lowerProfile(m_outer) + (m_id - m_start) :
00847                 m_outer + (m_id - m_start) + 1;
00848         ;
00849     }
00850 
00851     inline Index row() const {
00852         return IsRowMajor ? m_outer : index();
00853     }
00854 
00855     inline Index col() const {
00856         return IsRowMajor ? index() : m_outer;
00857     }
00858 
00859     inline size_t size() const {
00860         return m_matrix.m_data.lowerProfile(m_outer);
00861     }
00862 
00863     inline operator bool() const {
00864         return (m_id < m_end) && (m_id >= m_start);
00865     }
00866 
00867 protected:
00868     const SkylineMatrix& m_matrix;
00869     const Index m_outer;
00870     Index m_id;
00871     const Index m_start;
00872     const Index m_end;
00873 };
00874 
00875 } // end namespace Eigen
00876 
00877 #endif // EIGEN_SkylineMatrix_H