1 /////////////////////////////////////////////////////////////////////////////// 2 // // 3 // The Template Matrix/Vector Library for C++ was created by Mike Jarvis // 4 // Copyright (C) 1998 - 2016 // 5 // All rights reserved // 6 // // 7 // The project is hosted at https://code.google.com/p/tmv-cpp/ // 8 // where you can find the current version and current documention. // 9 // // 10 // For concerns or problems with the software, Mike may be contacted at // 11 // mike_jarvis17 [at] gmail. // 12 // // 13 // This software is licensed under a FreeBSD license. The file // 14 // TMV_LICENSE should have bee included with this distribution. // 15 // It not, you can get a copy from https://code.google.com/p/tmv-cpp/. // 16 // // 17 // Essentially, you can use this software however you want provided that // 18 // you include the TMV_LICENSE file in any distribution that uses it. // 19 // // 20 /////////////////////////////////////////////////////////////////////////////// 21 22 23 //--------------------------------------------------------------------------- 24 // 25 // This file defines the TMV DiagMatrix class. 26 // 27 // The DiagMatrix class is provided for efficient storage of a diagonal 28 // matrix. You can do most of the things that you can do with a 29 // regular Matrix, but it will do them more efficiently. 30 // 31 // Constructors: 32 // 33 // DiagMatrix<T,A>(int size) 34 // Makes a DiagMatrix with column size and row size = size 35 // with _uninitialized_ values 36 // 37 // DiagMatrix<T,A>(int size, T x) 38 // Makes a DiagMatrix of size n with all values = x 39 // 40 // DiagMatrix<T,A>(const Vector<T>& vv) 41 // Make a DiagMatrix which copies the elements of vv. 42 // 43 // 44 // Access Functions 45 // 46 // int colsize() const 47 // int rowsize() const 48 // int size() const 49 // Return the dimensions of the DiagMatrix 50 // 51 // T& operator()(int i) 52 // T operator()(int i) const 53 // T& operator()(int i, int j) 54 // T operator()(int i, int j) const 55 // Return the (i,j) element of the DiagMatrix 56 // For the single paramter version, j=i 57 // 58 // VectorView diag() 59 // ConstVectorView diag() const 60 // Return the diagonal of the DiagMatrix as a VectorView 61 // 62 // 63 // Modifying Functions - The same as the regular Matrix counterparts 64 // 65 // DiagMatrix& setZero() 66 // DiagMatrix& setAllTo(T x) 67 // DiagMatrix& addToAll(T x) 68 // DiagMatrix<T>& transposeSelf() 69 // (Does nothing.) 70 // DiagMatrix& conjugateSelf() 71 // DiagMatrix& setToIdentity(x = 1) 72 // void Swap(DiagMatrix& m1, DiagMatrix& m2) 73 // 74 // 75 // subDiagMatrix: 76 // 77 // DiagMatrixView subDiagMatrix(int i1, int i2, int istep=1) 78 // Returns a Sub-DiagMatrix which extends from i1 to i2 (step istep) 79 // which refers to the same physical elements as the original. 80 // As usual, i2 is the "one past the end" element. 81 // 82 // DiagMatrixView realPart() 83 // DiagMatrixView imagPart() 84 // Returns a view to the real/imag elements of a complex DiagMatrix. 85 // 86 // DiagMatrixView view() 87 // DiagMatrixView transpose() 88 // Returns a view of a DiagMatrix. 89 // 90 // DiagMatrixView conjugate(m) 91 // DiagMatrixView adjoint(m) 92 // Returns a view to the conjugate of a DiagMatrix. 93 // 94 // Functions of DiagMatrices - Same as for regular Matrices: 95 // 96 // m.det() or Det(m) 97 // m.logDet() or m.logDet(T* sign) or LogDet(m) 98 // m.trace() or Trace(m) 99 // m.sumElements() or SumElements(m) 100 // m.sumAbsElements() or SumAbsElements(m) 101 // m.sumAbs2Elements() or SumAbs2Elements(m) 102 // m.norm() or m.normF() or Norm(m) or NormF(m) 103 // m.normSq() or NormSq(m) 104 // m.norm1() or Norm1(m) 105 // m.norm2() or Norm2(m) 106 // m.normInf() or NormInf(m) 107 // m.maxAbsElement() or MaxAbsElements(m) 108 // m.maxAbs2Element() or MaxAbs2Elements(m) 109 // (Note - for diagonal matrices, 110 // norm1 = norm2 = normInf = maxAbsElement.) 111 // 112 // m.inverse() or Inverse(m) 113 // m.invertSelf() 114 // m.makeInverse(minv) (takes either Matrix or DiagMatrix argument) 115 // m.makeInverseATA(invata) (takes either Matrix or DiagMatrix argument) 116 // 117 // I/O: 118 // 119 // os << d 120 // Writes d to ostream os as a full matrix 121 // 122 // os << CompactIO() << d 123 // Writes only the diagonal Vector to os 124 // 125 // is >> d 126 // is >> CompactIO() >> d 127 // Reads in d in either format 128 // 129 // 130 131 #ifndef TMV_DiagMatrix_H 132 #define TMV_DiagMatrix_H 133 134 #include "tmv/TMV_BaseDiagMatrix.h" 135 #include "tmv/TMV_BaseTriMatrix.h" 136 #include "tmv/TMV_Vector.h" 137 #include "tmv/TMV_TriMatrix.h" 138 #include "tmv/TMV_Matrix.h" 139 #include <vector> 140 141 namespace tmv { 142 143 template <typename T> 144 class GenDiagMatrix : 145 virtual public AssignableToDiagMatrix<T>, 146 virtual public AssignableToUpperTriMatrix<T>, 147 virtual public AssignableToLowerTriMatrix<T>, 148 public BaseMatrix<T> 149 { 150 public: 151 152 typedef TMV_RealType(T) RT; 153 typedef TMV_ComplexType(T) CT; 154 typedef T value_type; 155 typedef RT real_type; 156 typedef CT complex_type; 157 typedef GenDiagMatrix<T> type; 158 typedef DiagMatrix<T> copy_type; 159 typedef ConstDiagMatrixView<T> const_view_type; 160 typedef const_view_type const_transpose_type; 161 typedef const_view_type const_conjugate_type; 162 typedef const_view_type const_adjoint_type; 163 typedef ConstDiagMatrixView<RT> const_realpart_type; 164 typedef ConstVectorView<T> const_vec_type; 165 typedef DiagMatrixView<T> nonconst_type; 166 typedef typename const_vec_type::const_iterator const_iterator; 167 168 // 169 // Constructors 170 // 171 172 inline GenDiagMatrix<T>() {} GenDiagMatrix(const type &)173 inline GenDiagMatrix(const type&) {} ~GenDiagMatrix()174 virtual inline ~GenDiagMatrix() {} 175 176 // 177 // Access Functions 178 // 179 size()180 inline ptrdiff_t size() const { return diag().size(); } colsize()181 inline ptrdiff_t colsize() const { return size(); } rowsize()182 inline ptrdiff_t rowsize() const { return size(); } dt()183 inline DiagType dt() const { return NonUnitDiag; } 184 operator()185 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 186 { 187 TMVAssert(i>=0 && i<size()); 188 TMVAssert(j>=0 && j<size()); 189 if (i==j) return diag()(i); 190 else return T(0); 191 } 192 operator()193 inline T operator()(ptrdiff_t i) const 194 { 195 TMVAssert(i>=0 && i<size()); 196 return diag()(i); 197 } 198 diag()199 inline const_vec_type diag() const { return cdiag(); } 200 201 template <typename T2> isSameAs(const BaseMatrix<T2> &)202 inline bool isSameAs(const BaseMatrix<T2>& ) const 203 { return false; } 204 isSameAs(const GenDiagMatrix<T> & m2)205 inline bool isSameAs(const GenDiagMatrix<T>& m2) const 206 { 207 if (this == &m2) return true; 208 else return (diag().isSameAs(m2.diag())); 209 } 210 assignToM(MatrixView<RT> m2)211 inline void assignToM(MatrixView<RT> m2) const 212 { 213 TMVAssert(m2.colsize() == size()); 214 TMVAssert(m2.rowsize() == size()); 215 TMVAssert(isReal(T())); 216 m2.diag() = diag(); 217 m2.upperTri().offDiag().setZero(); 218 m2.lowerTri().offDiag().setZero(); 219 } assignToM(MatrixView<CT> m2)220 inline void assignToM(MatrixView<CT> m2) const 221 { 222 TMVAssert(m2.colsize() == size()); 223 TMVAssert(m2.rowsize() == size()); 224 m2.diag() = diag(); 225 m2.upperTri().offDiag().setZero(); 226 m2.lowerTri().offDiag().setZero(); 227 } 228 assignToU(UpperTriMatrixView<RT> m2)229 inline void assignToU(UpperTriMatrixView<RT> m2) const 230 { 231 TMVAssert(m2.size() == size()); 232 TMVAssert(isReal(T())); 233 m2.diag() = diag(); 234 m2.offDiag().setZero(); 235 } assignToU(UpperTriMatrixView<CT> m2)236 inline void assignToU(UpperTriMatrixView<CT> m2) const 237 { 238 TMVAssert(m2.size() == size()); 239 m2.diag() = diag(); 240 m2.offDiag().setZero(); 241 } 242 assignToL(LowerTriMatrixView<RT> m2)243 inline void assignToL(LowerTriMatrixView<RT> m2) const 244 { 245 TMVAssert(m2.size() == size()); 246 TMVAssert(isReal(T())); 247 m2.diag() = diag(); 248 m2.offDiag().setZero(); 249 } assignToL(LowerTriMatrixView<CT> m2)250 inline void assignToL(LowerTriMatrixView<CT> m2) const 251 { 252 TMVAssert(m2.size() == size()); 253 m2.diag() = diag(); 254 m2.offDiag().setZero(); 255 } 256 assignToD(DiagMatrixView<RT> m2)257 inline void assignToD(DiagMatrixView<RT> m2) const 258 { 259 TMVAssert(m2.size() == size()); 260 TMVAssert(isReal(T())); 261 m2.diag() = diag(); 262 } assignToD(DiagMatrixView<CT> m2)263 inline void assignToD(DiagMatrixView<CT> m2) const 264 { 265 TMVAssert(m2.size() == size()); 266 m2.diag() = diag(); 267 } 268 269 // 270 // subDiagMatrix 271 // 272 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)273 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 274 { return const_view_type(diag().cSubVector(i1,i2)); } 275 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)276 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 277 { 278 TMVAssert(diag().hasSubVector(i1,i2,1)); 279 return cSubDiagMatrix(i1,i2); 280 } 281 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)282 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 283 { return const_view_type(diag().cSubVector(i1,i2,istep)); } 284 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)285 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 286 { 287 TMVAssert(diag().hasSubVector(i1,i2,istep)); 288 return cSubDiagMatrix(i1,i2,istep); 289 } 290 realPart()291 inline const_realpart_type realPart() const 292 { return const_realpart_type(diag().realPart()); } 293 imagPart()294 inline const_realpart_type imagPart() const 295 { return const_realpart_type(diag().imagPart()); } 296 view()297 inline const_view_type view() const 298 { return const_view_type(diag()); } 299 transpose()300 inline const_view_type transpose() const 301 { return view(); } 302 conjugate()303 inline const_view_type conjugate() const 304 { return const_view_type(diag().conjugate()); } 305 adjoint()306 inline const_view_type adjoint() const 307 { return conjugate(); } 308 nonConst()309 inline nonconst_type nonConst() const 310 { return nonconst_type(diag().nonConst()); } 311 312 313 // 314 // Functions of DiagMatrix 315 // 316 317 T det() const; 318 319 RT logDet(T* sign=0) const; 320 trace()321 inline T trace() const 322 { return diag().sumElements(); } 323 sumElements()324 inline T sumElements() const 325 { return diag().sumElements(); } 326 sumAbsElements()327 inline RT sumAbsElements() const 328 { return diag().sumAbsElements(); } 329 sumAbs2Elements()330 inline RT sumAbs2Elements() const 331 { return diag().sumAbs2Elements(); } 332 norm()333 inline RT norm() const 334 { return normF(); } 335 normF()336 inline RT normF() const 337 { return diag().norm(); } 338 339 inline RT normSq(RT scale = RT(1)) const 340 { return diag().normSq(scale); } 341 norm1()342 inline RT norm1() const 343 { return diag().maxAbsElement(); } 344 doNorm2()345 inline RT doNorm2() const 346 { return diag().maxAbsElement(); } 347 doCondition()348 inline RT doCondition() const 349 { return diag().maxAbsElement()/diag().minAbsElement(); } 350 norm2()351 inline RT norm2() const 352 { return doNorm2(); } 353 condition()354 inline RT condition() const 355 { return doCondition(); } 356 normInf()357 inline RT normInf() const 358 { return diag().maxAbsElement(); } 359 maxAbsElement()360 inline RT maxAbsElement() const 361 { return diag().maxAbsElement(); } 362 maxAbs2Element()363 inline RT maxAbs2Element() const 364 { return diag().maxAbs2Element(); } 365 isSingular()366 inline bool isSingular() const 367 { return diag().minAbsElement() == RT(0); } 368 369 template <typename T1> 370 void doMakeInverse(MatrixView<T1> minv) const; 371 372 template <typename T1> 373 void doMakeInverse(DiagMatrixView<T1> minv) const; 374 375 void doMakeInverseATA(MatrixView<T> ata) const; 376 void doMakeInverseATA(DiagMatrixView<T> ata) const; 377 makeInverse(MatrixView<T> minv)378 inline void makeInverse(MatrixView<T> minv) const 379 { 380 TMVAssert(minv.colsize() == size()); 381 TMVAssert(minv.rowsize() == size()); 382 doMakeInverse(minv); 383 } 384 385 template <typename T1> makeInverse(MatrixView<T1> minv)386 inline void makeInverse(MatrixView<T1> minv) const 387 { 388 TMVAssert(minv.colsize() == size()); 389 TMVAssert(minv.rowsize() == size()); 390 doMakeInverse(minv); 391 } 392 393 template <typename T1> makeInverse(DiagMatrixView<T1> minv)394 inline void makeInverse(DiagMatrixView<T1> minv) const 395 { 396 TMVAssert(minv.size() == size()); 397 doMakeInverse(minv); 398 } 399 400 QuotXD<T,T> QInverse() const; inverse()401 inline QuotXD<T,T> inverse() const 402 { return QInverse(); } 403 makeInverseATA(MatrixView<T> ata)404 inline void makeInverseATA(MatrixView<T> ata) const 405 { 406 TMVAssert(ata.colsize() == size()); 407 TMVAssert(ata.rowsize() == size()); 408 doMakeInverseATA(ata); 409 } makeInverseATA(DiagMatrixView<T> ata)410 inline void makeInverseATA(DiagMatrixView<T> ata) const 411 { 412 TMVAssert(ata.size() == size()); 413 doMakeInverseATA(ata); 414 } 415 416 template <typename T1, int A> makeInverse(DiagMatrix<T1,A> & minv)417 inline void makeInverse(DiagMatrix<T1,A>& minv) const 418 { makeInverse(minv.view()); } 419 420 template <typename T1, int A> makeInverse(Matrix<T1,A> & minv)421 inline void makeInverse(Matrix<T1,A>& minv) const 422 { makeInverse(minv.view()); } 423 424 template <int A> makeInverseATA(DiagMatrix<T,A> & minv)425 inline void makeInverseATA(DiagMatrix<T,A>& minv) const 426 { makeInverseATA(minv.view()); } 427 428 template <int A> makeInverseATA(Matrix<T,A> & minv)429 inline void makeInverseATA(Matrix<T,A>& minv) const 430 { makeInverseATA(minv.view()); } 431 432 433 // 434 // I/O 435 // 436 437 void write(const TMV_Writer& writer) const; 438 439 // 440 // Arithmetic Helpers 441 // 442 443 template <typename T1> 444 void doLDivEq(VectorView<T1> v) const; 445 446 template <typename T1, typename T0> 447 void doLDiv(const GenVector<T1>& v1, VectorView<T0> v0) const; 448 449 template <typename T1> 450 void doLDivEq(MatrixView<T1> m) const; 451 452 template <typename T1, typename T0> 453 void doLDiv(const GenMatrix<T1>& m1, MatrixView<T0> m0) const; 454 455 template <typename T1> LDivEq(VectorView<T1> v)456 inline void LDivEq(VectorView<T1> v) const 457 { 458 TMVAssert(v.size() == size()); 459 doLDivEq(v); 460 } 461 462 template <typename T1, typename T0> LDiv(const GenVector<T1> & v1,VectorView<T0> v0)463 inline void LDiv( 464 const GenVector<T1>& v1, VectorView<T0> v0) const 465 { 466 TMVAssert(v1.size() == size()); 467 TMVAssert(v0.size() == size()); 468 doLDiv(v1,v0); 469 } 470 471 template <typename T1> RDivEq(VectorView<T1> v)472 inline void RDivEq(VectorView<T1> v) const 473 { LDivEq(v); } 474 475 template <typename T1, typename T0> RDiv(const GenVector<T1> & v1,VectorView<T0> v0)476 inline void RDiv( 477 const GenVector<T1>& v1, VectorView<T0> v0) const 478 { LDiv(v1,v0); } 479 480 template <typename T1> LDivEq(MatrixView<T1> m)481 inline void LDivEq(MatrixView<T1> m) const 482 { 483 TMVAssert(m.colsize() == size()); 484 doLDivEq(m); 485 } 486 template <typename T1, typename T0> LDiv(const GenMatrix<T1> & m1,MatrixView<T0> m0)487 inline void LDiv( 488 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 489 { 490 TMVAssert(m1.colsize() == size()); 491 TMVAssert(m0.colsize() == size()); 492 TMVAssert(m1.rowsize() == m0.rowsize()); 493 doLDiv(m1,m0); 494 } 495 template <typename T1> RDivEq(MatrixView<T1> m)496 inline void RDivEq(MatrixView<T1> m) const 497 { LDivEq(m.transpose()); } 498 499 template <typename T1, typename T0> RDiv(const GenMatrix<T1> & m1,MatrixView<T0> m0)500 void RDiv(const GenMatrix<T1>& m1, MatrixView<T0> m0) const 501 { LDiv(m1.transpose(),m0.transpose()); } 502 503 template <typename T1> DivEq(DiagMatrixView<T1> m)504 void DivEq(DiagMatrixView<T1> m) const 505 { LDivEq(m.diag()); } 506 507 template <typename T1, typename T0> Div(const GenDiagMatrix<T1> & m1,DiagMatrixView<T0> m0)508 void Div( 509 const GenDiagMatrix<T1>& m1, DiagMatrixView<T0> m0) const 510 { LDiv(m1.diag(),m0.diag()); } 511 512 // For easier compatibility with regular matrices: divideInPlace()513 inline void divideInPlace() const {} saveDiv()514 inline void saveDiv() const {} setDiv()515 inline void setDiv() const {} unsetDiv()516 inline void unsetDiv() const {} resetDiv()517 inline void resetDiv() const {} divIsSet()518 inline bool divIsSet() const { return true; } divideUsing(DivType TMV_DEBUGPARAM (dt))519 inline void divideUsing(DivType TMV_DEBUGPARAM(dt)) const 520 { TMVAssert(dt == LU); } 521 inline bool checkDecomp(std::ostream* fout=0) const { return true; } 522 inline bool checkDecomp( 523 const BaseMatrix<T>& m2, std::ostream* fout=0) const 524 { return true; } 525 cref(ptrdiff_t i,ptrdiff_t j)526 inline T cref(ptrdiff_t i, ptrdiff_t j) const 527 { return i==j ? cdiag()(i) : 0; } cptr()528 inline const T* cptr() const 529 { return cdiag().cptr(); } step()530 inline ptrdiff_t step() const 531 { return cdiag().step(); } isconj()532 inline bool isconj() const 533 { return cdiag().isconj(); } 534 begin()535 inline const_iterator begin() const 536 { return diag().begin(); } end()537 inline const_iterator end() const 538 { return diag().end(); } 539 540 protected : 541 542 virtual ConstVectorView<T> cdiag() const = 0; 543 544 private : 545 546 type& operator=(const type&); 547 548 }; // GenDiagMatrix 549 550 template <typename T, int A> 551 class ConstDiagMatrixView : public GenDiagMatrix<T> 552 { 553 public : 554 555 typedef GenDiagMatrix<T> base; 556 typedef ConstDiagMatrixView<T,A> type; 557 typedef ConstVectorView<T> const_vec_type; 558 ConstDiagMatrixView(const type & rhs)559 inline ConstDiagMatrixView(const type& rhs) : itsdiag(rhs.cdiag()) 560 { TMVAssert(Attrib<A>::viewok); } 561 ConstDiagMatrixView(const base & rhs)562 inline ConstDiagMatrixView(const base& rhs) : itsdiag(rhs.diag()) 563 { TMVAssert(Attrib<A>::viewok); } 564 ConstDiagMatrixView(const GenVector<T> & v)565 explicit inline ConstDiagMatrixView(const GenVector<T>& v) : itsdiag(v) 566 { TMVAssert(Attrib<A>::viewok); } 567 ~ConstDiagMatrixView()568 virtual inline ~ConstDiagMatrixView() {} 569 570 protected : 571 572 ConstVectorView<T> itsdiag; cdiag()573 inline ConstVectorView<T> cdiag() const { return itsdiag; } 574 575 private : 576 577 type& operator=(const type&); 578 579 }; // ConstDiagMatrixView 580 581 template <typename T> 582 class ConstDiagMatrixView<T,FortranStyle> : 583 public ConstDiagMatrixView<T,CStyle> 584 { 585 public : 586 587 typedef TMV_RealType(T) RT; 588 typedef GenDiagMatrix<T> base; 589 typedef ConstDiagMatrixView<T,FortranStyle> type; 590 typedef ConstVectorView<T,FortranStyle> const_vec_type; 591 typedef ConstDiagMatrixView<T,FortranStyle> const_view_type; 592 typedef const_view_type const_transpose_type; 593 typedef const_view_type const_conjugate_type; 594 typedef const_view_type const_adjoint_type; 595 typedef ConstDiagMatrixView<RT,FortranStyle> const_realpart_type; 596 typedef ConstDiagMatrixView<T,CStyle> c_type; 597 typedef DiagMatrixView<T,FortranStyle> nonconst_type; 598 ConstDiagMatrixView(const type & rhs)599 inline ConstDiagMatrixView(const type& rhs) : c_type(rhs) {} 600 ConstDiagMatrixView(const base & rhs)601 inline ConstDiagMatrixView(const base& rhs) : c_type(rhs) {} 602 ConstDiagMatrixView(const GenVector<T> & v)603 explicit inline ConstDiagMatrixView(const GenVector<T>& v) : 604 c_type(v) {} 605 ~ConstDiagMatrixView()606 virtual inline ~ConstDiagMatrixView() {} 607 608 // 609 // Access Functions 610 // 611 operator()612 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 613 { 614 TMVAssert(i>0 && i<=c_type::size()); 615 TMVAssert(j>0 && j<=c_type::size()); 616 if (i==j) return diag()(i); 617 else return T(0); 618 } 619 operator()620 inline T operator()(ptrdiff_t i) const 621 { 622 TMVAssert(i>0 && i<=c_type::size()); 623 return diag()(i); 624 } 625 diag()626 inline const_vec_type diag() const 627 { return const_vec_type(c_type::diag()); } 628 629 // 630 // subDiagMatrix 631 // 632 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)633 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 634 { 635 TMVAssert(diag().hasSubVector(i1,i2,1)); 636 return base::cSubDiagMatrix(i1-1,i2); 637 } 638 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)639 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 640 { 641 TMVAssert(diag().hasSubVector(i1,i2,istep)); 642 return base::cSubDiagMatrix(i1-1,i2-1+istep,istep); 643 } 644 realPart()645 inline const_realpart_type realPart() const 646 { return const_realpart_type(diag().realPart()); } 647 imagPart()648 inline const_realpart_type imagPart() const 649 { return const_realpart_type(diag().imagPart()); } 650 view()651 inline const_view_type view() const 652 { return const_view_type(diag()); } 653 transpose()654 inline const_view_type transpose() const 655 { return view(); } 656 conjugate()657 inline const_view_type conjugate() const 658 { return const_view_type(diag().conjugate()); } 659 adjoint()660 inline const_view_type adjoint() const 661 { return conjugate(); } 662 nonConst()663 inline nonconst_type nonConst() const 664 { return nonconst_type(diag().nonConst()); } 665 666 type& operator=(const type&); 667 668 }; // ConstDiagMatrixView - FortranStyle 669 670 template <typename T, int A> 671 class DiagMatrixView : public GenDiagMatrix<T> 672 { 673 public: 674 675 typedef TMV_RealType(T) RT; 676 typedef TMV_ComplexType(T) CT; 677 typedef GenDiagMatrix<T> base; 678 typedef DiagMatrixView<T,A> type; 679 typedef DiagMatrixView<T,A> view_type; 680 typedef view_type transpose_type; 681 typedef view_type conjugate_type; 682 typedef view_type adjoint_type; 683 typedef DiagMatrixView<RT,A> realpart_type; 684 typedef VectorView<T,A> vec_type; 685 typedef ConstDiagMatrixView<T,A> const_view_type; 686 typedef const_view_type const_transpose_type; 687 typedef const_view_type const_conjugate_type; 688 typedef const_view_type const_adjoint_type; 689 typedef ConstDiagMatrixView<RT,A> const_realpart_type; 690 typedef ConstVectorView<T,A> const_vec_type; 691 typedef typename RefHelper<T>::reference reference; 692 typedef typename vec_type::iterator iterator; 693 typedef typename const_vec_type::const_iterator const_iterator; 694 695 // 696 // Constructors 697 // 698 DiagMatrixView(const type & rhs)699 inline DiagMatrixView(const type& rhs) : itsdiag(rhs.itsdiag) 700 { TMVAssert(Attrib<A>::diagmatrixok); } 701 DiagMatrixView(vec_type _diag)702 explicit inline DiagMatrixView(vec_type _diag) : itsdiag(_diag) 703 { TMVAssert(Attrib<A>::diagmatrixok); } 704 ~DiagMatrixView()705 virtual inline ~DiagMatrixView() {} 706 707 // 708 // Op= 709 // 710 711 inline type& operator=(const type& m2) 712 { m2.assignToD(*this); return *this; } 713 714 inline type& operator=(const GenDiagMatrix<RT>& m2) 715 { m2.assignToD(*this); return *this; } 716 717 inline type& operator=(const GenDiagMatrix<CT>& m2) 718 { m2.assignToD(*this); return *this; } 719 720 template <typename T2> 721 inline type& operator=(const GenDiagMatrix<T2>& m2) 722 { itsdiag = m2.diag(); return *this; } 723 724 inline type& operator=(const T& x) 725 { return setToIdentity(x); } 726 727 inline type& operator=(const AssignableToDiagMatrix<RT>& m2) 728 { 729 TMVAssert(size() == m2.size()); 730 m2.assignToD(*this); 731 return *this; 732 } 733 734 inline type& operator=(const AssignableToDiagMatrix<CT>& m2) 735 { 736 TMVAssert(size() == m2.size()); 737 m2.assignToD(*this); 738 return *this; 739 } 740 741 typedef ListAssigner<T,iterator> MyListAssigner; 742 inline MyListAssigner operator<<(const T& x) 743 { return MyListAssigner(begin(),size(),x); } 744 745 746 // 747 // Access 748 // 749 operator()750 inline reference operator()(ptrdiff_t i) 751 { 752 TMVAssert(i>=0 && i<size()); 753 return diag()(i); 754 } operator()755 inline reference operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j)) 756 { 757 TMVAssert(i>=0 && i<size()); 758 TMVAssert(i==j); 759 return diag()(i); 760 } 761 diag()762 inline vec_type diag() { return itsdiag; } 763 764 // Repeat const versions operator()765 inline T operator()(ptrdiff_t i) const 766 { return base::operator()(i); } operator()767 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 768 { return base::operator()(i,j); } diag()769 inline const_vec_type diag() const 770 { return base::diag(); } 771 772 // 773 // Modifying Functions 774 // 775 setZero()776 inline type& setZero() 777 { diag().setZero(); return *this; } 778 setAllTo(const T & x)779 inline type& setAllTo(const T& x) 780 { diag().setAllTo(x); return *this; } 781 addToAll(const T & x)782 inline type& addToAll(const T& x) 783 { diag().addToAll(x); return *this; } 784 clip(RT thresh)785 inline type& clip(RT thresh) 786 { diag().clip(thresh); return *this; } 787 transposeSelf()788 inline type& transposeSelf() 789 { return *this; } 790 conjugateSelf()791 inline type& conjugateSelf() 792 { diag().conjugateSelf(); return *this; } 793 794 type& invertSelf(); 795 796 inline type& setToIdentity(const T& x=T(1)) 797 { return setAllTo(x); } 798 799 // 800 // subDiagMatrix 801 // 802 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)803 inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) 804 { return view_type(itsdiag.cSubVector(i1,i2)); } 805 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)806 inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) 807 { 808 TMVAssert(diag().hasSubVector(i1,i2,1)); 809 return cSubDiagMatrix(i1,i2); 810 } 811 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)812 inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 813 { return view_type(itsdiag.cSubVector(i1,i2,istep)); } 814 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)815 inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 816 { 817 TMVAssert(diag().hasSubVector(i1,i2,istep)); 818 return cSubDiagMatrix(i1,i2,istep); 819 } 820 realPart()821 inline realpart_type realPart() 822 { return realpart_type(diag().realPart()); } 823 imagPart()824 inline realpart_type imagPart() 825 { return realpart_type(diag().imagPart()); } 826 view()827 inline view_type view() 828 { return *this; } 829 transpose()830 inline view_type transpose() 831 { return *this; } 832 conjugate()833 inline view_type conjugate() 834 { return view_type(diag().conjugate()); } 835 adjoint()836 inline view_type adjoint() 837 { return conjugate(); } 838 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)839 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 840 { return base::cSubDiagMatrix(i1,i2); } subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)841 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 842 { return base::subDiagMatrix(i1,i2); } cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)843 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 844 { return base::cSubDiagMatrix(i1,i2,istep); } subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)845 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 846 { return base::subDiagMatrix(i1,i2,istep); } realPart()847 inline const_realpart_type realPart() const 848 { return base::realPart(); } imagPart()849 inline const_realpart_type imagPart() const 850 { return base::imagPart(); } view()851 inline const_view_type view() const 852 { return base::view(); } transpose()853 inline const_view_type transpose() const 854 { return base::transpose(); } conjugate()855 inline const_view_type conjugate() const 856 { return base::conjugate(); } adjoint()857 inline const_view_type adjoint() const 858 { return base::adjoint(); } 859 860 // 861 // I/O 862 // 863 864 void read(const TMV_Reader& reader); 865 866 using GenDiagMatrix<T>::size; 867 begin()868 inline iterator begin() 869 { return diag().begin(); } end()870 inline iterator end() 871 { return diag().end(); } 872 begin()873 inline const_iterator begin() const 874 { return diag().begin(); } end()875 inline const_iterator end() const 876 { return diag().end(); } 877 878 protected: 879 880 VectorView<T> itsdiag; cdiag()881 inline ConstVectorView<T> cdiag() const { return itsdiag; } 882 883 }; // DiagMatrixView 884 885 template <typename T> 886 class DiagMatrixView<T,FortranStyle> : public DiagMatrixView<T,CStyle> 887 { 888 public: 889 890 typedef TMV_RealType(T) RT; 891 typedef TMV_ComplexType(T) CT; 892 typedef GenDiagMatrix<T> base; 893 typedef DiagMatrixView<T,FortranStyle> type; 894 typedef DiagMatrixView<T,CStyle> c_type; 895 typedef DiagMatrixView<T,FortranStyle> view_type; 896 typedef view_type transpose_type; 897 typedef view_type conjugate_type; 898 typedef view_type adjoint_type; 899 typedef DiagMatrixView<RT,FortranStyle> realpart_type; 900 typedef VectorView<T,FortranStyle> vec_type; 901 typedef ConstDiagMatrixView<T,FortranStyle> const_view_type; 902 typedef const_view_type const_transpose_type; 903 typedef const_view_type const_conjugate_type; 904 typedef const_view_type const_adjoint_type; 905 typedef ConstDiagMatrixView<RT,FortranStyle> const_realpart_type; 906 typedef ConstVectorView<T,FortranStyle> const_vec_type; 907 typedef typename RefHelper<T>::reference reference; 908 909 // 910 // Constructors 911 // 912 DiagMatrixView(const type & rhs)913 inline DiagMatrixView(const type& rhs) : c_type(rhs) {} 914 DiagMatrixView(const c_type & rhs)915 inline DiagMatrixView(const c_type& rhs) : c_type(rhs) {} 916 DiagMatrixView(VectorView<T> _diag)917 explicit inline DiagMatrixView(VectorView<T> _diag) : 918 c_type(_diag) {} 919 ~DiagMatrixView()920 virtual inline ~DiagMatrixView() {} 921 922 // 923 // Op= 924 // 925 926 inline type& operator=(const type& m2) 927 { c_type::operator=(m2); return *this; } 928 929 inline type& operator=(const c_type& m2) 930 { c_type::operator=(m2); return *this; } 931 932 inline type& operator=(const GenDiagMatrix<RT>& m2) 933 { c_type::operator=(m2); return *this; } 934 935 inline type& operator=(const GenDiagMatrix<CT>& m2) 936 { c_type::operator=(m2); return *this; } 937 938 template <typename T2> 939 inline type& operator=(const GenDiagMatrix<T2>& m2) 940 { c_type::operator=(m2); return *this; } 941 942 inline type& operator=(const T& x) 943 { c_type::operator=(x); return *this; } 944 945 inline type& operator=(const AssignableToDiagMatrix<RT>& m2) 946 { c_type::operator=(m2); return *this; } 947 948 inline type& operator=(const AssignableToDiagMatrix<CT>& m2) 949 { c_type::operator=(m2); return *this; } 950 951 typedef typename c_type::MyListAssigner MyListAssigner; 952 inline MyListAssigner operator<<(const T& x) 953 { return c_type::operator<<(x); } 954 955 956 // 957 // Access 958 // 959 operator()960 inline reference operator()(ptrdiff_t i) 961 { 962 TMVAssert(i>0 && i<=c_type::size()); 963 return diag()(i); 964 } operator()965 inline reference operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j)) 966 { 967 TMVAssert(i==j); 968 TMVAssert(i>0 && i<=c_type::size()); 969 return diag()(i); 970 } 971 diag()972 inline vec_type diag() 973 { return c_type::diag(); } 974 operator()975 inline T operator()(ptrdiff_t i) const 976 { 977 TMVAssert(i>0 && i<=c_type::size()); 978 return diag()(i); 979 } operator()980 inline T operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j)) const 981 { 982 TMVAssert(i==j); 983 TMVAssert(i>0 && i<=c_type::size()); 984 return diag()(i); 985 } 986 diag()987 inline const_vec_type diag() const 988 { return c_type::diag(); } 989 990 // 991 // Modifying Functions 992 // 993 setZero()994 inline type& setZero() 995 { diag().setZero(); return *this; } 996 setAllTo(const T & x)997 inline type& setAllTo(const T& x) 998 { diag().setAllTo(x); return *this; } 999 addToAll(const T & x)1000 inline type& addToAll(const T& x) 1001 { diag().addToAll(x); return *this; } 1002 clip(RT thresh)1003 inline type& clip(RT thresh) 1004 { diag().clip(thresh); return *this; } 1005 transposeSelf()1006 inline type& transposeSelf() 1007 { return *this; } 1008 conjugateSelf()1009 inline type& conjugateSelf() 1010 { diag().conjugateSelf(); return *this; } 1011 invertSelf()1012 inline type& invertSelf() 1013 { return c_type::invertSelf(); } 1014 1015 inline type& setToIdentity(const T& x=T(1)) 1016 { diag().setAllTo(x); return *this; } 1017 1018 1019 // 1020 // subDiagMatrix 1021 // 1022 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1023 inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) 1024 { 1025 TMVAssert(diag().hasSubVector(i1,i2,1)); 1026 return c_type::cSubDiagMatrix(i1-1,i2); 1027 } 1028 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1029 inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 1030 { 1031 TMVAssert(diag().hasSubVector(i1,i2,istep)); 1032 return c_type::cSubDiagMatrix(i1-1,i2-1+istep,istep); 1033 } 1034 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1035 inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 1036 { 1037 TMVAssert(diag().hasSubVector(i1,i2,istep)); 1038 return view_type(diag().subVector(i1,i2,istep)); 1039 } 1040 realPart()1041 inline realpart_type realPart() 1042 { return realpart_type(diag().realPart()); } 1043 imagPart()1044 inline realpart_type imagPart() 1045 { return realpart_type(diag().imagPart()); } 1046 view()1047 inline view_type view() 1048 { return *this; } 1049 transpose()1050 inline view_type transpose() 1051 { return *this; } 1052 conjugate()1053 inline view_type conjugate() 1054 { return view_type(diag().conjugate()); } 1055 adjoint()1056 inline view_type adjoint() 1057 { return conjugate(); } 1058 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1059 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1060 { 1061 TMVAssert(diag().hasSubVector(i1,i2,1)); 1062 return c_type::cSubDiagMatrix(i1-1,i2); 1063 } 1064 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1065 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1066 { 1067 TMVAssert(diag().hasSubVector(i1,i2,istep)); 1068 return c_type::cSubDiagMatrix(i1-1,i2-1+istep,istep); 1069 } 1070 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1071 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1072 { 1073 TMVAssert(diag().hasSubVector(i1,i2,istep)); 1074 return view_type(diag().subVector(i1,i2,istep)); 1075 } 1076 realPart()1077 inline const_realpart_type realPart() const 1078 { return realpart_type(diag().realPart()); } 1079 imagPart()1080 inline const_realpart_type imagPart() const 1081 { return realpart_type(diag().imagPart()); } 1082 view()1083 inline const_view_type view() const 1084 { return *this; } 1085 transpose()1086 inline const_view_type transpose() const 1087 { return *this; } 1088 conjugate()1089 inline const_view_type conjugate() const 1090 { return view_type(diag().conjugate()); } 1091 adjoint()1092 inline const_view_type adjoint() const 1093 { return conjugate(); } 1094 1095 }; // FortranStyle DiagMatrixView 1096 1097 1098 template <typename T, int A> 1099 class DiagMatrix : public GenDiagMatrix<T> 1100 { 1101 public: 1102 1103 typedef TMV_RealType(T) RT; 1104 typedef TMV_ComplexType(T) CT; 1105 typedef GenDiagMatrix<T> base; 1106 typedef DiagMatrix<T,A> type; 1107 typedef DiagMatrixView<T,A> view_type; 1108 typedef view_type transpose_type; 1109 typedef view_type conjugate_type; 1110 typedef view_type adjoint_type; 1111 typedef DiagMatrixView<RT,A> realpart_type; 1112 typedef VectorView<T,A> vec_type; 1113 typedef ConstDiagMatrixView<T,A> const_view_type; 1114 typedef const_view_type const_transpose_type; 1115 typedef const_view_type const_conjugate_type; 1116 typedef const_view_type const_adjoint_type; 1117 typedef ConstDiagMatrixView<RT,A> const_realpart_type; 1118 typedef ConstVectorView<T,A> const_vec_type; 1119 typedef T& reference; 1120 typedef typename const_vec_type::const_iterator const_iterator; 1121 typedef typename vec_type::iterator iterator; 1122 1123 // 1124 // Constructors 1125 // 1126 itsdiag(_size)1127 inline explicit DiagMatrix(ptrdiff_t _size=0) : itsdiag(_size) 1128 { TMVAssert(_size >= 0); } 1129 DiagMatrix(ptrdiff_t _size,const T & x)1130 inline DiagMatrix(ptrdiff_t _size, const T& x) : itsdiag(_size,x) 1131 { TMVAssert(_size >= 0); } 1132 DiagMatrix(const GenVector<T> & rhs)1133 inline explicit DiagMatrix(const GenVector<T>& rhs) : itsdiag(rhs) {} 1134 DiagMatrix(const GenMatrix<T> & m)1135 inline explicit DiagMatrix(const GenMatrix<T>& m) : itsdiag(m.diag()) {} 1136 DiagMatrix(const type & rhs)1137 inline DiagMatrix(const type& rhs) : itsdiag(rhs.diag()) {} 1138 DiagMatrix(const GenDiagMatrix<RT> & rhs)1139 inline DiagMatrix(const GenDiagMatrix<RT>& rhs) : itsdiag(rhs.size()) 1140 { rhs.assignToD(view()); } 1141 DiagMatrix(const GenDiagMatrix<CT> & rhs)1142 inline DiagMatrix(const GenDiagMatrix<CT>& rhs) : itsdiag(rhs.size()) 1143 { 1144 TMVAssert(isComplex(T())); 1145 rhs.assignToD(view()); 1146 } 1147 1148 template <typename T2> DiagMatrix(const GenDiagMatrix<T2> & rhs)1149 inline DiagMatrix(const GenDiagMatrix<T2>& rhs) : itsdiag(rhs.diag()) {} 1150 DiagMatrix(const AssignableToDiagMatrix<RT> & m2)1151 inline DiagMatrix(const AssignableToDiagMatrix<RT>& m2) : 1152 itsdiag(m2.colsize()) 1153 { 1154 TMVAssert(m2.colsize() == m2.rowsize()); 1155 m2.assignToD(view()); 1156 } 1157 DiagMatrix(const AssignableToDiagMatrix<CT> & m2)1158 inline DiagMatrix(const AssignableToDiagMatrix<CT>& m2) : 1159 itsdiag(m2.colsize()) 1160 { 1161 TMVAssert(m2.colsize() == m2.rowsize()); 1162 TMVAssert(isComplex(T())); 1163 m2.assignToD(view()); 1164 } 1165 ~DiagMatrix()1166 virtual inline ~DiagMatrix() {} 1167 1168 1169 // 1170 // Op= 1171 // 1172 1173 inline type& operator=(const type& m2) 1174 { 1175 TMVAssert(m2.size() == size()); 1176 m2.assignToD(view()); 1177 return *this; 1178 } 1179 1180 inline type& operator=(const GenDiagMatrix<RT>& m2) 1181 { 1182 TMVAssert(m2.size() == size()); 1183 m2.assignToD(view()); 1184 return *this; 1185 } 1186 1187 inline type& operator=(const GenDiagMatrix<CT>& m2) 1188 { 1189 TMVAssert(m2.size() == size()); 1190 TMVAssert(isComplex(T())); 1191 m2.assignToD(view()); 1192 return *this; 1193 } 1194 1195 template <typename T2> 1196 inline type& operator=(const GenDiagMatrix<T2>& m2) 1197 { 1198 TMVAssert(m2.size() == size()); 1199 view() = m2; 1200 return *this; 1201 } 1202 1203 inline type& operator=(const T& x) { view() = x; return *this; } 1204 1205 inline type& operator=(const AssignableToDiagMatrix<RT>& m2) 1206 { 1207 TMVAssert(m2.size() == size()); 1208 m2.assignToD(view()); 1209 return *this; 1210 } 1211 1212 inline type& operator=(const AssignableToDiagMatrix<CT>& m2) 1213 { 1214 TMVAssert(m2.size() == size()); 1215 TMVAssert(isComplex(T())); 1216 m2.assignToD(view()); 1217 return *this; 1218 } 1219 1220 typedef ListAssigner<T,iterator> MyListAssigner; 1221 inline MyListAssigner operator<<(const T& x) 1222 { return MyListAssigner(begin(),size(),x); } 1223 1224 1225 // 1226 // Access 1227 // 1228 operator()1229 inline T& operator()(ptrdiff_t i) 1230 { 1231 if (A==CStyle) { 1232 TMVAssert(i>=0 && i<size()); 1233 return itsdiag(i); 1234 } else { 1235 TMVAssert(i>0 && i<=size()); 1236 return itsdiag(i-1); 1237 } 1238 } 1239 operator()1240 inline T& operator()(ptrdiff_t i, ptrdiff_t TMV_DEBUGPARAM(j)) 1241 { 1242 if (A==CStyle) { 1243 TMVAssert(i>=0 && i<size()); 1244 TMVAssert(j>=0 && j<size()); 1245 } else { 1246 TMVAssert(i>0 && i<=size()); 1247 TMVAssert(j>0 && j<=size()); 1248 } 1249 TMVAssert(i==j); 1250 return operator()(i); 1251 } 1252 operator()1253 inline T operator()(ptrdiff_t i) const 1254 { 1255 if (A==CStyle) { 1256 TMVAssert(i>=0 && i<size()); 1257 return itsdiag(i); 1258 } else { 1259 TMVAssert(i>0 && i<=size()); 1260 return itsdiag(i-1); 1261 } 1262 } 1263 operator()1264 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 1265 { 1266 if (A==CStyle) { 1267 TMVAssert(i>=0 && i<size()); 1268 TMVAssert(j>=0 && j<size()); 1269 } else { 1270 TMVAssert(i>0 && i<=size()); 1271 TMVAssert(j>0 && j<=size()); 1272 } 1273 if (i==j) return operator()(i); 1274 else return T(0); 1275 } 1276 diag()1277 inline vec_type diag() { return itsdiag.view(); } 1278 diag()1279 inline const_vec_type diag() const { return itsdiag.view(); } 1280 1281 // 1282 // Modifying Functions 1283 // 1284 setZero()1285 inline type& setZero() 1286 { itsdiag.setZero(); return *this; } 1287 setAllTo(const T & x)1288 inline type& setAllTo(const T& x) 1289 { itsdiag.setAllTo(x); return *this; } 1290 addToAll(const T & x)1291 inline type& addToAll(const T& x) 1292 { itsdiag.addToAll(x); return *this; } 1293 clip(RT thresh)1294 inline type& clip(RT thresh) 1295 { diag().clip(thresh); return *this; } 1296 transposeSelf()1297 inline type& transposeSelf() 1298 { return *this; } 1299 conjugateSelf()1300 inline type& conjugateSelf() 1301 { itsdiag.conjugateSelf(); return *this; } 1302 invertSelf()1303 inline type& invertSelf() 1304 { view().invertSelf(); return *this; } 1305 1306 inline type& setToIdentity(const T& x=T(1)) 1307 { itsdiag.setAllTo(x); return *this; } 1308 1309 // 1310 // subDiagMatrix 1311 // 1312 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1313 inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) 1314 { return view_type(itsdiag.cSubVector(i1,i2)); } 1315 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1316 inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) 1317 { 1318 TMVAssert(diag().hasSubVector(i1,i2,1)); 1319 if (A == int(FortranStyle)) { --i1; } 1320 return cSubDiagMatrix(i1,i2); 1321 } 1322 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1323 inline view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 1324 { return view_type(itsdiag.cSubVector(i1,i2,istep)); } 1325 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1326 inline view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 1327 { 1328 TMVAssert(diag().hasSubVector(i1,i2,istep)); 1329 if (A == int(FortranStyle)) { --i1; i2 += istep-1; } 1330 return cSubDiagMatrix(i1,i2,istep); 1331 } 1332 realPart()1333 inline realpart_type realPart() 1334 { return realpart_type(diag().realPart()); } 1335 imagPart()1336 inline realpart_type imagPart() 1337 { return realpart_type(diag().imagPart()); } 1338 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1339 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1340 { return const_view_type(itsdiag.cSubVector(i1,i2)); } 1341 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2)1342 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1343 { 1344 TMVAssert(diag().hasSubVector(i1,i2,1)); 1345 if (A == int(FortranStyle)) { --i1; } 1346 return cSubDiagMatrix(i1,i2); 1347 } 1348 cSubDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1349 inline const_view_type cSubDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1350 { return const_view_type(itsdiag.cSubVector(i1,i2,istep)); } 1351 subDiagMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep)1352 inline const_view_type subDiagMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1353 { 1354 TMVAssert(diag().hasSubVector(i1,i2,istep)); 1355 if (A == int(FortranStyle)) { --i1; i2 += istep-1; } 1356 return cSubDiagMatrix(i1,i2,istep); 1357 } 1358 realPart()1359 inline const_realpart_type realPart() const 1360 { return const_realpart_type(diag().realPart()); } 1361 imagPart()1362 inline const_realpart_type imagPart() const 1363 { return const_realpart_type(diag().imagPart()); } 1364 view()1365 inline const_view_type view() const 1366 { return const_view_type(itsdiag); } 1367 transpose()1368 inline const_view_type transpose() const 1369 { return view(); } 1370 conjugate()1371 inline const_view_type conjugate() const 1372 { return const_view_type(itsdiag.conjugate()); } 1373 adjoint()1374 inline const_view_type adjoint() const 1375 { return conjugate(); } 1376 view()1377 inline view_type view() 1378 { return view_type(itsdiag.view()); } 1379 transpose()1380 inline view_type transpose() 1381 { return view(); } 1382 conjugate()1383 inline view_type conjugate() 1384 { return view_type(itsdiag.conjugate()); } 1385 adjoint()1386 inline view_type adjoint() 1387 { return conjugate(); } 1388 1389 // 1390 // I/O 1391 // 1392 1393 void read(const TMV_Reader& reader); 1394 1395 using base::size; 1396 resize(ptrdiff_t n)1397 inline void resize(ptrdiff_t n) 1398 { 1399 TMVAssert(n >= 0); 1400 itsdiag.resize(n); 1401 } 1402 begin()1403 inline iterator begin() 1404 { return diag().begin(); } end()1405 inline iterator end() 1406 { return diag().end(); } 1407 begin()1408 inline const_iterator begin() const 1409 { return diag().begin(); } end()1410 inline const_iterator end() const 1411 { return diag().end(); } 1412 1413 protected : 1414 1415 Vector<T> itsdiag; cdiag()1416 inline ConstVectorView<T> cdiag() const { return itsdiag.view(); } 1417 Swap(DiagMatrix<T,A> & m1,DiagMatrix<T,A> & m2)1418 friend void Swap(DiagMatrix<T,A>& m1, DiagMatrix<T,A>& m2) 1419 { 1420 TMVAssert(m1.size() == m2.size()); 1421 Swap(m1.itsdiag,m2.itsdiag); 1422 } 1423 1424 }; // DiagMatrix 1425 1426 //--------------------------------------------------------------------------- 1427 1428 // 1429 // Special Creators: 1430 // DiagMatrixViewOf(v) 1431 // DiagMatrixViewOf(T* v, n) 1432 // DiagMatrixViewOf(T* v, n, step) 1433 // 1434 1435 template <typename T> DiagMatrixViewOf(const GenVector<T> & v)1436 inline ConstDiagMatrixView<T> DiagMatrixViewOf(const GenVector<T>& v) 1437 { return ConstDiagMatrixView<T>(v); } 1438 1439 template <typename T, int A> DiagMatrixViewOf(const ConstVectorView<T,A> & v)1440 inline ConstDiagMatrixView<T,A> DiagMatrixViewOf( 1441 const ConstVectorView<T,A>& v) 1442 { return ConstDiagMatrixView<T,A>(v); } 1443 1444 template <typename T, int A> DiagMatrixViewOf(const Vector<T,A> & v)1445 inline ConstDiagMatrixView<T,A> DiagMatrixViewOf(const Vector<T,A>& v) 1446 { return ConstDiagMatrixView<T,A>(v.view()); } 1447 1448 template <typename T, int A> DiagMatrixViewOf(VectorView<T,A> v)1449 inline DiagMatrixView<T,A> DiagMatrixViewOf(VectorView<T,A> v) 1450 { return DiagMatrixView<T,A>(v); } 1451 1452 template <typename T, int A> DiagMatrixViewOf(Vector<T,A> & v)1453 inline DiagMatrixView<T,A> DiagMatrixViewOf(Vector<T,A>& v) 1454 { return DiagMatrixView<T,A>(v.view()); } 1455 1456 template <typename T> DiagMatrixViewOf(const T * m,ptrdiff_t size)1457 inline ConstDiagMatrixView<T> DiagMatrixViewOf(const T* m, ptrdiff_t size) 1458 { 1459 TMVAssert(size >= 0); 1460 return ConstDiagMatrixView<T>(VectorViewOf(m,size)); 1461 } 1462 1463 template <typename T> DiagMatrixViewOf(T * m,ptrdiff_t size)1464 inline DiagMatrixView<T> DiagMatrixViewOf(T* m, ptrdiff_t size) 1465 { 1466 TMVAssert(size >= 0); 1467 return DiagMatrixView<T>(VectorViewOf(m,size)); 1468 } 1469 1470 template <typename T> DiagMatrixViewOf(const T * m,ptrdiff_t size,ptrdiff_t step)1471 inline ConstDiagMatrixView<T> DiagMatrixViewOf( 1472 const T* m, ptrdiff_t size, ptrdiff_t step) 1473 { 1474 TMVAssert(size >= 0); 1475 return ConstDiagMatrixView<T>(VectorViewOf(m,size,step)); 1476 } 1477 1478 template <typename T> DiagMatrixViewOf(T * m,ptrdiff_t size,ptrdiff_t step)1479 inline DiagMatrixView<T> DiagMatrixViewOf(T* m, ptrdiff_t size, ptrdiff_t step) 1480 { 1481 TMVAssert(size >= 0); 1482 return DiagMatrixView<T>(VectorViewOf(m,size,step)); 1483 } 1484 1485 1486 // 1487 // Swap Matrices 1488 // 1489 1490 template <typename T> Swap(DiagMatrixView<T> m1,DiagMatrixView<T> m2)1491 inline void Swap(DiagMatrixView<T> m1, DiagMatrixView<T> m2) 1492 { Swap(m1.diag(),m2.diag()); } 1493 1494 template <typename T, int A> Swap(DiagMatrix<T,A> & m1,DiagMatrixView<T> m2)1495 inline void Swap(DiagMatrix<T,A>& m1, DiagMatrixView<T> m2) 1496 { Swap(m1.diag(),m2.diag()); } 1497 1498 template <typename T, int A> Swap(DiagMatrixView<T> m1,DiagMatrix<T,A> & m2)1499 inline void Swap(DiagMatrixView<T> m1, DiagMatrix<T,A>& m2) 1500 { Swap(m1.diag(),m2.diag()); } 1501 1502 // 1503 // Views: 1504 // 1505 1506 template <typename T> Transpose(const GenDiagMatrix<T> & m)1507 inline ConstDiagMatrixView<T> Transpose(const GenDiagMatrix<T>& m) 1508 { return m.transpose(); } 1509 1510 template <typename T, int A> Transpose(const ConstDiagMatrixView<T,A> & m)1511 inline ConstDiagMatrixView<T,A> Transpose(const ConstDiagMatrixView<T,A>& m) 1512 { return m.transpose(); } 1513 1514 template <typename T, int A> Transpose(const DiagMatrix<T,A> & m)1515 inline ConstDiagMatrixView<T,A> Transpose(const DiagMatrix<T,A>& m) 1516 { return m.transpose(); } 1517 1518 template <typename T, int A> Transpose(DiagMatrixView<T,A> m)1519 inline DiagMatrixView<T,A> Transpose(DiagMatrixView<T,A> m) 1520 { return m.transpose(); } 1521 1522 template <typename T, int A> Transpose(DiagMatrix<T,A> & m)1523 inline DiagMatrixView<T,A> Transpose(DiagMatrix<T,A>& m) 1524 { return m.transpose(); } 1525 1526 template <typename T> Conjugate(const GenDiagMatrix<T> & m)1527 inline ConstDiagMatrixView<T> Conjugate(const GenDiagMatrix<T>& m) 1528 { return m.conjugate(); } 1529 1530 template <typename T, int A> Conjugate(const ConstDiagMatrixView<T,A> & m)1531 inline ConstDiagMatrixView<T,A> Conjugate(const ConstDiagMatrixView<T,A>& m) 1532 { return m.conjugate(); } 1533 1534 template <typename T, int A> Conjugate(const DiagMatrix<T,A> & m)1535 inline ConstDiagMatrixView<T,A> Conjugate(const DiagMatrix<T,A>& m) 1536 { return m.conjugate(); } 1537 1538 template <typename T, int A> Conjugate(DiagMatrixView<T,A> m)1539 inline DiagMatrixView<T,A> Conjugate(DiagMatrixView<T,A> m) 1540 { return m.conjugate(); } 1541 1542 template <typename T, int A> Conjugate(DiagMatrix<T,A> & m)1543 inline DiagMatrixView<T,A> Conjugate(DiagMatrix<T,A>& m) 1544 { return m.conjugate(); } 1545 1546 template <typename T> Adjoint(const GenDiagMatrix<T> & m)1547 inline ConstDiagMatrixView<T> Adjoint(const GenDiagMatrix<T>& m) 1548 { return m.adjoint(); } 1549 1550 template <typename T, int A> Adjoint(const ConstDiagMatrixView<T,A> & m)1551 inline ConstDiagMatrixView<T,A> Adjoint(const ConstDiagMatrixView<T,A>& m) 1552 { return m.adjoint(); } 1553 1554 template <typename T, int A> Adjoint(const DiagMatrix<T,A> & m)1555 inline ConstDiagMatrixView<T,A> Adjoint(const DiagMatrix<T,A>& m) 1556 { return m.adjoint(); } 1557 1558 template <typename T, int A> Adjoint(DiagMatrixView<T,A> m)1559 inline DiagMatrixView<T,A> Adjoint(DiagMatrixView<T,A> m) 1560 { return m.adjoint(); } 1561 1562 template <typename T, int A> Adjoint(DiagMatrix<T,A> & m)1563 inline DiagMatrixView<T,A> Adjoint(DiagMatrix<T,A>& m) 1564 { return m.adjoint(); } 1565 1566 template <typename T> Inverse(const GenDiagMatrix<T> & m)1567 inline QuotXD<T,T> Inverse(const GenDiagMatrix<T>& m) 1568 { return m.inverse(); } 1569 1570 1571 1572 // 1573 // DiagMatrix ==, != DiagMatrix 1574 // 1575 1576 template <typename T1, typename T2> 1577 inline bool operator==( 1578 const GenDiagMatrix<T1>& m1, const GenDiagMatrix<T2>& m2) 1579 { return m1.diag() == m2.diag(); } 1580 1581 template <typename T1, typename T2> 1582 inline bool operator!=( 1583 const GenDiagMatrix<T1>& m1, const GenDiagMatrix<T2>& m2) 1584 { return !(m1 == m2); } 1585 1586 template <typename T1, typename T2> 1587 inline bool operator==( 1588 const GenDiagMatrix<T1>& m1, const GenMatrix<T2>& m2) 1589 { 1590 return 1591 m1.diag() == m2.diag() && 1592 m2.upperTri().offDiag().maxAbs2Element() == T2(0) && 1593 m2.lowerTri().offDiag().maxAbs2Element() == T2(0); 1594 } 1595 1596 template <typename T1, typename T2> 1597 inline bool operator==( 1598 const GenMatrix<T1>& m1, const GenDiagMatrix<T2>& m2) 1599 { return m2 == m1; } 1600 1601 template <typename T1, typename T2> 1602 inline bool operator!=( 1603 const GenDiagMatrix<T1>& m1, const GenMatrix<T2>& m2) 1604 { return !(m1 == m2); } 1605 1606 template <typename T1, typename T2> 1607 inline bool operator!=( 1608 const GenMatrix<T1>& m1, const GenDiagMatrix<T2>& m2) 1609 { return !(m1 == m2); } 1610 1611 // 1612 // I/O 1613 // 1614 1615 template <typename T> 1616 inline std::istream& operator>>( 1617 const TMV_Reader& reader, DiagMatrixView<T> m) 1618 { m.read(reader); return reader.getis(); } 1619 1620 template <typename T, int A> 1621 inline std::istream& operator>>( 1622 const TMV_Reader& reader, DiagMatrix<T,A>& m) 1623 { m.read(reader); return reader.getis(); } 1624 1625 template <typename T> 1626 inline std::istream& operator>>( 1627 std::istream& is, DiagMatrixView<T> m) 1628 { return is >> IOStyle() >> m; } 1629 1630 template <typename T, int A> 1631 inline std::istream& operator>>(std::istream& is, DiagMatrix<T,A>& m) 1632 { return is >> IOStyle() >> m; } 1633 1634 } // namespace tmv 1635 1636 #endif 1637