1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_SKYLINEMATRIX_H 11 #define EIGEN_SKYLINEMATRIX_H 12 13 #include "SkylineStorage.h" 14 #include "SkylineMatrixBase.h" 15 16 namespace Eigen { 17 18 /** \ingroup Skyline_Module 19 * 20 * \class SkylineMatrix 21 * 22 * \brief The main skyline matrix class 23 * 24 * This class implements a skyline matrix using the very uncommon storage 25 * scheme. 26 * 27 * \param _Scalar the scalar type, i.e. the type of the coefficients 28 * \param _Options Union of bit flags controlling the storage scheme. Currently the only possibility 29 * is RowMajor. The default is 0 which means column-major. 30 * 31 * 32 */ 33 namespace internal { 34 template<typename _Scalar, int _Options> 35 struct traits<SkylineMatrix<_Scalar, _Options> > { 36 typedef _Scalar Scalar; 37 typedef Sparse StorageKind; 38 39 enum { 40 RowsAtCompileTime = Dynamic, 41 ColsAtCompileTime = Dynamic, 42 MaxRowsAtCompileTime = Dynamic, 43 MaxColsAtCompileTime = Dynamic, 44 Flags = SkylineBit | _Options, 45 CoeffReadCost = NumTraits<Scalar>::ReadCost, 46 }; 47 }; 48 } 49 50 template<typename _Scalar, int _Options> 51 class SkylineMatrix 52 : public SkylineMatrixBase<SkylineMatrix<_Scalar, _Options> > { 53 public: 54 EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(SkylineMatrix) 55 EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, +=) 56 EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, -=) 57 58 using Base::IsRowMajor; 59 60 protected: 61 62 typedef SkylineMatrix<Scalar, (Flags&~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0) > TransposedSkylineMatrix; 63 64 Index m_outerSize; 65 Index m_innerSize; 66 67 public: 68 Index* m_colStartIndex; 69 Index* m_rowStartIndex; 70 SkylineStorage<Scalar> m_data; 71 72 public: 73 74 inline Index rows() const { 75 return IsRowMajor ? m_outerSize : m_innerSize; 76 } 77 78 inline Index cols() const { 79 return IsRowMajor ? m_innerSize : m_outerSize; 80 } 81 82 inline Index innerSize() const { 83 return m_innerSize; 84 } 85 86 inline Index outerSize() const { 87 return m_outerSize; 88 } 89 90 inline Index upperNonZeros() const { 91 return m_data.upperSize(); 92 } 93 94 inline Index lowerNonZeros() const { 95 return m_data.lowerSize(); 96 } 97 98 inline Index upperNonZeros(Index j) const { 99 return m_colStartIndex[j + 1] - m_colStartIndex[j]; 100 } 101 102 inline Index lowerNonZeros(Index j) const { 103 return m_rowStartIndex[j + 1] - m_rowStartIndex[j]; 104 } 105 106 inline const Scalar* _diagPtr() const { 107 return &m_data.diag(0); 108 } 109 110 inline Scalar* _diagPtr() { 111 return &m_data.diag(0); 112 } 113 114 inline const Scalar* _upperPtr() const { 115 return &m_data.upper(0); 116 } 117 118 inline Scalar* _upperPtr() { 119 return &m_data.upper(0); 120 } 121 122 inline const Scalar* _lowerPtr() const { 123 return &m_data.lower(0); 124 } 125 126 inline Scalar* _lowerPtr() { 127 return &m_data.lower(0); 128 } 129 130 inline const Index* _upperProfilePtr() const { 131 return &m_data.upperProfile(0); 132 } 133 134 inline Index* _upperProfilePtr() { 135 return &m_data.upperProfile(0); 136 } 137 138 inline const Index* _lowerProfilePtr() const { 139 return &m_data.lowerProfile(0); 140 } 141 142 inline Index* _lowerProfilePtr() { 143 return &m_data.lowerProfile(0); 144 } 145 146 inline Scalar coeff(Index row, Index col) const { 147 const Index outer = IsRowMajor ? row : col; 148 const Index inner = IsRowMajor ? col : row; 149 150 eigen_assert(outer < outerSize()); 151 eigen_assert(inner < innerSize()); 152 153 if (outer == inner) 154 return this->m_data.diag(outer); 155 156 if (IsRowMajor) { 157 if (inner > outer) //upper matrix 158 { 159 const Index minOuterIndex = inner - m_data.upperProfile(inner); 160 if (outer >= minOuterIndex) 161 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); 162 else 163 return Scalar(0); 164 } 165 if (inner < outer) //lower matrix 166 { 167 const Index minInnerIndex = outer - m_data.lowerProfile(outer); 168 if (inner >= minInnerIndex) 169 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); 170 else 171 return Scalar(0); 172 } 173 return m_data.upper(m_colStartIndex[inner] + outer - inner); 174 } else { 175 if (outer > inner) //upper matrix 176 { 177 const Index maxOuterIndex = inner + m_data.upperProfile(inner); 178 if (outer <= maxOuterIndex) 179 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); 180 else 181 return Scalar(0); 182 } 183 if (outer < inner) //lower matrix 184 { 185 const Index maxInnerIndex = outer + m_data.lowerProfile(outer); 186 187 if (inner <= maxInnerIndex) 188 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); 189 else 190 return Scalar(0); 191 } 192 } 193 } 194 195 inline Scalar& coeffRef(Index row, Index col) { 196 const Index outer = IsRowMajor ? row : col; 197 const Index inner = IsRowMajor ? col : row; 198 199 eigen_assert(outer < outerSize()); 200 eigen_assert(inner < innerSize()); 201 202 if (outer == inner) 203 return this->m_data.diag(outer); 204 205 if (IsRowMajor) { 206 if (col > row) //upper matrix 207 { 208 const Index minOuterIndex = inner - m_data.upperProfile(inner); 209 eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage"); 210 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); 211 } 212 if (col < row) //lower matrix 213 { 214 const Index minInnerIndex = outer - m_data.lowerProfile(outer); 215 eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage"); 216 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); 217 } 218 } else { 219 if (outer > inner) //upper matrix 220 { 221 const Index maxOuterIndex = inner + m_data.upperProfile(inner); 222 eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage"); 223 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); 224 } 225 if (outer < inner) //lower matrix 226 { 227 const Index maxInnerIndex = outer + m_data.lowerProfile(outer); 228 eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage"); 229 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); 230 } 231 } 232 } 233 234 inline Scalar coeffDiag(Index idx) const { 235 eigen_assert(idx < outerSize()); 236 eigen_assert(idx < innerSize()); 237 return this->m_data.diag(idx); 238 } 239 240 inline Scalar coeffLower(Index row, Index col) const { 241 const Index outer = IsRowMajor ? row : col; 242 const Index inner = IsRowMajor ? col : row; 243 244 eigen_assert(outer < outerSize()); 245 eigen_assert(inner < innerSize()); 246 eigen_assert(inner != outer); 247 248 if (IsRowMajor) { 249 const Index minInnerIndex = outer - m_data.lowerProfile(outer); 250 if (inner >= minInnerIndex) 251 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); 252 else 253 return Scalar(0); 254 255 } else { 256 const Index maxInnerIndex = outer + m_data.lowerProfile(outer); 257 if (inner <= maxInnerIndex) 258 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); 259 else 260 return Scalar(0); 261 } 262 } 263 264 inline Scalar coeffUpper(Index row, Index col) const { 265 const Index outer = IsRowMajor ? row : col; 266 const Index inner = IsRowMajor ? col : row; 267 268 eigen_assert(outer < outerSize()); 269 eigen_assert(inner < innerSize()); 270 eigen_assert(inner != outer); 271 272 if (IsRowMajor) { 273 const Index minOuterIndex = inner - m_data.upperProfile(inner); 274 if (outer >= minOuterIndex) 275 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); 276 else 277 return Scalar(0); 278 } else { 279 const Index maxOuterIndex = inner + m_data.upperProfile(inner); 280 if (outer <= maxOuterIndex) 281 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); 282 else 283 return Scalar(0); 284 } 285 } 286 287 inline Scalar& coeffRefDiag(Index idx) { 288 eigen_assert(idx < outerSize()); 289 eigen_assert(idx < innerSize()); 290 return this->m_data.diag(idx); 291 } 292 293 inline Scalar& coeffRefLower(Index row, Index col) { 294 const Index outer = IsRowMajor ? row : col; 295 const Index inner = IsRowMajor ? col : row; 296 297 eigen_assert(outer < outerSize()); 298 eigen_assert(inner < innerSize()); 299 eigen_assert(inner != outer); 300 301 if (IsRowMajor) { 302 const Index minInnerIndex = outer - m_data.lowerProfile(outer); 303 eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage"); 304 return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); 305 } else { 306 const Index maxInnerIndex = outer + m_data.lowerProfile(outer); 307 eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage"); 308 return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer)); 309 } 310 } 311 312 inline bool coeffExistLower(Index row, Index col) { 313 const Index outer = IsRowMajor ? row : col; 314 const Index inner = IsRowMajor ? col : row; 315 316 eigen_assert(outer < outerSize()); 317 eigen_assert(inner < innerSize()); 318 eigen_assert(inner != outer); 319 320 if (IsRowMajor) { 321 const Index minInnerIndex = outer - m_data.lowerProfile(outer); 322 return inner >= minInnerIndex; 323 } else { 324 const Index maxInnerIndex = outer + m_data.lowerProfile(outer); 325 return inner <= maxInnerIndex; 326 } 327 } 328 329 inline Scalar& coeffRefUpper(Index row, Index col) { 330 const Index outer = IsRowMajor ? row : col; 331 const Index inner = IsRowMajor ? col : row; 332 333 eigen_assert(outer < outerSize()); 334 eigen_assert(inner < innerSize()); 335 eigen_assert(inner != outer); 336 337 if (IsRowMajor) { 338 const Index minOuterIndex = inner - m_data.upperProfile(inner); 339 eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage"); 340 return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); 341 } else { 342 const Index maxOuterIndex = inner + m_data.upperProfile(inner); 343 eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage"); 344 return this->m_data.upper(m_colStartIndex[inner] + (outer - inner)); 345 } 346 } 347 348 inline bool coeffExistUpper(Index row, Index col) { 349 const Index outer = IsRowMajor ? row : col; 350 const Index inner = IsRowMajor ? col : row; 351 352 eigen_assert(outer < outerSize()); 353 eigen_assert(inner < innerSize()); 354 eigen_assert(inner != outer); 355 356 if (IsRowMajor) { 357 const Index minOuterIndex = inner - m_data.upperProfile(inner); 358 return outer >= minOuterIndex; 359 } else { 360 const Index maxOuterIndex = inner + m_data.upperProfile(inner); 361 return outer <= maxOuterIndex; 362 } 363 } 364 365 366 protected: 367 368 public: 369 class InnerUpperIterator; 370 class InnerLowerIterator; 371 372 class OuterUpperIterator; 373 class OuterLowerIterator; 374 375 /** Removes all non zeros */ 376 inline void setZero() { 377 m_data.clear(); 378 memset(m_colStartIndex, 0, (m_outerSize + 1) * sizeof (Index)); 379 memset(m_rowStartIndex, 0, (m_outerSize + 1) * sizeof (Index)); 380 } 381 382 /** \returns the number of non zero coefficients */ 383 inline Index nonZeros() const { 384 return m_data.diagSize() + m_data.upperSize() + m_data.lowerSize(); 385 } 386 387 /** Preallocates \a reserveSize non zeros */ 388 inline void reserve(Index reserveSize, Index reserveUpperSize, Index reserveLowerSize) { 389 m_data.reserve(reserveSize, reserveUpperSize, reserveLowerSize); 390 } 391 392 /** \returns a reference to a novel non zero coefficient with coordinates \a row x \a col. 393 394 * 395 * \warning This function can be extremely slow if the non zero coefficients 396 * are not inserted in a coherent order. 397 * 398 * After an insertion session, you should call the finalize() function. 399 */ 400 EIGEN_DONT_INLINE Scalar & insert(Index row, Index col) { 401 const Index outer = IsRowMajor ? row : col; 402 const Index inner = IsRowMajor ? col : row; 403 404 eigen_assert(outer < outerSize()); 405 eigen_assert(inner < innerSize()); 406 407 if (outer == inner) 408 return m_data.diag(col); 409 410 if (IsRowMajor) { 411 if (outer < inner) //upper matrix 412 { 413 Index minOuterIndex = 0; 414 minOuterIndex = inner - m_data.upperProfile(inner); 415 416 if (outer < minOuterIndex) //The value does not yet exist 417 { 418 const Index previousProfile = m_data.upperProfile(inner); 419 420 m_data.upperProfile(inner) = inner - outer; 421 422 423 const Index bandIncrement = m_data.upperProfile(inner) - previousProfile; 424 //shift data stored after this new one 425 const Index stop = m_colStartIndex[cols()]; 426 const Index start = m_colStartIndex[inner]; 427 428 429 for (Index innerIdx = stop; innerIdx >= start; innerIdx--) { 430 m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx); 431 } 432 433 for (Index innerIdx = cols(); innerIdx > inner; innerIdx--) { 434 m_colStartIndex[innerIdx] += bandIncrement; 435 } 436 437 //zeros new data 438 memset(this->_upperPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar)); 439 440 return m_data.upper(m_colStartIndex[inner]); 441 } else { 442 return m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner))); 443 } 444 } 445 446 if (outer > inner) //lower matrix 447 { 448 const Index minInnerIndex = outer - m_data.lowerProfile(outer); 449 if (inner < minInnerIndex) //The value does not yet exist 450 { 451 const Index previousProfile = m_data.lowerProfile(outer); 452 m_data.lowerProfile(outer) = outer - inner; 453 454 const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile; 455 //shift data stored after this new one 456 const Index stop = m_rowStartIndex[rows()]; 457 const Index start = m_rowStartIndex[outer]; 458 459 460 for (Index innerIdx = stop; innerIdx >= start; innerIdx--) { 461 m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx); 462 } 463 464 for (Index innerIdx = rows(); innerIdx > outer; innerIdx--) { 465 m_rowStartIndex[innerIdx] += bandIncrement; 466 } 467 468 //zeros new data 469 memset(this->_lowerPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar)); 470 return m_data.lower(m_rowStartIndex[outer]); 471 } else { 472 return m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer))); 473 } 474 } 475 } else { 476 if (outer > inner) //upper matrix 477 { 478 const Index maxOuterIndex = inner + m_data.upperProfile(inner); 479 if (outer > maxOuterIndex) //The value does not yet exist 480 { 481 const Index previousProfile = m_data.upperProfile(inner); 482 m_data.upperProfile(inner) = outer - inner; 483 484 const Index bandIncrement = m_data.upperProfile(inner) - previousProfile; 485 //shift data stored after this new one 486 const Index stop = m_rowStartIndex[rows()]; 487 const Index start = m_rowStartIndex[inner + 1]; 488 489 for (Index innerIdx = stop; innerIdx >= start; innerIdx--) { 490 m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx); 491 } 492 493 for (Index innerIdx = inner + 1; innerIdx < outerSize() + 1; innerIdx++) { 494 m_rowStartIndex[innerIdx] += bandIncrement; 495 } 496 memset(this->_upperPtr() + m_rowStartIndex[inner] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar)); 497 return m_data.upper(m_rowStartIndex[inner] + m_data.upperProfile(inner)); 498 } else { 499 return m_data.upper(m_rowStartIndex[inner] + (outer - inner)); 500 } 501 } 502 503 if (outer < inner) //lower matrix 504 { 505 const Index maxInnerIndex = outer + m_data.lowerProfile(outer); 506 if (inner > maxInnerIndex) //The value does not yet exist 507 { 508 const Index previousProfile = m_data.lowerProfile(outer); 509 m_data.lowerProfile(outer) = inner - outer; 510 511 const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile; 512 //shift data stored after this new one 513 const Index stop = m_colStartIndex[cols()]; 514 const Index start = m_colStartIndex[outer + 1]; 515 516 for (Index innerIdx = stop; innerIdx >= start; innerIdx--) { 517 m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx); 518 } 519 520 for (Index innerIdx = outer + 1; innerIdx < outerSize() + 1; innerIdx++) { 521 m_colStartIndex[innerIdx] += bandIncrement; 522 } 523 memset(this->_lowerPtr() + m_colStartIndex[outer] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar)); 524 return m_data.lower(m_colStartIndex[outer] + m_data.lowerProfile(outer)); 525 } else { 526 return m_data.lower(m_colStartIndex[outer] + (inner - outer)); 527 } 528 } 529 } 530 } 531 532 /** Must be called after inserting a set of non zero entries. 533 */ 534 inline void finalize() { 535 if (IsRowMajor) { 536 if (rows() > cols()) 537 m_data.resize(cols(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1); 538 else 539 m_data.resize(rows(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1); 540 541 // eigen_assert(rows() == cols() && "memory reorganisatrion only works with suare matrix"); 542 // 543 // Scalar* newArray = new Scalar[m_colStartIndex[cols()] + 1 + m_rowStartIndex[rows()] + 1]; 544 // Index dataIdx = 0; 545 // for (Index row = 0; row < rows(); row++) { 546 // 547 // const Index nbLowerElts = m_rowStartIndex[row + 1] - m_rowStartIndex[row]; 548 // // std::cout << "nbLowerElts" << nbLowerElts << std::endl; 549 // memcpy(newArray + dataIdx, m_data.m_lower + m_rowStartIndex[row], nbLowerElts * sizeof (Scalar)); 550 // m_rowStartIndex[row] = dataIdx; 551 // dataIdx += nbLowerElts; 552 // 553 // const Index nbUpperElts = m_colStartIndex[row + 1] - m_colStartIndex[row]; 554 // memcpy(newArray + dataIdx, m_data.m_upper + m_colStartIndex[row], nbUpperElts * sizeof (Scalar)); 555 // m_colStartIndex[row] = dataIdx; 556 // dataIdx += nbUpperElts; 557 // 558 // 559 // } 560 // //todo : don't access m_data profile directly : add an accessor from SkylineMatrix 561 // m_rowStartIndex[rows()] = m_rowStartIndex[rows()-1] + m_data.lowerProfile(rows()-1); 562 // m_colStartIndex[cols()] = m_colStartIndex[cols()-1] + m_data.upperProfile(cols()-1); 563 // 564 // delete[] m_data.m_lower; 565 // delete[] m_data.m_upper; 566 // 567 // m_data.m_lower = newArray; 568 // m_data.m_upper = newArray; 569 } else { 570 if (rows() > cols()) 571 m_data.resize(cols(), rows(), cols(), m_rowStartIndex[cols()] + 1, m_colStartIndex[cols()] + 1); 572 else 573 m_data.resize(rows(), rows(), cols(), m_rowStartIndex[rows()] + 1, m_colStartIndex[rows()] + 1); 574 } 575 } 576 577 inline void squeeze() { 578 finalize(); 579 m_data.squeeze(); 580 } 581 582 void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar > ()) { 583 //TODO 584 } 585 586 /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero 587 * \sa resizeNonZeros(Index), reserve(), setZero() 588 */ 589 void resize(size_t rows, size_t cols) { 590 const Index diagSize = rows > cols ? cols : rows; 591 m_innerSize = IsRowMajor ? cols : rows; 592 593 eigen_assert(rows == cols && "Skyline matrix must be square matrix"); 594 595 if (diagSize % 2) { // diagSize is odd 596 const Index k = (diagSize - 1) / 2; 597 598 m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols, 599 2 * k * k + k + 1, 600 2 * k * k + k + 1); 601 602 } else // diagSize is even 603 { 604 const Index k = diagSize / 2; 605 m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols, 606 2 * k * k - k + 1, 607 2 * k * k - k + 1); 608 } 609 610 if (m_colStartIndex && m_rowStartIndex) { 611 delete[] m_colStartIndex; 612 delete[] m_rowStartIndex; 613 } 614 m_colStartIndex = new Index [cols + 1]; 615 m_rowStartIndex = new Index [rows + 1]; 616 m_outerSize = diagSize; 617 618 m_data.reset(); 619 m_data.clear(); 620 621 m_outerSize = diagSize; 622 memset(m_colStartIndex, 0, (cols + 1) * sizeof (Index)); 623 memset(m_rowStartIndex, 0, (rows + 1) * sizeof (Index)); 624 } 625 626 void resizeNonZeros(Index size) { 627 m_data.resize(size); 628 } 629 630 inline SkylineMatrix() 631 : m_outerSize(-1), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { 632 resize(0, 0); 633 } 634 635 inline SkylineMatrix(size_t rows, size_t cols) 636 : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { 637 resize(rows, cols); 638 } 639 640 template<typename OtherDerived> 641 inline SkylineMatrix(const SkylineMatrixBase<OtherDerived>& other) 642 : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { 643 *this = other.derived(); 644 } 645 646 inline SkylineMatrix(const SkylineMatrix & other) 647 : Base(), m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) { 648 *this = other.derived(); 649 } 650 651 inline void swap(SkylineMatrix & other) { 652 //EIGEN_DBG_SKYLINE(std::cout << "SkylineMatrix:: swap\n"); 653 std::swap(m_colStartIndex, other.m_colStartIndex); 654 std::swap(m_rowStartIndex, other.m_rowStartIndex); 655 std::swap(m_innerSize, other.m_innerSize); 656 std::swap(m_outerSize, other.m_outerSize); 657 m_data.swap(other.m_data); 658 } 659 660 inline SkylineMatrix & operator=(const SkylineMatrix & other) { 661 std::cout << "SkylineMatrix& operator=(const SkylineMatrix& other)\n"; 662 if (other.isRValue()) { 663 swap(other.const_cast_derived()); 664 } else { 665 resize(other.rows(), other.cols()); 666 memcpy(m_colStartIndex, other.m_colStartIndex, (m_outerSize + 1) * sizeof (Index)); 667 memcpy(m_rowStartIndex, other.m_rowStartIndex, (m_outerSize + 1) * sizeof (Index)); 668 m_data = other.m_data; 669 } 670 return *this; 671 } 672 673 template<typename OtherDerived> 674 inline SkylineMatrix & operator=(const SkylineMatrixBase<OtherDerived>& other) { 675 const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); 676 if (needToTranspose) { 677 // TODO 678 // return *this; 679 } else { 680 // there is no special optimization 681 return SkylineMatrixBase<SkylineMatrix>::operator=(other.derived()); 682 } 683 } 684 685 friend std::ostream & operator <<(std::ostream & s, const SkylineMatrix & m) { 686 687 EIGEN_DBG_SKYLINE( 688 std::cout << "upper elements : " << std::endl; 689 for (Index i = 0; i < m.m_data.upperSize(); i++) 690 std::cout << m.m_data.upper(i) << "\t"; 691 std::cout << std::endl; 692 std::cout << "upper profile : " << std::endl; 693 for (Index i = 0; i < m.m_data.upperProfileSize(); i++) 694 std::cout << m.m_data.upperProfile(i) << "\t"; 695 std::cout << std::endl; 696 std::cout << "lower startIdx : " << std::endl; 697 for (Index i = 0; i < m.m_data.upperProfileSize(); i++) 698 std::cout << (IsRowMajor ? m.m_colStartIndex[i] : m.m_rowStartIndex[i]) << "\t"; 699 std::cout << std::endl; 700 701 702 std::cout << "lower elements : " << std::endl; 703 for (Index i = 0; i < m.m_data.lowerSize(); i++) 704 std::cout << m.m_data.lower(i) << "\t"; 705 std::cout << std::endl; 706 std::cout << "lower profile : " << std::endl; 707 for (Index i = 0; i < m.m_data.lowerProfileSize(); i++) 708 std::cout << m.m_data.lowerProfile(i) << "\t"; 709 std::cout << std::endl; 710 std::cout << "lower startIdx : " << std::endl; 711 for (Index i = 0; i < m.m_data.lowerProfileSize(); i++) 712 std::cout << (IsRowMajor ? m.m_rowStartIndex[i] : m.m_colStartIndex[i]) << "\t"; 713 std::cout << std::endl; 714 ); 715 for (Index rowIdx = 0; rowIdx < m.rows(); rowIdx++) { 716 for (Index colIdx = 0; colIdx < m.cols(); colIdx++) { 717 s << m.coeff(rowIdx, colIdx) << "\t"; 718 } 719 s << std::endl; 720 } 721 return s; 722 } 723 724 /** Destructor */ 725 inline ~SkylineMatrix() { 726 delete[] m_colStartIndex; 727 delete[] m_rowStartIndex; 728 } 729 730 /** Overloaded for performance */ 731 Scalar sum() const; 732 }; 733 734 template<typename Scalar, int _Options> 735 class SkylineMatrix<Scalar, _Options>::InnerUpperIterator { 736 public: 737 738 InnerUpperIterator(const SkylineMatrix& mat, Index outer) 739 : m_matrix(mat), m_outer(outer), 740 m_id(_Options == RowMajor ? mat.m_colStartIndex[outer] : mat.m_rowStartIndex[outer] + 1), 741 m_start(m_id), 742 m_end(_Options == RowMajor ? mat.m_colStartIndex[outer + 1] : mat.m_rowStartIndex[outer + 1] + 1) { 743 } 744 745 inline InnerUpperIterator & operator++() { 746 m_id++; 747 return *this; 748 } 749 750 inline InnerUpperIterator & operator+=(Index shift) { 751 m_id += shift; 752 return *this; 753 } 754 755 inline Scalar value() const { 756 return m_matrix.m_data.upper(m_id); 757 } 758 759 inline Scalar* valuePtr() { 760 return const_cast<Scalar*> (&(m_matrix.m_data.upper(m_id))); 761 } 762 763 inline Scalar& valueRef() { 764 return const_cast<Scalar&> (m_matrix.m_data.upper(m_id)); 765 } 766 767 inline Index index() const { 768 return IsRowMajor ? m_outer - m_matrix.m_data.upperProfile(m_outer) + (m_id - m_start) : 769 m_outer + (m_id - m_start) + 1; 770 } 771 772 inline Index row() const { 773 return IsRowMajor ? index() : m_outer; 774 } 775 776 inline Index col() const { 777 return IsRowMajor ? m_outer : index(); 778 } 779 780 inline size_t size() const { 781 return m_matrix.m_data.upperProfile(m_outer); 782 } 783 784 inline operator bool() const { 785 return (m_id < m_end) && (m_id >= m_start); 786 } 787 788 protected: 789 const SkylineMatrix& m_matrix; 790 const Index m_outer; 791 Index m_id; 792 const Index m_start; 793 const Index m_end; 794 }; 795 796 template<typename Scalar, int _Options> 797 class SkylineMatrix<Scalar, _Options>::InnerLowerIterator { 798 public: 799 800 InnerLowerIterator(const SkylineMatrix& mat, Index outer) 801 : m_matrix(mat), 802 m_outer(outer), 803 m_id(_Options == RowMajor ? mat.m_rowStartIndex[outer] : mat.m_colStartIndex[outer] + 1), 804 m_start(m_id), 805 m_end(_Options == RowMajor ? mat.m_rowStartIndex[outer + 1] : mat.m_colStartIndex[outer + 1] + 1) { 806 } 807 808 inline InnerLowerIterator & operator++() { 809 m_id++; 810 return *this; 811 } 812 813 inline InnerLowerIterator & operator+=(Index shift) { 814 m_id += shift; 815 return *this; 816 } 817 818 inline Scalar value() const { 819 return m_matrix.m_data.lower(m_id); 820 } 821 822 inline Scalar* valuePtr() { 823 return const_cast<Scalar*> (&(m_matrix.m_data.lower(m_id))); 824 } 825 826 inline Scalar& valueRef() { 827 return const_cast<Scalar&> (m_matrix.m_data.lower(m_id)); 828 } 829 830 inline Index index() const { 831 return IsRowMajor ? m_outer - m_matrix.m_data.lowerProfile(m_outer) + (m_id - m_start) : 832 m_outer + (m_id - m_start) + 1; 833 ; 834 } 835 836 inline Index row() const { 837 return IsRowMajor ? m_outer : index(); 838 } 839 840 inline Index col() const { 841 return IsRowMajor ? index() : m_outer; 842 } 843 844 inline size_t size() const { 845 return m_matrix.m_data.lowerProfile(m_outer); 846 } 847 848 inline operator bool() const { 849 return (m_id < m_end) && (m_id >= m_start); 850 } 851 852 protected: 853 const SkylineMatrix& m_matrix; 854 const Index m_outer; 855 Index m_id; 856 const Index m_start; 857 const Index m_end; 858 }; 859 860 } // end namespace Eigen 861 862 #endif // EIGEN_SkylineMatrix_H 863