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 TriMatrix class. 26 // 27 // There are two TriMatrix classes: UpperTriMatrix<T> and LowerTriMatrix<T> 28 // For these notes, I will just write TriMatrix, but for all uses, 29 // you need to write Upper or Lower before the Tri. 30 // 31 // As usual, the first template parameter is the type of the data, 32 // and the optional second template parameter specifies the known 33 // attributes. The valid attributs for a SymMatrix are: 34 // - ColMajor or RoMajor 35 // - CStyle or FortranStyle 36 // - NonUnitDiag or UnitDiag 37 // The defaults are (ColMajor,CStyle,NonUnitDiag) if you do not specify 38 // otherwise. 39 // 40 // The storage and index style follow the same meaning as for regular 41 // Matrices. 42 // NonUnitDiag or UnitDiag indicates whether or not the diagonal elements 43 // should be taken to be all 1's, or whether to use the actual values 44 // in memory. 45 // 46 // Constructors: 47 // 48 // TriMatrix<T,A>(int n) 49 // Makes a Triangular Matrix with column size = row size = n 50 // with _uninitialized_ values. 51 // 52 // TriMatrix<T,A>(int n, T x) 53 // Makes a Triangular Matrix with column size = row size = n 54 // with all values = x 55 // 56 // TriMatrix<T,A>(const Matrix<T>& m) 57 // TriMatrix<T,A>(const TriMatrix<T>& m) 58 // Makes a TriMatrix which copies the corresponding elements of m. 59 // 60 // ConstUpperTriMatrixView<T> UpperTriMatrixViewOf( 61 // const T* m, size, stor) 62 // UpperTriMatrixView<T> UpperTriMatrixViewOf(T* m, size, stor) 63 // ConstLowerTriMatrixView<T> LowerTriMatrixViewOf( 64 // const T* m, size, stor) 65 // LowerTriMatrixView<T> LowerTriMatrixViewOf(T* m, size, stor) 66 // Make a TriMatrix view of the specific location in memory 67 // specified by m. The size and storage are parameters. 68 // 69 // ConstUpperTriMatrixView<T> UnitUpperTriMatrixViewOf( 70 // const T* m, size, stor) 71 // UpperTriMatrixView<T> UnitUpperTriMatrixViewOf(T* m, size, stor) 72 // ConstLowerTriMatrixView<T> UnitLowerTriMatrixViewOf( 73 // const T* m, size, stor) 74 // LowerTriMatrixView<T> UnitLowerTriMatrixViewOf(T* m, size, stor) 75 // Same as above, but with unit-diagonal. 76 // 77 // Access Functions 78 // 79 // int colsize() const 80 // int rowsize() const 81 // int size() const 82 // Return the dimensions of the TriMatrix 83 // 84 // T& operator()(int i, int j) 85 // T operator()(int i, int j) const 86 // Return the (i,j) element of the TriMatrix 87 // 88 // Vector& row(int i, int j1, int j2) 89 // Return a portion of the ith row 90 // This range must be a valid range for the requested row. 91 // 92 // Vector& col(int j, int i1, int i2) 93 // Return a portion of the jth column 94 // This range must be a valid range for the requested column. 95 // 96 // Vector& diag() 97 // Return the main diagonal 98 // The TriMatrix must be NonUnitDiag. 99 // 100 // Vector& diag(int i, int j1, int j2) 101 // Vector& diag(int i) 102 // Return the super- or sub-diagonal i 103 // If i > 0 return the super diagonal starting at m_0i 104 // If i < 0 return the sub diagonal starting at m_|i|0 105 // If j1,j2 are given, it returns the diagonal subVector 106 // either from m_j1,i+j1 to m_j2,i+j2 (for i>0) 107 // or from m_|i|+j1,j1 to m_|i|+j2,j2 (for i<0) 108 // i>0 will give an error for a LowerTriMatrix 109 // i<0 will give an error for an UpperTriMatrix 110 // i=0 will give an error for a UnitDiag TriMatrix 111 // 112 // Modifying Functions: 113 // 114 // setZero() 115 // setAllTo(T x) 116 // addToAll(T x) 117 // conjugateSelf() 118 // setToIdentity(x = 1) 119 // void Swap(TriMatrix& m1, TriMatrix& m2) 120 // The TriMatrices must be the same size and shape (Upper or Lower). 121 // 122 // Views of a TriMatrix: 123 // 124 // subMatrix(int i1, int i2, int j1, int j2, int istep=1, int jstep=1) 125 // This member function will return a submatrix using rows i1 to i2 126 // and columns j1 to j2 which refers 127 // to the same physical elements as the original. 128 // The submatrix must be completely contained within the TriMatrix. 129 // 130 // subVector(int i, int j, int istep, int jstep, int size) 131 // Returns a subVector which starts at position (i,j) in the 132 // matrix, moves in the directions (istep,jstep) and has a length 133 // of size. 134 // 135 // subTriMatrix(int i1, int i2, int istep) 136 // Returns the TriMatrix which runs from i1 to i2 along the diagonal 137 // (not including i2) with an optional step, and includes the 138 // off diagonal in the same rows/cols. 139 // 140 // For example, with an UpperTriMatrix of size 10, the x's below 141 // are the original data, the O's are the subTriMatrix returned 142 // with the command subTriMatrix(3,11,2), and the #'s are the 143 // subTriMatrix returned with subTriMatrix(0,3) 144 // 145 // ###xxxxxxx 146 // ##xxxxxxx 147 // #xxxxxxx 148 // OxOxOxO 149 // xxxxxx 150 // OxOxO 151 // xxxx 152 // OxO 153 // xx 154 // O 155 // 156 // offDiag() 157 // Returns the (NonUnitDiag) TriMatrix of all the off-diagonal 158 // elements of a TriMatrix. 159 // 160 // realPart(), imagPart() 161 // For a complex TriMatrix, returns the real or imaginary part 162 // as a real TriMatrix. 163 // 164 // view() 165 // transpose() 166 // adjoint() 167 // Note that the transpose or adjoint of an UpperTriMatrix returns 168 // a view which is a LowerTriMatrix, and vice versa. 169 // conjugate() 170 // 171 // viewAsUnitDiag() 172 // Returns a UnitDiag view of a TriMatrix. 173 // i.e. a TriMatrix of the same elements, but treating the diagonl 174 // as all 1's. 175 // 176 // 177 // Functions of Matrixs: 178 // 179 // m.det() or Det(m) 180 // m.logDet() or m.logDet(T* sign) or LogDet(m) 181 // m.trace() or Trace(m) 182 // m.sumElements() or SumElements(m) 183 // m.sumAbsElements() or SumAbsElements(m) 184 // m.sumAbs2Elements() or SumAbs2Elements(m) 185 // m.norm() or m.normF() or Norm(m) or NormF(m) 186 // m.normSq() or NormSq(m) 187 // m.norm1() or Norm1(m) 188 // m.norm2() or Norm2(m) 189 // m.normInf() or NormInf(m) 190 // m.maxAbsElement() or MaxAbsElements(m) 191 // m.maxAbs2Element() or MaxAbs2Elements(m) 192 // 193 // m.inverse() or Inverse(m) 194 // m.makeInverse(minv) // Takes either a TriMatrix or Matrix argument 195 // m.makeInverseATA(invata) 196 // 197 // 198 // I/O: 199 // 200 // os << m 201 // Writes m to ostream os in the usual Matrix format 202 // 203 // os << CompactIO() << m 204 // Writes m to ostream os in the following compact format: 205 // For an UpperTriMatrix: 206 // size 207 // ( m(0,0) m(0,1) ... m(0,size) ) 208 // ( m(1,1) .. m(1,size) ) 209 // ... 210 // ( m(size,size) ) 211 // 212 // For a LowerTriMatrix: 213 // size 214 // ( m(0,0) ) 215 // ( m(1,0) m(1,1) ) 216 // ... 217 // ( m(size,0) ... m(size,size) ) 218 // 219 // is >> m 220 // is >> CompactIO() >> m 221 // reads m from istream in either format 222 // (Note: if the DiagType for the TriMatrix is UnitDiag, then 223 // all of the diagonals read in must be = 1.) 224 // 225 // 226 // Division Control Functions: 227 // 228 // Most of the point of using TriMatrixes is that they are easy 229 // to divide using either forward substitution or back substitution. 230 // Therefore, the only division available for TriMatrixes is 231 // this variety. To do something else (like SVD), you need to 232 // copy it to a regular matrix. 233 // 234 235 236 #ifndef TMV_TriMatrix_H 237 #define TMV_TriMatrix_H 238 239 #include "tmv/TMV_BaseTriMatrix.h" 240 #include "tmv/TMV_BaseDiagMatrix.h" 241 #include "tmv/TMV_Vector.h" 242 #include "tmv/TMV_Matrix.h" 243 #include "tmv/TMV_DiagMatrix.h" 244 #include "tmv/TMV_Array.h" 245 #include <vector> 246 247 namespace tmv { 248 249 template <typename T, bool C> 250 struct TriRefHelper // C = false 251 { 252 typedef T& reference; TriRefHelperTriRefHelper253 TriRefHelper(T& _r, bool ) : itsref(_r) {} TriRefHelperTriRefHelper254 TriRefHelper(const TriRefHelper& rhs) : itsref(rhs.itsref) {} 255 reference itsref; 256 }; 257 258 template <typename T> 259 struct TriRefHelper<T,true> // T is real 260 { 261 typedef T& reference; 262 TriRefHelper(T& _r, bool ) : itsref(_r) {} 263 TriRefHelper(const TriRefHelper& rhs) : itsref(rhs.itsref) {} 264 reference itsref; 265 }; 266 267 template <typename T> 268 struct TriRefHelper<std::complex<T>,true> // complex and maybe conj 269 { 270 typedef VarConjRef<std::complex<T> > reference; 271 TriRefHelper(std::complex<T>& _r, bool _c) : 272 itsref(_r,_c?Conj:NonConj) {} 273 TriRefHelper(const TriRefHelper& rhs) : itsref(rhs.itsref) {} 274 reference itsref; 275 }; 276 277 // A helper struct to provide a valid reference for TriMatrix, correctly 278 // dealing with the possibility of UnitDiag. 279 // It is very similar to VarConjRef, but also has a boolean isunit 280 // which indicates whether the reference is on the diagonal of a 281 // UnitDiag TriMatrix. If so, it can be an rhs with the value 1, 282 // but not a lhs. 283 // If C is true, then the underlying reference is a VarCojRef. 284 // Otherwise it is a simple T&. 285 template <typename T, bool C> 286 class TriRef 287 { 288 public: 289 typedef typename TriRefHelper<T,C>::reference reference; 290 explicit inline TriRef(bool _u, T& _v, bool _c) : 291 isunit(_u), helper(_v,_c) {} 292 inline TriRef(const TriRef<T,C>& rhs) : 293 isunit(rhs.isunit), helper(rhs.helper) {} 294 inline ~TriRef() {} 295 296 inline operator T() const { return val(); } 297 inline reference getRef() { return ref(); } 298 inline T operator-() const { return -val(); } 299 300 inline TriRef<T,C>& operator=(const TriRef<T,C>& rhs) 301 { assign(rhs.val()); return *this; } 302 template <typename T2> 303 inline TriRef<T,C>& operator=(const TriRef<T2,C>& rhs) 304 { assign(rhs.val()); return *this; } 305 template <typename T2> 306 inline TriRef<T,C>& operator=(T2 rhs) 307 { assign(rhs); return *this; } 308 309 template <typename T2> 310 inline TriRef<T,C>& operator+=(const TriRef<T2,C>& x2) 311 { assign(val() + x2.val()); return *this; } 312 template <typename T2> 313 inline TriRef<T,C>& operator+=(T2 x2) 314 { assign(val() + x2); return *this; } 315 template <typename T2> 316 inline typename Traits2<T,T2>::type operator+(const TriRef<T2,C>& x2) 317 { return val() + x2.val(); } 318 template <typename T2> 319 inline friend typename Traits2<T,T2>::type operator+( 320 const TriRef<T,C>& x1, T2 x2) 321 { return x1.val()+x2; } 322 template <typename T2> 323 inline friend typename Traits2<T,T2>::type operator+( 324 T2 x1, const TriRef<T,C>& x2) 325 { return x1+x2.val(); } 326 template <typename T2> 327 inline friend T2& operator+=(T2& x1, const TriRef<T,C>& x2) 328 { return x1 += x2.val(); } 329 330 template <typename T2> 331 inline TriRef<T,C>& operator-=(const TriRef<T2,C>& x2) 332 { assign(val() - x2.val()); return *this; } 333 template <typename T2> 334 inline TriRef<T,C>& operator-=(T2 x2) 335 { assign(val() - x2); return *this; } 336 template <typename T2> 337 inline typename Traits2<T,T2>::type operator-(const TriRef<T2,C>& x2) 338 { return val()-x2.val(); } 339 template <typename T2> 340 inline friend typename Traits2<T,T2>::type operator-( 341 const TriRef<T,C>& x1, T2 x2) 342 { return x1.val()-x2; } 343 template <typename T2> 344 inline friend typename Traits2<T,T2>::type operator-( 345 T2 x1, const TriRef<T,C>& x2) 346 { return x1-x2.val(); } 347 template <typename T2> 348 inline friend T2& operator-=(T2& x1, const TriRef<T,C>& x2) 349 { return x1 -= x2.val(); } 350 351 template <typename T2> 352 inline TriRef<T,C>& operator*=(const TriRef<T2,C>& x2) 353 { assign(x2.val() * val()); return *this; } 354 template <typename T2> 355 inline TriRef<T,C>& operator*=(T2 x2) 356 { assign(x2 * val()); return *this; } 357 template <typename T2> 358 inline typename Traits2<T,T2>::type operator*(const TriRef<T2,C> x2) 359 { return val()*x2.val(); } 360 template <typename T2> 361 inline friend typename Traits2<T,T2>::type operator*( 362 const TriRef<T,C>& x1, T2 x2) 363 { return x1.val()*x2; } 364 template <typename T2> 365 inline friend typename Traits2<T,T2>::type operator*( 366 T2 x1, const TriRef<T,C>& x2) 367 { return x1*x2.val(); } 368 template <typename T2> 369 inline friend T2& operator*=(T2& x1, const TriRef<T,C>& x2) 370 { 371 if (x2.isunit) return x1; 372 else return x1 *= x2.val(); 373 } 374 375 template <typename T2> 376 inline TriRef<T,C>& operator/=(const TriRef<T2,C>& x2) 377 { assign(val() / x2.val()); return *this; } 378 template <typename T2> 379 inline TriRef<T,C>& operator/=(T2 x2) 380 { assign(val() / x2); return *this; } 381 template <typename T2> 382 inline typename Traits2<T,T2>::type operator/(const TriRef<T2,C>& x2) 383 { return val()/x2.val(); } 384 template <typename T2> 385 inline friend typename Traits2<T,T2>::type operator/( 386 const TriRef<T,C>& x1, T2 x2) 387 { return x1.val()/x2; } 388 template <typename T2> 389 inline friend typename Traits2<T,T2>::type operator/( 390 T2 x1, const TriRef<T,C>& x2) 391 { return x1/x2.val(); } 392 template <typename T2> 393 inline friend T2& operator/=(T2& x1, const TriRef<T,C>& x2) 394 { 395 if (x2.isunit) return x1; 396 else return x1 /= x2.val(); 397 } 398 399 template <typename T2> 400 inline bool operator==(const TriRef<T2,C>& x2) const 401 { return val() == x2.val(); } 402 template <typename T2> 403 inline bool operator==(T2 x2) const 404 { return val() == x2; } 405 template <typename T2> 406 inline friend bool operator==(T2 x1, const TriRef<T,C>& x2) 407 { return x1 == x2.val(); } 408 template <typename T2> 409 inline bool operator!=(const TriRef<T2,C>& x2) const 410 { return !(operator==(x2)); } 411 template <typename T2> 412 inline bool operator!=(T2 x2) const 413 { return !(operator==(x2)); } 414 template <typename T2> 415 inline friend bool operator!=(T2 x1, const TriRef<T,C>& x2) 416 { return !(x2==x1); } 417 418 inline friend std::ostream& operator<<(std::ostream& os, TriRef<T,C> x) 419 { os << x.val(); return os; } 420 inline friend std::istream& operator>>(std::istream& is, TriRef<T,C> x) 421 { is >> x.ref(); return is; } 422 423 private: 424 425 reference ref() 426 { 427 TMVAssert(!isunit); 428 return helper.itsref; 429 } 430 template <typename T2> 431 void assign(T2 x) 432 { 433 TMVAssert(!isunit || x == T2(1)); 434 if (!isunit) helper.itsref = x; 435 } 436 T val() const 437 { return isunit ? T(1) : T(helper.itsref); } 438 439 const bool isunit; 440 TriRefHelper<T,C> helper; 441 }; 442 443 // Overload some functions to work with TriRef 444 template <typename T, bool C> 445 inline T TMV_CONJ(const TriRef<T,C>& x) 446 { return TMV_CONJ(T(x)); } 447 template <typename T, bool C> 448 inline typename Traits<T>::real_type TMV_NORM(const TriRef<T,C>& x) 449 { return TMV_NORM(T(x)); } 450 template <typename T, bool C> 451 inline typename Traits<T>::real_type TMV_ABS(const TriRef<T,C>& x) 452 { return TMV_ABS(T(x)); } 453 template <typename T, bool C> 454 inline T TMV_SQR(const TriRef<T,C>& x) 455 { return TMV_SQR(T(x)); } 456 template <typename T, bool C> 457 inline T TMV_SQRT(const TriRef<T,C>& x) 458 { return TMV_SQRT(T(x)); } 459 template <typename T, bool C> 460 inline typename Traits<T>::real_type TMV_REAL(const TriRef<T,C>& x) 461 { return TMV_REAL(T(x)); } 462 template <typename T, bool C> 463 inline typename Traits<T>::real_type TMV_IMAG(const TriRef<T,C>& x) 464 { return TMV_IMAG(T(x)); } 465 466 template <typename T, bool C1, bool C2> 467 inline void TMV_SWAP(TriRef<T,C1> x1, TriRef<T,C2> x2) 468 { TMV_SWAP(x1.getRef(),x2.getRef()); } 469 template <typename T, bool C2> 470 inline void TMV_SWAP(T& x1, TriRef<T,C2> x2) 471 { TMV_SWAP(x1,x2.getRef()); } 472 template <typename T, bool C1> 473 inline void TMV_SWAP(TriRef<T,C1> x1, T& x2) 474 { TMV_SWAP(x1.getRef(),x2); } 475 476 // Another helper class to deal with the case of the regular 477 // UpperTriMatrix and LowerTriMatrix classes (i.e. not views) 478 // Here, we can use a simple T& for the reference if D is NonUnitDiag 479 // Also, here C is always false. 480 template <typename T, int A> 481 struct TriRefHelper2 // D = NonUnitDiag 482 { 483 typedef T& reference; 484 static reference makeRef(bool, T& r) { return r; } 485 }; 486 template <typename T> 487 struct TriRefHelper2<T,UnitDiag> 488 { 489 typedef TriRef<T,false> reference; 490 static reference makeRef(bool ondiag, T& r) 491 { return reference(ondiag,r,false); } 492 }; 493 494 template <typename T> 495 class GenUpperTriMatrix : 496 virtual public AssignableToUpperTriMatrix<T>, 497 public BaseMatrix<T> 498 { 499 public: 500 501 typedef TMV_RealType(T) RT; 502 typedef TMV_ComplexType(T) CT; 503 typedef T value_type; 504 typedef RT real_type; 505 typedef CT complex_type; 506 typedef GenUpperTriMatrix<T> type; 507 typedef UpperTriMatrix<T> copy_type; 508 typedef ConstVectorView<T> const_vec_type; 509 typedef ConstMatrixView<T> const_rec_type; 510 typedef ConstUpperTriMatrixView<T> const_uppertri_type; 511 typedef ConstLowerTriMatrixView<T> const_lowertri_type; 512 typedef const_uppertri_type const_view_type; 513 typedef const_lowertri_type const_transpose_type; 514 typedef const_uppertri_type const_conjugate_type; 515 typedef const_lowertri_type const_adjoint_type; 516 typedef ConstUpperTriMatrixView<RT> const_realpart_type; 517 typedef UpperTriMatrixView<T> nonconst_type; 518 typedef CRMIt<type> const_rowmajor_iterator; 519 typedef CCMIt<type> const_colmajor_iterator; 520 521 // 522 // Constructors 523 // 524 525 inline GenUpperTriMatrix() {} 526 inline GenUpperTriMatrix(const type&) {} 527 virtual inline ~GenUpperTriMatrix() {} 528 529 // 530 // Access Functions 531 // 532 533 using AssignableToUpperTriMatrix<T>::size; 534 using AssignableToUpperTriMatrix<T>::dt; 535 inline ptrdiff_t colsize() const { return size(); } 536 inline ptrdiff_t rowsize() const { return size(); } 537 538 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 539 { 540 TMVAssert(i>=0 && i<size()); 541 TMVAssert(j>=0 && j<size()); 542 return cref(i,j); 543 } 544 545 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 546 { 547 TMVAssert(i>=0 && i<size()); 548 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 549 TMVAssert(j1==j2 || okij(i,j1)); 550 return const_vec_type( 551 cptr()+i*stepi()+j1*stepj(),j2-j1,stepj(),ct()); 552 } 553 554 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 555 { 556 TMVAssert(j>=0 && j<size()); 557 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 558 TMVAssert(i1==i2 || okij(i2-1,j)); 559 return const_vec_type( 560 cptr()+i1*stepi()+j*stepj(),i2-i1,stepi(),ct()); 561 } 562 563 inline const_vec_type diag() const 564 { 565 TMVAssert(!isunit()); 566 return const_vec_type(cptr(),size(),stepi()+stepj(),ct()); 567 } 568 569 inline const_vec_type diag(ptrdiff_t i) const 570 { 571 TMVAssert(isunit() ? i>0 : i>=0); 572 TMVAssert(i<=size()); 573 return const_vec_type( 574 cptr()+i*stepj(),size()-i,stepi()+stepj(),ct()); 575 } 576 577 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 578 { 579 TMVAssert(isunit() ? i>0 : i>=0); 580 TMVAssert(i<=size()); 581 TMVAssert(j1>=0 && j1 <= j2 && j2 <= size()-i); 582 const ptrdiff_t ds = stepi()+stepj(); 583 return const_vec_type(cptr()+i*stepj()+j1*ds,j2-j1,ds,ct()); 584 } 585 586 template <typename T2> 587 inline bool isSameAs(const BaseMatrix<T2>& ) const 588 { return false; } 589 590 inline bool isSameAs(const GenUpperTriMatrix<T>& m2) const 591 { 592 return 593 this == &m2 || 594 ( cptr()==m2.cptr() && size()==m2.size() && 595 dt() == m2.dt() && ct() == m2.ct() && 596 stepi()==m2.stepi() && stepj()==m2.stepj() 597 ); 598 } 599 600 inline void assignToM(MatrixView<RT> m2) const 601 { 602 TMVAssert(m2.colsize() == size()); 603 TMVAssert(m2.rowsize() == size()); 604 TMVAssert(isReal(T())); 605 assignToU(m2.upperTri(dt())); 606 if (isunit()) m2.diag().setAllTo(RT(1)); 607 if (size() > 0) m2.lowerTri().offDiag().setZero(); 608 } 609 610 inline void assignToM(MatrixView<CT> m2) const 611 { 612 TMVAssert(m2.colsize() == size()); 613 TMVAssert(m2.rowsize() == size()); 614 assignToU(m2.upperTri(dt())); 615 if (isunit()) m2.diag().setAllTo(T(1)); 616 if (size() > 0) m2.lowerTri().offDiag().setZero(); 617 } 618 619 inline void assignToU(UpperTriMatrixView<RT> m2) const 620 { 621 TMVAssert(m2.size() == size()); 622 TMVAssert(isunit() || !m2.isunit()); 623 TMVAssert(isReal(T())); 624 if (!isSameAs(m2)) Copy(*this,m2); 625 } 626 627 inline void assignToU(UpperTriMatrixView<CT> m2) const 628 { 629 TMVAssert(m2.size() == size()); 630 TMVAssert(isunit() || !m2.isunit()); 631 if (!isSameAs(m2)) Copy(*this,m2); 632 } 633 634 // 635 // subMatrix 636 // 637 638 bool hasSubMatrix( 639 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const; 640 641 bool hasSubVector( 642 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const; 643 644 bool hasSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const; 645 646 inline const_rec_type cSubMatrix( 647 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 648 { 649 return const_rec_type( 650 cptr()+i1*stepi()+j1*stepj(), 651 i2-i1, j2-j1, stepi(), stepj(), ct()); 652 } 653 654 inline const_rec_type subMatrix( 655 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 656 { 657 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 658 return cSubMatrix(i1,i2,j1,j2); 659 } 660 661 inline const_rec_type cSubMatrix( 662 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 663 { 664 return const_rec_type( 665 cptr()+i1*stepi()+j1*stepj(), 666 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(),ct()); 667 } 668 669 inline const_rec_type subMatrix( 670 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 671 { 672 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 673 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 674 } 675 676 inline const_vec_type cSubVector( 677 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 678 { 679 TMVAssert(size >= 0); 680 return const_vec_type( 681 cptr()+i*stepi()+j*stepj(),size, 682 istep*stepi()+jstep*stepj(),ct()); 683 } 684 685 inline const_vec_type subVector( 686 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 687 { 688 TMVAssert(size >= 0); 689 TMVAssert(hasSubVector(i,j,istep,jstep,size)); 690 return cSubVector(i,j,istep,jstep,size); 691 } 692 693 inline const_uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 694 { 695 return const_uppertri_type( 696 cptr()+i1*(stepi()+stepj()), 697 i2-i1,stepi(),stepj(),dt(),ct()); 698 } 699 700 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 701 { 702 TMVAssert(hasSubTriMatrix(i1,i2,1)); 703 return cSubTriMatrix(i1,i2); 704 } 705 706 inline const_uppertri_type cSubTriMatrix( 707 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 708 { 709 return const_uppertri_type( 710 cptr()+i1*(stepi()+stepj()), 711 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(),ct()); 712 } 713 714 inline const_uppertri_type subTriMatrix( 715 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 716 { 717 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 718 return cSubTriMatrix(i1,i2,istep); 719 } 720 721 inline const_uppertri_type offDiag(ptrdiff_t noff=1) const 722 { 723 TMVAssert(noff >= 1); 724 TMVAssert(noff <= size()); 725 return const_uppertri_type( 726 cptr()+noff*stepj(),size()-noff, 727 stepi(),stepj(),NonUnitDiag,ct()); 728 } 729 730 inline const_realpart_type realPart() const 731 { 732 return const_realpart_type( 733 reinterpret_cast<const RT*>(cptr()), size(), 734 isReal(T()) ? stepi() : 2*stepi(), 735 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj); 736 } 737 738 inline const_realpart_type imagPart() const 739 { 740 TMVAssert(isComplex(T())); 741 TMVAssert(!isunit()); 742 // Since Imag of a UnitDiag TriMatrix has 0's on diagonal. 743 return const_realpart_type( 744 reinterpret_cast<const RT*>(cptr())+1, size(), 745 2*stepi(), 2*stepj(), dt(), NonConj); 746 } 747 748 inline const_uppertri_type view() const 749 { 750 return const_uppertri_type( 751 cptr(),size(),stepi(),stepj(),dt(),ct()); 752 } 753 754 inline const_uppertri_type viewAsUnitDiag() const 755 { 756 return const_uppertri_type( 757 cptr(),size(),stepi(),stepj(),UnitDiag,ct()); 758 } 759 760 inline const_lowertri_type transpose() const 761 { 762 return const_lowertri_type( 763 cptr(),size(),stepj(),stepi(),dt(),ct()); 764 } 765 766 inline const_uppertri_type conjugate() const 767 { 768 return const_uppertri_type( 769 cptr(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,ct())); 770 } 771 772 inline const_lowertri_type adjoint() const 773 { 774 return const_lowertri_type( 775 cptr(),size(),stepj(),stepi(),dt(),TMV_ConjOf(T,ct())); 776 } 777 778 inline nonconst_type nonConst() const 779 { 780 const ptrdiff_t n=size(); 781 return nonconst_type( 782 const_cast<T*>(cptr()),n,stepi(),dt(),ct() 783 TMV_FIRSTLAST1(cptr(),row(n-1,n-1,n).end().getP())); 784 } 785 786 787 // 788 // Functions of Matrix 789 // 790 791 T det() const; 792 793 RT logDet(T* sign=0) const; 794 795 inline T trace() const 796 { return isunit() ? T(RT(size())) : diag().sumElements(); } 797 798 T sumElements() const; 799 800 RT sumAbsElements() const; 801 802 RT sumAbs2Elements() const; 803 804 RT norm() const 805 { return normF(); } 806 807 RT normF() const; 808 809 RT normSq(const RT scale = RT(1)) const; 810 811 RT norm1() const; 812 813 RT doNorm2() const; 814 inline RT norm2() const 815 { return doNorm2(); } 816 817 RT doCondition() const; 818 inline RT condition() const 819 { return doCondition(); } 820 821 RT normInf() const; 822 823 RT maxAbsElement() const; 824 RT maxAbs2Element() const; 825 826 bool isSingular() const { return det() == T(0); } 827 828 template <typename T1> 829 void doMakeInverse(UpperTriMatrixView<T1> minv) const; 830 template <typename T1> 831 void doMakeInverse(MatrixView<T1> minv) const; 832 void doMakeInverseATA(MatrixView<T> ata) const; 833 834 inline void makeInverse(MatrixView<T> minv) const 835 { 836 TMVAssert(minv.colsize() == size()); 837 TMVAssert(minv.rowsize() == size()); 838 doMakeInverse(minv); 839 } 840 841 template <typename T1> 842 inline void makeInverse(MatrixView<T1> minv) const 843 { 844 TMVAssert(minv.colsize() == size()); 845 TMVAssert(minv.rowsize() == size()); 846 doMakeInverse(minv); 847 } 848 849 template <typename T1> 850 inline void makeInverse(UpperTriMatrixView<T1> minv) const 851 { 852 TMVAssert(minv.size() == size()); 853 doMakeInverse(minv); 854 } 855 856 QuotXU<T,T> QInverse() const; 857 inline QuotXU<T,T> inverse() const 858 { return QInverse(); } 859 860 inline void makeInverseATA(MatrixView<T> ata) const 861 { 862 TMVAssert(ata.colsize() == size()); 863 TMVAssert(ata.rowsize() == size()); 864 doMakeInverseATA(ata); 865 } 866 867 template <typename T1, int A> 868 inline void makeInverse(UpperTriMatrix<T1,A>& minv) const 869 { 870 TMVAssert(dt()==NonUnitDiag || isunit()); 871 makeInverse(minv.view()); 872 } 873 874 template <typename T1, int A> 875 inline void makeInverse(Matrix<T1,A>& minv) const 876 { makeInverse(minv.view()); } 877 878 template <int A> 879 inline void makeInverseATA(Matrix<T,A>& minv) const 880 { makeInverseATA(minv.view()); } 881 882 // 883 // I/O 884 // 885 886 void write(const TMV_Writer& writer) const; 887 888 // 889 // Arithmetic Helpers 890 // 891 892 template <typename T1> 893 void doLDivEq(VectorView<T1> v) const; 894 template <typename T1, typename T0> 895 void doLDiv(const GenVector<T1>& v1, VectorView<T0> v0) const; 896 template <typename T1> 897 void doLDivEq(MatrixView<T1> m) const; 898 template <typename T1, typename T0> 899 void doLDiv(const GenMatrix<T1>& m1, MatrixView<T0> m0) const; 900 template <typename T1> 901 void doLDivEq(UpperTriMatrixView<T1> m) const; 902 template <typename T1, typename T0> 903 void doLDiv( 904 const GenUpperTriMatrix<T1>& m1, 905 UpperTriMatrixView<T0> m0) const; 906 907 template <typename T1> 908 inline void LDivEq(VectorView<T1> v) const 909 { 910 TMVAssert(v.size() == size()); 911 doLDivEq(v); 912 } 913 template <typename T1, typename T0> 914 inline void LDiv( 915 const GenVector<T1>& v1, VectorView<T0> v0) const 916 { 917 TMVAssert(v0.size() == size()); 918 TMVAssert(v1.size() == size()); 919 doLDiv(v1,v0); 920 } 921 template <typename T1> 922 inline void RDivEq(VectorView<T1> v) const 923 { transpose().LDivEq(v); } 924 template <typename T1, typename T0> 925 inline void RDiv( 926 const GenVector<T1>& v1, VectorView<T0> v0) const 927 { transpose().LDiv(v1,v0); } 928 929 template <typename T1> 930 inline void LDivEq(MatrixView<T1> m) const 931 { 932 TMVAssert(m.colsize() == size()); 933 doLDivEq(m); 934 } 935 template <typename T1, typename T0> 936 inline void LDiv( 937 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 938 { 939 TMVAssert(m0.colsize() == size()); 940 TMVAssert(m1.colsize() == size()); 941 TMVAssert(m0.rowsize() == m1.rowsize()); 942 doLDiv(m1,m0); 943 } 944 template <typename T1> 945 inline void RDivEq(MatrixView<T1> m) const 946 { transpose().LDivEq(m.transpose()); } 947 template <typename T1, typename T0> 948 inline void RDiv( 949 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 950 { transpose().LDiv(m1.transpose(),m0.transpose()); } 951 952 template <typename T1> 953 inline void LDivEq(UpperTriMatrixView<T1> m) const 954 { 955 TMVAssert(m.colsize() == size()); 956 doLDivEq(m); 957 } 958 template <typename T1, typename T0> 959 inline void LDiv( 960 const GenUpperTriMatrix<T1>& m1, 961 UpperTriMatrixView<T0> m0) const 962 { 963 TMVAssert(m0.size() == size()); 964 TMVAssert(m1.size() == size()); 965 doLDiv(m1,m0); 966 } 967 template <typename T1> 968 inline void RDivEq(UpperTriMatrixView<T1> m) const 969 { transpose().LDivEq(m.transpose()); } 970 template <typename T1, typename T0> 971 inline void RDiv( 972 const GenUpperTriMatrix<T1>& m1, 973 UpperTriMatrixView<T0> m0) const 974 { transpose().LDiv(m1.transpose(),m0.transpose()); } 975 976 // For easier compatibility with regular matrices: 977 inline void divideInPlace() const {} 978 inline void saveDiv() const {} 979 inline void setDiv() const {} 980 inline void unsetDiv() const {} 981 inline void resetDiv() const {} 982 inline bool divIsSet() const { return true; } 983 inline void divideUsing(DivType TMV_DEBUGPARAM(dt)) const 984 { TMVAssert(dt == LU); } 985 986 virtual const T* cptr() const = 0; 987 virtual ptrdiff_t stepi() const = 0; 988 virtual ptrdiff_t stepj() const = 0; 989 virtual ConjType ct() const = 0; 990 inline bool isrm() const { return stepj() == 1; } 991 inline bool iscm() const { return stepi() == 1; } 992 inline bool isunit() const { return dt() == UnitDiag; } 993 inline bool isconj() const 994 { 995 TMVAssert(isComplex(T()) || ct()==NonConj); 996 return isComplex(T()) && ct()==Conj; 997 } 998 999 virtual inline T cref(ptrdiff_t i, ptrdiff_t j) const 1000 { 1001 if (isunit() && i==j) return T(1); 1002 else if (i>j) return T(0); 1003 else { 1004 const T* mi = cptr() + i*stepi() + j*stepj(); 1005 return (isconj() ? TMV_CONJ(*mi) : *mi); 1006 } 1007 } 1008 1009 inline ptrdiff_t rowstart(ptrdiff_t i) const { return i; } 1010 inline ptrdiff_t rowend(ptrdiff_t ) const { return size(); } 1011 1012 inline ptrdiff_t colstart(ptrdiff_t ) const { return 0; } 1013 inline ptrdiff_t colend(ptrdiff_t j) const { return j+1; } 1014 1015 inline const_rowmajor_iterator rowmajor_begin() const 1016 { return const_rowmajor_iterator(this,0,0); } 1017 inline const_rowmajor_iterator rowmajor_end() const 1018 { return const_rowmajor_iterator(this,size(),size()); } 1019 1020 inline const_colmajor_iterator colmajor_begin() const 1021 { return const_colmajor_iterator(this,0,0); } 1022 inline const_colmajor_iterator colmajor_end() const 1023 { return const_colmajor_iterator(this,0,size()); } 1024 1025 protected : 1026 1027 inline bool okij(ptrdiff_t i, ptrdiff_t j) const 1028 { 1029 TMVAssert(i>=0 && i < size()); 1030 TMVAssert(j>=0 && j < size()); 1031 return isunit() ? i-j<0 : i-j<=0; 1032 } 1033 1034 private : 1035 1036 type& operator=(const type&); 1037 1038 }; // GenUpperTriMatrix 1039 1040 template <typename T> 1041 class GenLowerTriMatrix : 1042 virtual public AssignableToLowerTriMatrix<T>, 1043 public BaseMatrix<T> 1044 { 1045 public: 1046 1047 typedef TMV_RealType(T) RT; 1048 typedef TMV_ComplexType(T) CT; 1049 typedef T value_type; 1050 typedef RT real_type; 1051 typedef CT complex_type; 1052 typedef GenLowerTriMatrix<T> type; 1053 typedef LowerTriMatrix<T> copy_type; 1054 typedef ConstVectorView<T> const_vec_type; 1055 typedef ConstMatrixView<T> const_rec_type; 1056 typedef ConstUpperTriMatrixView<T> const_uppertri_type; 1057 typedef ConstLowerTriMatrixView<T> const_lowertri_type; 1058 typedef const_lowertri_type const_view_type; 1059 typedef const_uppertri_type const_transpose_type; 1060 typedef const_lowertri_type const_conjugate_type; 1061 typedef const_uppertri_type const_adjoint_type; 1062 typedef ConstLowerTriMatrixView<RT> const_realpart_type; 1063 typedef LowerTriMatrixView<T> nonconst_type; 1064 typedef CRMIt<type> const_rowmajor_iterator; 1065 typedef CCMIt<type> const_colmajor_iterator; 1066 1067 // 1068 // Constructors 1069 // 1070 1071 inline GenLowerTriMatrix() {} 1072 inline GenLowerTriMatrix(const type&) {} 1073 virtual inline ~GenLowerTriMatrix() {} 1074 1075 // 1076 // Access Functions 1077 // 1078 1079 using AssignableToLowerTriMatrix<T>::size; 1080 using AssignableToLowerTriMatrix<T>::dt; 1081 inline ptrdiff_t colsize() const { return size(); } 1082 inline ptrdiff_t rowsize() const { return size(); } 1083 1084 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 1085 { 1086 TMVAssert(i>=0 && i<size()); 1087 TMVAssert(j>=0 && j<size()); 1088 return cref(i,j); 1089 } 1090 1091 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1092 { 1093 TMVAssert(i>=0 && i<size()); 1094 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 1095 TMVAssert(j1==j2 || okij(i,j2-1)); 1096 return const_vec_type( 1097 cptr()+i*stepi()+j1*stepj(),j2-j1,stepj(),ct()); 1098 } 1099 1100 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 1101 { 1102 TMVAssert(j>=0 && j<size()); 1103 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 1104 TMVAssert(i1==i2 || okij(i1,j)); 1105 return const_vec_type( 1106 cptr()+i1*stepi()+j*stepj(),i2-i1,stepi(),ct()); 1107 } 1108 1109 inline const_vec_type diag() const 1110 { 1111 TMVAssert(!isunit()); 1112 return const_vec_type(cptr(),size(),stepi()+stepj(),ct()); 1113 } 1114 1115 inline const_vec_type diag(ptrdiff_t i) const 1116 { 1117 TMVAssert(i>=-size()); 1118 TMVAssert(isunit() ? i<0 : i<=0); 1119 return const_vec_type( 1120 cptr()-i*stepi(),size()-i,stepi()+stepj(),ct()); 1121 } 1122 1123 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1124 { 1125 TMVAssert(i>=-size()); 1126 TMVAssert(isunit() ? i<0 : i<=0); 1127 TMVAssert(j1>=0 && j1 <= j2 && j2 <= size()+i); 1128 const ptrdiff_t ds = stepi()+stepj(); 1129 return const_vec_type(cptr()-i*stepi()+j1*ds,j2-j1,ds,ct()); 1130 } 1131 1132 template <typename T2> 1133 inline bool isSameAs(const BaseMatrix<T2>& ) const 1134 { return false; } 1135 1136 inline bool isSameAs(const GenLowerTriMatrix<T>& m2) const 1137 { 1138 if (this == &m2) return true; 1139 else return (cptr()==m2.cptr() && size()==m2.size() && 1140 dt() == m2.dt() && ct() == m2.ct() && 1141 stepi()==m2.stepi() && stepj()==m2.stepj()); 1142 } 1143 1144 inline void assignToM(MatrixView<RT> m2) const 1145 { transpose().assignToM(m2.transpose()); } 1146 inline void assignToM(MatrixView<CT> m2) const 1147 { transpose().assignToM(m2.transpose()); } 1148 inline void assignToL(LowerTriMatrixView<RT> m2) const 1149 { transpose().assignToU(m2.transpose()); } 1150 inline void assignToL(LowerTriMatrixView<CT> m2) const 1151 { transpose().assignToU(m2.transpose()); } 1152 1153 // 1154 // subMatrix 1155 // 1156 1157 inline bool hasSubMatrix( 1158 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1159 { return transpose().hasSubMatrix(j1,j2,i1,i2,jstep,istep); } 1160 1161 inline bool hasSubVector( 1162 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 1163 { return transpose().hasSubVector(j,i,jstep,istep,size); } 1164 1165 inline bool hasSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1166 { return transpose().hasSubTriMatrix(i1,i2,istep); } 1167 1168 inline const_rec_type cSubMatrix( 1169 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 1170 { 1171 return const_rec_type( 1172 cptr()+i1*stepi()+j1*stepj(), 1173 i2-i1, j2-j1, stepi(), stepj(), ct()); 1174 } 1175 1176 inline const_rec_type subMatrix( 1177 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 1178 { 1179 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 1180 return cSubMatrix(i1,i2,j1,j2); 1181 } 1182 1183 inline const_rec_type cSubMatrix( 1184 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1185 { 1186 return const_rec_type( 1187 cptr()+i1*stepi()+j1*stepj(), 1188 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(),ct()); 1189 } 1190 1191 inline const_rec_type subMatrix( 1192 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1193 { 1194 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 1195 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 1196 } 1197 1198 inline const_vec_type cSubVector( 1199 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 1200 { 1201 return const_vec_type( 1202 cptr()+i*stepi()+j*stepj(),size, 1203 istep*stepi()+jstep*stepj(),ct()); 1204 } 1205 1206 inline const_vec_type subVector( 1207 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 1208 { 1209 TMVAssert(size >= 0); 1210 TMVAssert(hasSubVector(i,j,istep,jstep,size)); 1211 return cSubVector(i,j,istep,jstep,size); 1212 } 1213 1214 inline const_lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1215 { 1216 return const_lowertri_type( 1217 cptr()+i1*(stepi()+stepj()), 1218 i2-i1,stepi(),stepj(),dt(),ct()); 1219 } 1220 1221 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1222 { 1223 TMVAssert(hasSubTriMatrix(i1,i2,1)); 1224 return cSubTriMatrix(i1,i2); 1225 } 1226 1227 inline const_lowertri_type cSubTriMatrix( 1228 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1229 { 1230 return const_lowertri_type( 1231 cptr()+i1*(stepi()+stepj()), 1232 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(),ct()); 1233 } 1234 1235 inline const_lowertri_type subTriMatrix( 1236 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1237 { 1238 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 1239 return cSubTriMatrix(i1,i2,istep); 1240 } 1241 1242 inline const_lowertri_type offDiag(ptrdiff_t noff=1) const 1243 { 1244 TMVAssert(noff >= 1); 1245 TMVAssert(noff <= size()); 1246 return const_lowertri_type( 1247 cptr()+noff*stepi(),size()-noff, 1248 stepi(),stepj(),NonUnitDiag,ct()); 1249 } 1250 1251 inline const_realpart_type realPart() const 1252 { 1253 return const_realpart_type( 1254 reinterpret_cast<const RT*>(cptr()), size(), 1255 isReal(T()) ? stepi() : 2*stepi(), 1256 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj); 1257 } 1258 1259 inline const_realpart_type imagPart() const 1260 { 1261 TMVAssert(isComplex(T())); 1262 TMVAssert(!isunit()); 1263 return const_realpart_type( 1264 reinterpret_cast<const RT*>(cptr())+1, size(), 1265 2*stepi(), 2*stepj(), dt(), NonConj); 1266 } 1267 1268 inline const_lowertri_type view() const 1269 { 1270 return const_lowertri_type( 1271 cptr(),size(),stepi(),stepj(),dt(),ct()); 1272 } 1273 1274 inline const_lowertri_type viewAsUnitDiag() const 1275 { 1276 return const_lowertri_type( 1277 cptr(),size(),stepi(),stepj(),UnitDiag,ct()); 1278 } 1279 1280 inline const_uppertri_type transpose() const 1281 { 1282 return const_uppertri_type( 1283 cptr(),size(),stepj(),stepi(),dt(),ct()); 1284 } 1285 1286 inline const_lowertri_type conjugate() const 1287 { 1288 return const_lowertri_type( 1289 cptr(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,ct())); 1290 } 1291 1292 inline const_uppertri_type adjoint() const 1293 { 1294 return const_uppertri_type( 1295 cptr(),size(),stepj(),stepi(),dt(),TMV_ConjOf(T,ct())); 1296 } 1297 1298 inline nonconst_type nonConst() const 1299 { 1300 const ptrdiff_t n=size(); 1301 return nonconst_type( 1302 const_cast<T*>(cptr()),n,stepi(),dt(),ct() 1303 TMV_FIRSTLAST1(cptr(),row(n-1,n-1,n).end().getP())); 1304 } 1305 1306 // 1307 // Functions of Matrix 1308 // 1309 1310 T det() const 1311 { return transpose().det(); } 1312 1313 RT logDet(T* sign=0) const 1314 { return transpose().logDet(sign); } 1315 1316 inline T trace() const 1317 { return isunit() ? T(RT(size())) : diag().sumElements(); } 1318 1319 inline T sumElements() const 1320 { return transpose().sumElements(); } 1321 1322 inline RT sumAbsElements() const 1323 { return transpose().sumAbsElements(); } 1324 1325 inline RT sumAbs2Elements() const 1326 { return transpose().sumAbs2Elements(); } 1327 1328 inline RT norm() const 1329 { return transpose().normF(); } 1330 1331 inline RT normF() const 1332 { return transpose().normF(); } 1333 1334 inline RT normSq(const RT scale = RT(1)) const 1335 { return transpose().normSq(scale); } 1336 1337 inline RT norm1() const 1338 { return transpose().normInf(); } 1339 1340 inline RT doNorm2() const 1341 { return transpose().doNorm2(); } 1342 1343 inline RT doCondition() const 1344 { return transpose().doCondition(); } 1345 1346 inline RT norm2() const 1347 { return transpose().norm2(); } 1348 1349 inline RT condition() const 1350 { return transpose().condition(); } 1351 1352 inline RT normInf() const 1353 { return transpose().norm1(); } 1354 1355 inline RT maxAbsElement() const 1356 { return transpose().maxAbsElement(); } 1357 1358 inline RT maxAbs2Element() const 1359 { return transpose().maxAbs2Element(); } 1360 1361 bool isSingular() const { return det() == T(0); } 1362 1363 QuotXL<T,T> QInverse() const; 1364 inline QuotXL<T,T> inverse() const 1365 { return QInverse(); } 1366 1367 template <typename T1> 1368 inline void makeInverse(LowerTriMatrixView<T1> minv) const 1369 { 1370 TMVAssert(minv.size() == size()); 1371 transpose().makeInverse(minv.transpose()); 1372 } 1373 1374 inline void makeInverse(MatrixView<T> minv) const 1375 { 1376 TMVAssert(minv.colsize() == size()); 1377 TMVAssert(minv.rowsize() == size()); 1378 transpose().makeInverse(minv.transpose()); 1379 } 1380 1381 template <typename T1> 1382 inline void makeInverse(MatrixView<T1> minv) const 1383 { 1384 TMVAssert(minv.colsize() == size()); 1385 TMVAssert(minv.rowsize() == size()); 1386 transpose().makeInverse(minv.transpose()); 1387 } 1388 1389 void doMakeInverseATA(MatrixView<T> minv) const; 1390 1391 inline void makeInverseATA(MatrixView<T> ata) const 1392 { 1393 TMVAssert(ata.colsize() == size()); 1394 TMVAssert(ata.rowsize() == size()); 1395 doMakeInverseATA(ata); 1396 } 1397 1398 template <typename T1, int A> 1399 inline void makeInverse(LowerTriMatrix<T1,A>& minv) const 1400 { 1401 TMVAssert(dt()==NonUnitDiag || isunit()); 1402 makeInverse(minv.view()); 1403 } 1404 1405 template <typename T1, int A> 1406 inline void makeInverse(Matrix<T1,A>& minv) const 1407 { makeInverse(minv.view()); } 1408 1409 template <int A> 1410 inline void makeInverseATA(Matrix<T,A>& minv) const 1411 { makeInverseATA(minv.view()); } 1412 1413 // 1414 // I/O 1415 // 1416 1417 void write(const TMV_Writer& writer) const; 1418 1419 // 1420 // Arithmetic Helpers 1421 // 1422 1423 template <typename T1> 1424 void doLDivEq(VectorView<T1> v) const; 1425 template <typename T1, typename T0> 1426 void doLDiv(const GenVector<T1>& v1, VectorView<T0> v0) const; 1427 template <typename T1> 1428 void doLDivEq(MatrixView<T1> m) const; 1429 template <typename T1, typename T0> 1430 void doLDiv(const GenMatrix<T1>& m1, MatrixView<T0> m0) const; 1431 template <typename T1> 1432 void doLDivEq(LowerTriMatrixView<T1> m) const; 1433 template <typename T1, typename T0> 1434 void doLDiv( 1435 const GenLowerTriMatrix<T1>& m1, 1436 LowerTriMatrixView<T0> m0) const; 1437 1438 template <typename T1> 1439 inline void LDivEq(VectorView<T1> v) const 1440 { 1441 TMVAssert(v.size() == size()); 1442 doLDivEq(v); 1443 } 1444 template <typename T1, typename T0> 1445 inline void LDiv(const GenVector<T1>& v1, VectorView<T0> v0) const 1446 { 1447 TMVAssert(v0.size() == size()); 1448 TMVAssert(v1.size() == size()); 1449 doLDiv(v1,v0); 1450 } 1451 template <typename T1> 1452 inline void RDivEq(VectorView<T1> v) const 1453 { transpose().LDivEq(v); } 1454 template <typename T1, typename T0> 1455 inline void RDiv(const GenVector<T1>& v1, VectorView<T0> v0) const 1456 { transpose().LDiv(v1,v0); } 1457 1458 template <typename T1> 1459 inline void LDivEq(MatrixView<T1> m) const 1460 { 1461 TMVAssert(m.colsize() == size()); 1462 doLDivEq(m); 1463 } 1464 template <typename T1, typename T0> 1465 inline void LDiv( 1466 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 1467 { 1468 TMVAssert(m0.colsize() == size()); 1469 TMVAssert(m1.colsize() == size()); 1470 TMVAssert(m0.rowsize() == m1.rowsize()); 1471 doLDiv(m1,m0); 1472 } 1473 template <typename T1> 1474 inline void RDivEq(MatrixView<T1> m) const 1475 { transpose().LDivEq(m.transpose()); } 1476 template <typename T1, typename T0> 1477 inline void RDiv( 1478 const GenMatrix<T1>& m1, MatrixView<T0> m0) const 1479 { transpose().LDiv(m1.transpose(),m0.transpose()); } 1480 1481 template <typename T1> 1482 inline void LDivEq(LowerTriMatrixView<T1> m) const 1483 { 1484 TMVAssert(m.colsize() == size()); 1485 doLDivEq(m); 1486 } 1487 template <typename T1, typename T0> 1488 inline void LDiv( 1489 const GenLowerTriMatrix<T1>& m1, 1490 LowerTriMatrixView<T0> m0) const 1491 { 1492 TMVAssert(m0.size() == size()); 1493 TMVAssert(m1.size() == size()); 1494 doLDiv(m1,m0); 1495 } 1496 template <typename T1> 1497 inline void RDivEq(LowerTriMatrixView<T1> m) const 1498 { transpose().LDivEq(m.transpose()); } 1499 template <typename T1, typename T0> 1500 inline void RDiv( 1501 const GenLowerTriMatrix<T1>& m1, 1502 LowerTriMatrixView<T0> m0) const 1503 { transpose().LDiv(m1.transpose(),m0.transpose()); } 1504 1505 1506 // For easier compatibility with regular matrices: 1507 inline void divideInPlace() const {} 1508 inline void saveDiv() const {} 1509 inline void setDiv() const {} 1510 inline void unsetDiv() const {} 1511 inline void resetDiv() const {} 1512 inline bool divIsSet() const { return true; } 1513 inline void divideUsing(DivType TMV_DEBUGPARAM(dt)) const 1514 { TMVAssert(dt == LU); } 1515 inline bool checkDecomp(std::ostream* fout=0) const { return true; } 1516 inline bool checkDecomp( 1517 const BaseMatrix<T>& m2, std::ostream* fout=0) const 1518 { return true; } 1519 1520 virtual const T* cptr() const = 0; 1521 virtual ptrdiff_t stepi() const = 0; 1522 virtual ptrdiff_t stepj() const = 0; 1523 virtual ConjType ct() const = 0; 1524 inline bool isrm() const { return stepj() == 1; } 1525 inline bool iscm() const { return stepi() == 1; } 1526 inline bool isunit() const { return dt() == UnitDiag; } 1527 inline bool isconj() const 1528 { 1529 TMVAssert(isComplex(T()) || ct()==NonConj); 1530 return isComplex(T()) && ct()==Conj; 1531 } 1532 1533 virtual inline T cref(ptrdiff_t i, ptrdiff_t j) const 1534 { 1535 if (isunit() && i==j) return T(1); 1536 else if (i<j) return T(0); 1537 else { 1538 const T* mi = cptr() + i*stepi() + j*stepj(); 1539 return (isconj() ? TMV_CONJ(*mi) : *mi); 1540 } 1541 } 1542 1543 inline ptrdiff_t rowstart(ptrdiff_t ) const { return 0; } 1544 inline ptrdiff_t rowend(ptrdiff_t i) const { return i+1; } 1545 1546 inline ptrdiff_t colstart(ptrdiff_t j) const { return j; } 1547 inline ptrdiff_t colend(ptrdiff_t ) const { return size(); } 1548 1549 inline const_rowmajor_iterator rowmajor_begin() const 1550 { return const_rowmajor_iterator(this,0,0); } 1551 inline const_rowmajor_iterator rowmajor_end() const 1552 { return const_rowmajor_iterator(this,size(),0); } 1553 1554 inline const_colmajor_iterator colmajor_begin() const 1555 { return const_colmajor_iterator(this,0,0); } 1556 inline const_colmajor_iterator colmajor_end() const 1557 { return const_colmajor_iterator(this,size(),size()); } 1558 1559 1560 protected : 1561 1562 inline bool okij(ptrdiff_t i, ptrdiff_t j) const 1563 { 1564 TMVAssert(i>=0 && i < size()); 1565 TMVAssert(j>=0 && j < size()); 1566 return isunit() ? i-j>0 : i-j>=0; 1567 } 1568 1569 private : 1570 1571 type& operator=(const type&); 1572 1573 }; // GenLowerTriMatrix 1574 1575 template <typename T, int A> 1576 class ConstUpperTriMatrixView : public GenUpperTriMatrix<T> 1577 { 1578 public : 1579 1580 typedef GenUpperTriMatrix<T> base; 1581 typedef ConstUpperTriMatrixView<T,A> type; 1582 1583 inline ConstUpperTriMatrixView(const type& rhs) : 1584 itsm(rhs.itsm), itss(rhs.itss), itssi(rhs.itssi), itssj(rhs.itssj), 1585 itsdiag(rhs.itsdiag), itsct(rhs.itsct) {} 1586 1587 inline ConstUpperTriMatrixView(const base& rhs) : 1588 itsm(rhs.cptr()), itss(rhs.size()), 1589 itssi(rhs.stepi()), itssj(rhs.stepj()), 1590 itsdiag(rhs.dt()), itsct(rhs.ct()) {} 1591 1592 inline ConstUpperTriMatrixView( 1593 const T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 1594 DiagType _dt, ConjType _ct) : 1595 itsm(_m), itss(_s), itssi(_si), itssj(_sj), 1596 itsdiag(_dt), itsct(_ct) {} 1597 1598 virtual inline ~ConstUpperTriMatrixView() 1599 { 1600 #ifdef TMV_EXTRA_DEBUG 1601 const_cast<const T*&>(itsm) = 0; 1602 #endif 1603 } 1604 1605 inline ptrdiff_t size() const { return itss; } 1606 inline const T* cptr() const { return itsm; } 1607 inline ptrdiff_t stepi() const { return itssi; } 1608 inline ptrdiff_t stepj() const { return itssj; } 1609 inline DiagType dt() const { return itsdiag; } 1610 inline ConjType ct() const { return itsct; } 1611 1612 protected : 1613 1614 const T*const itsm; 1615 const ptrdiff_t itss; 1616 const ptrdiff_t itssi; 1617 const ptrdiff_t itssj; 1618 DiagType itsdiag; 1619 ConjType itsct; 1620 1621 private : 1622 1623 type& operator=(const type&); 1624 1625 }; // ConstUpperTriMatrixView 1626 1627 template <typename T, int A> 1628 class ConstLowerTriMatrixView : public GenLowerTriMatrix<T> 1629 { 1630 public : 1631 1632 typedef GenLowerTriMatrix<T> base; 1633 typedef ConstLowerTriMatrixView<T,A> type; 1634 1635 inline ConstLowerTriMatrixView(const type& rhs) : 1636 itsm(rhs.itsm), itss(rhs.itss), itssi(rhs.itssi), itssj(rhs.itssj), 1637 itsdiag(rhs.itsdiag), itsct(rhs.itsct) {} 1638 1639 inline ConstLowerTriMatrixView(const GenLowerTriMatrix<T>& rhs) : 1640 itsm(rhs.cptr()), itss(rhs.size()), 1641 itssi(rhs.stepi()), itssj(rhs.stepj()), 1642 itsdiag(rhs.dt()), itsct(rhs.ct()) {} 1643 1644 inline ConstLowerTriMatrixView( 1645 const T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 1646 DiagType _dt, ConjType _ct) : 1647 itsm(_m), itss(_s), itssi(_si), itssj(_sj), 1648 itsdiag(_dt), itsct(_ct) {} 1649 1650 virtual inline ~ConstLowerTriMatrixView() 1651 { 1652 #ifdef TMV_EXTRA_DEBUG 1653 const_cast<const T*&>(itsm) = 0; 1654 #endif 1655 } 1656 1657 inline ptrdiff_t size() const { return itss; } 1658 inline const T* cptr() const { return itsm; } 1659 inline ptrdiff_t stepi() const { return itssi; } 1660 inline ptrdiff_t stepj() const { return itssj; } 1661 inline DiagType dt() const { return itsdiag; } 1662 inline ConjType ct() const { return itsct; } 1663 1664 protected : 1665 1666 const T*const itsm; 1667 const ptrdiff_t itss; 1668 const ptrdiff_t itssi; 1669 const ptrdiff_t itssj; 1670 DiagType itsdiag; 1671 ConjType itsct; 1672 1673 private : 1674 1675 type& operator=(const type&); 1676 1677 }; // ConstLowerTriMatrixView 1678 1679 template <typename T> 1680 class ConstUpperTriMatrixView<T,FortranStyle> : 1681 public ConstUpperTriMatrixView<T,CStyle> 1682 { 1683 public : 1684 1685 typedef TMV_RealType(T) RT; 1686 typedef GenUpperTriMatrix<T> base; 1687 typedef ConstUpperTriMatrixView<T,FortranStyle> type; 1688 typedef ConstUpperTriMatrixView<T,CStyle> c_type; 1689 typedef ConstVectorView<T,FortranStyle> const_vec_type; 1690 typedef ConstMatrixView<T,FortranStyle> const_rec_type; 1691 typedef ConstLowerTriMatrixView<T,FortranStyle> const_lowertri_type; 1692 typedef ConstUpperTriMatrixView<T,FortranStyle> const_uppertri_type; 1693 typedef const_uppertri_type const_view_type; 1694 typedef const_lowertri_type const_transpose_type; 1695 typedef const_uppertri_type const_conjugate_type; 1696 typedef const_lowertri_type const_adjoint_type; 1697 typedef ConstUpperTriMatrixView<RT,FortranStyle> const_realpart_type; 1698 1699 inline ConstUpperTriMatrixView(const type& rhs) : c_type(rhs) {} 1700 1701 inline ConstUpperTriMatrixView(const c_type& rhs) : c_type(rhs) {} 1702 1703 inline ConstUpperTriMatrixView(const base& rhs) : c_type(rhs) {} 1704 1705 inline ConstUpperTriMatrixView( 1706 const T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 1707 DiagType _dt, ConjType _ct) : 1708 c_type(_m,_s,_si,_sj,_dt,_ct) 1709 {} 1710 1711 virtual inline ~ConstUpperTriMatrixView() {} 1712 1713 // 1714 // Access Functions 1715 // 1716 1717 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 1718 { 1719 TMVAssert(i>0 && i<=c_type::size()); 1720 TMVAssert(j>0 && j<=c_type::size()); 1721 return c_type::cref(i-1,j-1); 1722 } 1723 1724 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1725 { 1726 TMVAssert(i>0 && i<=c_type::size()); 1727 TMVAssert(j1>0 && j1<=j2 && j2<=c_type::size()); 1728 return base::row(i-1,j1-1,j2); 1729 } 1730 1731 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 1732 { 1733 TMVAssert(j>0 && j<=c_type::size()); 1734 TMVAssert(i1>0 && i1<=i2 && i2<=c_type::size()); 1735 return base::col(j-1,i1-1,i2); 1736 } 1737 1738 inline const_vec_type diag() const 1739 { return base::diag(); } 1740 1741 inline const_vec_type diag(ptrdiff_t i) const 1742 { return base::diag(i); } 1743 1744 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1745 { 1746 TMVAssert(j1>0); 1747 return base::diag(i,j1-1,j2); 1748 } 1749 1750 // 1751 // subMatrix 1752 // 1753 1754 bool hasSubMatrix( 1755 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const; 1756 1757 bool hasSubVector( 1758 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const; 1759 1760 bool hasSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const; 1761 1762 inline const_rec_type subMatrix( 1763 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 1764 { 1765 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 1766 return base::cSubMatrix(i1-1,i2,j1-1,j2); 1767 } 1768 1769 inline const_rec_type subMatrix( 1770 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1771 { 1772 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 1773 return base::cSubMatrix( 1774 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 1775 } 1776 1777 inline const_vec_type subVector( 1778 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 1779 { 1780 TMVAssert(size >= 0); 1781 TMVAssert(hasSubVector(i,j,istep,jstep,size)); 1782 return base::cSubVector(i-1,j-1,istep,jstep,size); 1783 } 1784 1785 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1786 { 1787 TMVAssert(hasSubTriMatrix(i1,i2,1)); 1788 return base::cSubTriMatrix(i1-1,i2); 1789 } 1790 1791 inline const_uppertri_type subTriMatrix( 1792 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1793 { 1794 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 1795 return base::cSubTriMatrix(i1-1,i2-1+istep,istep); 1796 } 1797 1798 inline const_uppertri_type offDiag(ptrdiff_t noff=1) const 1799 { return base::offDiag(noff); } 1800 1801 inline const_realpart_type realPart() const 1802 { return base::realPart(); } 1803 1804 inline const_realpart_type imagPart() const 1805 { return base::imagPart(); } 1806 1807 inline const_uppertri_type view() const 1808 { return base::view(); } 1809 1810 inline const_uppertri_type viewAsUnitDiag() const 1811 { return base::viewAsUnitDiag(); } 1812 1813 inline const_lowertri_type transpose() const 1814 { return base::transpose(); } 1815 1816 inline const_uppertri_type conjugate() const 1817 { return base::conjugate(); } 1818 1819 inline const_lowertri_type adjoint() const 1820 { return base::adjoint(); } 1821 1822 private : 1823 1824 type& operator=(const type&); 1825 1826 }; // FortranStyle ConstUpperTriMatrixView 1827 1828 template <typename T> 1829 class ConstLowerTriMatrixView<T,FortranStyle> : 1830 public ConstLowerTriMatrixView<T,CStyle> 1831 { 1832 public : 1833 1834 typedef TMV_RealType(T) RT; 1835 typedef GenLowerTriMatrix<T> base; 1836 typedef ConstLowerTriMatrixView<T,FortranStyle> type; 1837 typedef ConstLowerTriMatrixView<T,CStyle> c_type; 1838 typedef ConstVectorView<T,FortranStyle> const_vec_type; 1839 typedef ConstMatrixView<T,FortranStyle> const_rec_type; 1840 typedef ConstUpperTriMatrixView<T,FortranStyle> const_uppertri_type; 1841 typedef ConstLowerTriMatrixView<T,FortranStyle> const_lowertri_type; 1842 typedef const_lowertri_type const_view_type; 1843 typedef const_uppertri_type const_transpose_type; 1844 typedef const_lowertri_type const_conjugate_type; 1845 typedef const_uppertri_type const_adjoint_type; 1846 typedef ConstLowerTriMatrixView<RT,FortranStyle> const_realpart_type; 1847 1848 inline ConstLowerTriMatrixView(const type& rhs) : c_type(rhs) {} 1849 1850 inline ConstLowerTriMatrixView(const c_type& rhs) : c_type(rhs) {} 1851 1852 inline ConstLowerTriMatrixView(const base& rhs) : c_type(rhs) {} 1853 1854 inline ConstLowerTriMatrixView( 1855 const T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 1856 DiagType _dt, ConjType _ct) : 1857 c_type(_m,_s,_si,_sj,_dt,_ct) 1858 {} 1859 1860 virtual inline ~ConstLowerTriMatrixView() {} 1861 1862 // 1863 // Access Functions 1864 // 1865 1866 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 1867 { 1868 TMVAssert(i>0 && i<=c_type::size()); 1869 TMVAssert(j>0 && j<=c_type::size()); 1870 return c_type::cref(i-1,j-1); 1871 } 1872 1873 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1874 { 1875 TMVAssert(i>0 && i<=c_type::size()); 1876 TMVAssert(j1>0 && j1<=j2 && j2<=c_type::size()); 1877 return base::row(i-1,j1-1,j2); 1878 } 1879 1880 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 1881 { 1882 TMVAssert(j>0 && j<=c_type::size()); 1883 TMVAssert(i1>0 && i1<=i2 && i2<=c_type::size()); 1884 return base::col(j-1,i1-1,i2); 1885 } 1886 1887 inline const_vec_type diag() const 1888 { return base::diag(); } 1889 1890 inline const_vec_type diag(ptrdiff_t i) const 1891 { return base::diag(i); } 1892 1893 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1894 { 1895 TMVAssert(j1>0); 1896 return base::diag(i,j1-1,j2); 1897 } 1898 1899 // 1900 // subMatrix 1901 // 1902 1903 inline bool hasSubMatrix( 1904 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1905 { return transpose().hasSubMatrix(j1,j2,i1,i2,jstep,istep); } 1906 1907 inline bool hasSubVector( 1908 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 1909 { return transpose().hasSubVector(j,i,jstep,istep,size); } 1910 1911 inline bool hasSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1912 { return transpose().hasSubTriMatrix(i1,i2,istep); } 1913 1914 inline const_rec_type subMatrix( 1915 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 1916 { 1917 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 1918 return base::cSubMatrix(i1-1,i2,j1-1,j2); 1919 } 1920 1921 inline const_rec_type subMatrix( 1922 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1923 { 1924 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 1925 return base::cSubMatrix( 1926 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 1927 } 1928 1929 inline const_vec_type subVector( 1930 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 1931 { 1932 TMVAssert(size >= 0); 1933 TMVAssert(hasSubVector(i,j,istep,jstep,size)); 1934 return base::cSubVector(i-1,j-1,istep,jstep,size); 1935 } 1936 1937 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1938 { 1939 TMVAssert(hasSubTriMatrix(i1,i2,1)); 1940 return base::cSubTriMatrix(i1-1,i2); 1941 } 1942 1943 inline const_lowertri_type subTriMatrix( 1944 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1945 { 1946 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 1947 return base::cSubTriMatrix(i1-1,i2-1+istep,istep); 1948 } 1949 1950 inline const_lowertri_type offDiag(ptrdiff_t noff=1) const 1951 { return base::offDiag(noff); } 1952 1953 inline const_realpart_type realPart() const 1954 { return base::realPart(); } 1955 1956 inline const_realpart_type imagPart() const 1957 { return base::imagPart(); } 1958 1959 inline const_lowertri_type view() const 1960 { return base::view(); } 1961 1962 inline const_lowertri_type viewAsUnitDiag() const 1963 { return base::viewAsUnitDiag(); } 1964 1965 inline const_uppertri_type transpose() const 1966 { return base::transpose(); } 1967 1968 inline const_lowertri_type conjugate() const 1969 { return base::conjugate(); } 1970 1971 inline const_uppertri_type adjoint() const 1972 { return base::adjoint(); } 1973 1974 private : 1975 1976 type& operator=(const type&); 1977 1978 }; // FortranStyle ConstLowerTriMatrixView 1979 1980 template <typename T, int A> 1981 class UpperTriMatrixView : public GenUpperTriMatrix<T> 1982 { 1983 public: 1984 1985 typedef TMV_RealType(T) RT; 1986 typedef TMV_ComplexType(T) CT; 1987 typedef GenUpperTriMatrix<T> base; 1988 typedef UpperTriMatrixView<T,A> type; 1989 typedef LowerTriMatrixView<T,A> lowertri_type; 1990 typedef UpperTriMatrixView<T,A> uppertri_type; 1991 typedef uppertri_type view_type; 1992 typedef lowertri_type transpose_type; 1993 typedef uppertri_type conjugate_type; 1994 typedef lowertri_type adjoint_type; 1995 typedef MatrixView<T,A> rec_type; 1996 typedef VectorView<T,A> vec_type; 1997 typedef UpperTriMatrixView<RT,A> realpart_type; 1998 typedef ConstLowerTriMatrixView<T,A> const_lowertri_type; 1999 typedef ConstUpperTriMatrixView<T,A> const_uppertri_type; 2000 typedef const_uppertri_type const_view_type; 2001 typedef const_lowertri_type const_transpose_type; 2002 typedef const_uppertri_type const_conjugate_type; 2003 typedef const_lowertri_type const_adjoint_type; 2004 typedef ConstMatrixView<T,A> const_rec_type; 2005 typedef ConstVectorView<T,A> const_vec_type; 2006 typedef ConstUpperTriMatrixView<RT,A> const_realpart_type; 2007 typedef TriRef<T,true> reference; 2008 typedef RMIt<type> rowmajor_iterator; 2009 typedef CMIt<type> colmajor_iterator; 2010 typedef RMIt<const type> const_rowmajor_iterator; 2011 typedef CMIt<const type> const_colmajor_iterator; 2012 2013 // 2014 // Constructors 2015 // 2016 2017 inline UpperTriMatrixView(const type& rhs) : 2018 itsm(rhs.itsm), itss(rhs.itss), itssi(rhs.itssi), itssj(rhs.itssj), 2019 itsdiag(rhs.itsdiag), itsct(rhs.itsct) 2020 TMV_DEFFIRSTLAST(rhs._first,rhs._last) {} 2021 2022 inline UpperTriMatrixView( 2023 T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 2024 DiagType _dt, ConjType _ct 2025 TMV_PARAMFIRSTLAST(T) ) : 2026 itsm(_m), itss(_s), itssi(_si), itssj(_sj), 2027 itsdiag(_dt), itsct(_ct) TMV_DEFFIRSTLAST(_first,_last) {} 2028 2029 virtual inline ~UpperTriMatrixView() 2030 { 2031 TMV_SETFIRSTLAST(0,0); 2032 #ifdef TMV_EXTRA_DEBUG 2033 const_cast<T*&>(itsm) = 0; 2034 #endif 2035 } 2036 2037 // 2038 // Op= 2039 // 2040 2041 inline type& operator=(const type& m2) 2042 { 2043 TMVAssert(size() == m2.size()); 2044 m2.assignToU(*this); 2045 return *this; 2046 } 2047 2048 inline type& operator=(const GenUpperTriMatrix<RT>& m2) 2049 { 2050 TMVAssert(size() == m2.size()); 2051 m2.assignToU(*this); 2052 return *this; 2053 } 2054 2055 inline type& operator=(const GenUpperTriMatrix<CT>& m2) 2056 { 2057 TMVAssert(size() == m2.size()); 2058 TMVAssert(isComplex(T())); 2059 m2.assignToU(*this); 2060 return *this; 2061 } 2062 2063 inline type& operator=(const GenDiagMatrix<RT>& m2) 2064 { 2065 TMVAssert(size() == m2.size()); 2066 m2.assignToD(DiagMatrixViewOf(diag())); 2067 offDiag().setZero(); 2068 return *this; 2069 } 2070 2071 inline type& operator=(const GenDiagMatrix<CT>& m2) 2072 { 2073 TMVAssert(size() == m2.size()); 2074 TMVAssert(isComplex(T())); 2075 m2.assignToD(DiagMatrixViewOf(diag())); 2076 offDiag().setZero(); 2077 return *this; 2078 } 2079 2080 template <typename T2> 2081 inline type& operator=(const GenUpperTriMatrix<T2>& m2) 2082 { 2083 TMVAssert(size() == m2.size()); 2084 TMVAssert(isReal(T2()) || isComplex(T())); 2085 TMVAssert(!isunit() || m2.isunit()); 2086 Copy(m2,*this); 2087 return *this; 2088 } 2089 2090 inline type& operator=(const T& x) 2091 { TMVAssert(!isunit() || x==T(1)); return setToIdentity(x); } 2092 2093 inline type& operator=(const AssignableToUpperTriMatrix<RT>& m2) 2094 { 2095 TMVAssert(size() == m2.size()); 2096 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 2097 m2.assignToU(view()); 2098 return *this; 2099 } 2100 2101 inline type& operator=(const AssignableToUpperTriMatrix<CT>& m2) 2102 { 2103 TMVAssert(size() == m2.size()); 2104 TMVAssert(isComplex(T())); 2105 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 2106 m2.assignToU(view()); 2107 return *this; 2108 } 2109 2110 typedef ListAssigner<T,rowmajor_iterator> MyListAssigner; 2111 inline MyListAssigner operator<<(const T& x) 2112 { return MyListAssigner(rowmajor_begin(),size()*(size()+1)/2,x); } 2113 2114 2115 // 2116 // Access 2117 // 2118 2119 inline reference operator()(ptrdiff_t i,ptrdiff_t j) 2120 { 2121 TMVAssert(i>=0 && i<size()); 2122 TMVAssert(j>=0 && j<size()); 2123 TMVAssert(i<=j); 2124 return ref(i,j); 2125 } 2126 2127 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 2128 { 2129 TMVAssert(i>=0 && i<size()); 2130 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 2131 TMVAssert(j1==j2 || okij(i,j1)); 2132 return vec_type( 2133 ptr()+i*stepi()+j1*stepj(),j2-j1,stepj(),ct() TMV_FIRSTLAST); 2134 } 2135 2136 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 2137 { 2138 TMVAssert(j>=0 && j<size()); 2139 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 2140 TMVAssert(i1==i2 || okij(i2-1,j)); 2141 return vec_type( 2142 ptr()+i1*stepi()+j*stepj(),i2-i1,stepi(),ct() TMV_FIRSTLAST); 2143 } 2144 2145 inline vec_type diag() 2146 { 2147 TMVAssert(!isunit()); 2148 return vec_type( 2149 ptr(),size(),stepi()+stepj(),ct() TMV_FIRSTLAST); 2150 } 2151 2152 inline vec_type diag(ptrdiff_t i) 2153 { 2154 TMVAssert(isunit() ? i>0 : i>=0); 2155 TMVAssert(i<=size()); 2156 return vec_type( 2157 ptr()+i*stepj(),size()-i,(stepi()+stepj()),ct() 2158 TMV_FIRSTLAST); 2159 } 2160 2161 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 2162 { 2163 TMVAssert(isunit() ? i>0 : i>=0); 2164 TMVAssert(i<=size()); 2165 TMVAssert(j1>=0 && j1 <= j2 && j2 <= size()-i); 2166 const ptrdiff_t ds = stepi()+stepj(); 2167 return vec_type( 2168 ptr()+i*stepj()+j1*ds,j2-j1,ds,ct() TMV_FIRSTLAST); 2169 } 2170 2171 // Repeat const versions 2172 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 2173 { return base::operator()(i,j); } 2174 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 2175 { return base::row(i,j1,j2); } 2176 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 2177 { return base::col(j,i1,i2); } 2178 inline const_vec_type diag() const 2179 { return base::diag(); } 2180 inline const_vec_type diag(ptrdiff_t i) const 2181 { return base::diag(i); } 2182 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 2183 { return base::diag(i,j1,j2); } 2184 2185 // 2186 // Modifying Functions 2187 // 2188 2189 type& setZero(); 2190 2191 type& setAllTo(const T& x); 2192 2193 type& addToAll(const T& x); 2194 2195 type& clip(RT thresh); 2196 2197 type& conjugateSelf(); 2198 2199 type& invertSelf(); 2200 2201 type& setToIdentity(const T& x=T(1)); 2202 2203 // 2204 // subMatrix 2205 // 2206 2207 using base::hasSubMatrix; 2208 using base::hasSubVector; 2209 using base::hasSubTriMatrix; 2210 2211 inline rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 2212 { 2213 return rec_type( 2214 ptr()+i1*stepi()+j1*stepj(), 2215 i2-i1, j2-j1, stepi(),stepj(),ct() TMV_FIRSTLAST ); 2216 } 2217 2218 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 2219 { 2220 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 2221 return cSubMatrix(i1,i2,j1,j2); 2222 } 2223 2224 inline rec_type cSubMatrix( 2225 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 2226 { 2227 return rec_type( 2228 ptr()+i1*stepi()+j1*stepj(), 2229 (i2-i1)/istep, (j2-j1)/jstep, istep*stepi(), jstep*stepj(), 2230 ct() TMV_FIRSTLAST ); 2231 } 2232 2233 inline rec_type subMatrix( 2234 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 2235 { 2236 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 2237 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 2238 } 2239 2240 inline vec_type cSubVector( 2241 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 2242 { 2243 return vec_type( 2244 ptr()+i*stepi()+j*stepj(),size, 2245 istep*stepi()+jstep*stepj(),ct() TMV_FIRSTLAST ); 2246 } 2247 2248 inline vec_type subVector( 2249 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 2250 { 2251 TMVAssert(size >= 0); 2252 TMVAssert(hasSubVector(i,j,istep,jstep,size)); 2253 return cSubVector(i,j,istep,jstep,size); 2254 } 2255 2256 inline uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 2257 { 2258 return uppertri_type( 2259 ptr()+i1*(stepi()+stepj()),i2-i1, 2260 stepi(),stepj(),dt(),ct() TMV_FIRSTLAST); 2261 } 2262 2263 inline uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 2264 { 2265 TMVAssert(hasSubTriMatrix(i1,i2,1)); 2266 return cSubTriMatrix(i1,i2); 2267 } 2268 2269 inline uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 2270 { 2271 return uppertri_type( 2272 ptr()+i1*(stepi()+stepj()), 2273 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(),ct() 2274 TMV_FIRSTLAST); 2275 } 2276 2277 inline uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 2278 { 2279 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 2280 return cSubTriMatrix(i1,i2,istep); 2281 } 2282 2283 inline uppertri_type offDiag(ptrdiff_t noff=1) 2284 { 2285 TMVAssert(noff >= 1); 2286 TMVAssert(noff <= size()); 2287 return uppertri_type( 2288 ptr()+noff*stepj(),size()-noff, 2289 stepi(),stepj(),NonUnitDiag,ct() TMV_FIRSTLAST); 2290 } 2291 2292 inline realpart_type realPart() 2293 { 2294 return realpart_type( 2295 reinterpret_cast<RT*>(ptr()), size(), 2296 isReal(T()) ? stepi() : 2*stepi(), 2297 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj 2298 #ifdef TMVFLDEBUG 2299 ,reinterpret_cast<const RT*>(_first) 2300 ,reinterpret_cast<const RT*>(_last) 2301 #endif 2302 ); 2303 } 2304 2305 inline realpart_type imagPart() 2306 { 2307 TMVAssert(isComplex(T())); 2308 TMVAssert(!isunit()); 2309 return realpart_type( 2310 reinterpret_cast<RT*>(ptr())+1, size(), 2311 2*stepi(), 2*stepj(), dt(), NonConj 2312 #ifdef TMVFLDEBUG 2313 ,reinterpret_cast<const RT*>(_first)+1 2314 ,reinterpret_cast<const RT*>(_last)+1 2315 #endif 2316 ); 2317 } 2318 2319 inline uppertri_type view() 2320 { return *this; } 2321 2322 inline uppertri_type viewAsUnitDiag() 2323 { 2324 return uppertri_type( 2325 ptr(),size(),stepi(),stepj(),UnitDiag,ct() TMV_FIRSTLAST); 2326 } 2327 2328 inline lowertri_type transpose() 2329 { 2330 return lowertri_type( 2331 ptr(),size(),stepj(),stepi(),dt(),ct() TMV_FIRSTLAST); 2332 } 2333 2334 inline uppertri_type conjugate() 2335 { 2336 return uppertri_type( 2337 ptr(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,ct()) 2338 TMV_FIRSTLAST); 2339 } 2340 2341 inline lowertri_type adjoint() 2342 { 2343 return lowertri_type( 2344 ptr(),size(),stepj(),stepi(),dt(),TMV_ConjOf(T,ct()) 2345 TMV_FIRSTLAST); 2346 } 2347 2348 2349 inline const_rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 2350 { return base::cSubMatrix(i1,i2,j1,j2); } 2351 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 2352 { return base::subMatrix(i1,i2,j1,j2); } 2353 inline const_rec_type cSubMatrix( 2354 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 2355 { return base::cSubMatrix(i1,i2,j1,j2,istep,jstep); } 2356 inline const_rec_type subMatrix( 2357 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 2358 { return base::subMatrix(i1,i2,j1,j2,istep,jstep); } 2359 inline const_vec_type cSubVector( 2360 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 2361 { return base::cSubVector(i,j,istep,jstep,size); } 2362 inline const_vec_type subVector( 2363 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 2364 { return base::subVector(i,j,istep,jstep,size); } 2365 inline const_uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 2366 { return base::cSubTriMatrix(i1,i2); } 2367 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 2368 { return base::subTriMatrix(i1,i2); } 2369 inline const_uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 2370 { return base::cSubTriMatrix(i1,i2,istep); } 2371 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 2372 { return base::subTriMatrix(i1,i2,istep); } 2373 inline const_uppertri_type offDiag(ptrdiff_t noff=1) const 2374 { return base::offDiag(noff); } 2375 inline const_realpart_type realPart() const 2376 { return base::realPart(); } 2377 inline const_realpart_type imagPart() const 2378 { return base::imagPart(); } 2379 inline const_uppertri_type view() const 2380 { return base::view(); } 2381 inline const_uppertri_type viewAsUnitDiag() const 2382 { return base::viewAsUnitDiag(); } 2383 inline const_lowertri_type transpose() const 2384 { return base::transpose(); } 2385 inline const_uppertri_type conjugate() const 2386 { return base::conjugate(); } 2387 inline const_lowertri_type adjoint() const 2388 { return base::adjoint(); } 2389 2390 // 2391 // I/O 2392 // 2393 2394 void read(const TMV_Reader& reader); 2395 2396 inline ptrdiff_t size() const { return itss; } 2397 inline const T* cptr() const { return itsm; } 2398 inline T* ptr() { return itsm; } 2399 inline ptrdiff_t stepi() const { return itssi; } 2400 inline ptrdiff_t stepj() const { return itssj; } 2401 using base::isconj; 2402 using base::isrm; 2403 using base::iscm; 2404 using base::isunit; 2405 inline DiagType dt() const { return itsdiag; } 2406 inline ConjType ct() const { return itsct; } 2407 2408 inline reference ref(ptrdiff_t i, ptrdiff_t j) 2409 { 2410 T* mi = ptr() + i*stepi() + j*stepj(); 2411 return reference(isunit() && i==j,*mi,ct()); 2412 } 2413 2414 inline rowmajor_iterator rowmajor_begin() 2415 { return rowmajor_iterator(this,0,0); } 2416 inline rowmajor_iterator rowmajor_end() 2417 { return rowmajor_iterator(this,size(),size()); } 2418 2419 inline colmajor_iterator colmajor_begin() 2420 { return colmajor_iterator(this,0,0); } 2421 inline colmajor_iterator colmajor_end() 2422 { return colmajor_iterator(this,0,size()); } 2423 2424 inline const_rowmajor_iterator rowmajor_begin() const 2425 { return const_rowmajor_iterator(this,0,0); } 2426 inline const_rowmajor_iterator rowmajor_end() const 2427 { return const_rowmajor_iterator(this,size(),size()); } 2428 2429 inline const_colmajor_iterator colmajor_begin() const 2430 { return const_colmajor_iterator(this,0,0); } 2431 inline const_colmajor_iterator colmajor_end() const 2432 { return const_colmajor_iterator(this,0,size()); } 2433 2434 2435 protected : 2436 2437 T*const itsm; 2438 const ptrdiff_t itss; 2439 const ptrdiff_t itssi; 2440 const ptrdiff_t itssj; 2441 DiagType itsdiag; 2442 ConjType itsct; 2443 2444 using base::okij; 2445 2446 #ifdef TMVFLDEBUG 2447 public: 2448 const T* _first; 2449 const T* _last; 2450 #endif 2451 2452 }; // UpperTriMatrixView 2453 2454 template <typename T, int A> 2455 class LowerTriMatrixView : public GenLowerTriMatrix<T> 2456 { 2457 public: 2458 2459 typedef TMV_RealType(T) RT; 2460 typedef TMV_ComplexType(T) CT; 2461 typedef GenLowerTriMatrix<T> base; 2462 typedef LowerTriMatrixView<T,A> type; 2463 typedef UpperTriMatrixView<T,A> uppertri_type; 2464 typedef LowerTriMatrixView<T,A> lowertri_type; 2465 typedef lowertri_type view_type; 2466 typedef uppertri_type transpose_type; 2467 typedef lowertri_type conjugate_type; 2468 typedef uppertri_type adjoint_type; 2469 typedef MatrixView<T,A> rec_type; 2470 typedef VectorView<T,A> vec_type; 2471 typedef LowerTriMatrixView<RT,A> realpart_type; 2472 typedef ConstUpperTriMatrixView<T,A> const_uppertri_type; 2473 typedef ConstLowerTriMatrixView<T,A> const_lowertri_type; 2474 typedef const_lowertri_type const_view_type; 2475 typedef const_uppertri_type const_transpose_type; 2476 typedef const_lowertri_type const_conjugate_type; 2477 typedef const_uppertri_type const_adjoint_type; 2478 typedef ConstMatrixView<T,A> const_rec_type; 2479 typedef ConstVectorView<T,A> const_vec_type; 2480 typedef ConstLowerTriMatrixView<RT,A> const_realpart_type; 2481 typedef TriRef<T,true> reference; 2482 typedef RMIt<type> rowmajor_iterator; 2483 typedef CMIt<type> colmajor_iterator; 2484 typedef RMIt<const type> const_rowmajor_iterator; 2485 typedef CMIt<const type> const_colmajor_iterator; 2486 2487 // 2488 // Constructors 2489 // 2490 2491 inline LowerTriMatrixView(const type& rhs) : 2492 itsm(rhs.itsm), itss(rhs.itss), itssi(rhs.itssi), itssj(rhs.itssj), 2493 itsdiag(rhs.itsdiag), itsct(rhs.itsct) 2494 TMV_DEFFIRSTLAST(rhs._first,rhs._last) {} 2495 2496 inline LowerTriMatrixView( 2497 T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, DiagType _dt, ConjType _ct 2498 TMV_PARAMFIRSTLAST(T) ) : 2499 itsm(_m), itss(_s), itssi(_si), itssj(_sj), 2500 itsdiag(_dt), itsct(_ct) TMV_DEFFIRSTLAST(_first,_last) {} 2501 2502 virtual inline ~LowerTriMatrixView() 2503 { 2504 TMV_SETFIRSTLAST(0,0); 2505 #ifdef TMV_EXTRA_DEBUG 2506 const_cast<T*&>(itsm) = 0; 2507 #endif 2508 } 2509 2510 // 2511 // Op= 2512 // 2513 2514 inline type& operator=(const type& m2) 2515 { 2516 TMVAssert(size() == m2.size()); 2517 m2.assignToL(*this); 2518 return *this; 2519 } 2520 2521 inline type& operator=(const GenLowerTriMatrix<RT>& m2) 2522 { 2523 TMVAssert(size() == m2.size()); 2524 m2.assignToL(*this); 2525 return *this; 2526 } 2527 2528 inline type& operator=(const GenLowerTriMatrix<CT>& m2) 2529 { 2530 TMVAssert(size() == m2.size()); 2531 TMVAssert(isComplex(T())); 2532 m2.assignToL(*this); 2533 return *this; 2534 } 2535 2536 inline type& operator=(const GenDiagMatrix<RT>& m2) 2537 { 2538 TMVAssert(size() == m2.size()); 2539 transpose() = m2; 2540 return *this; 2541 } 2542 2543 inline type& operator=(const GenDiagMatrix<CT>& m2) 2544 { 2545 TMVAssert(size() == m2.size()); 2546 TMVAssert(isComplex(T())); 2547 transpose() = m2; 2548 return *this; 2549 } 2550 2551 template <typename T2> 2552 inline type& operator=(const GenLowerTriMatrix<T2>& m2) 2553 { 2554 TMVAssert(size() == m2.size()); 2555 TMVAssert(!isunit() || m2.isunit()); 2556 Copy(m2.transpose(),transpose()); 2557 return *this; 2558 } 2559 2560 inline type& operator=(const T& x) 2561 { 2562 TMVAssert(!isunit()); 2563 return setToIdentity(x); 2564 } 2565 2566 inline type& operator=(const AssignableToLowerTriMatrix<RT>& m2) 2567 { 2568 TMVAssert(size() == m2.size()); 2569 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 2570 m2.assignToL(view()); 2571 return *this; 2572 } 2573 2574 inline type& operator=(const AssignableToLowerTriMatrix<CT>& m2) 2575 { 2576 TMVAssert(size() == m2.size()); 2577 TMVAssert(isComplex(T())); 2578 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 2579 m2.assignToL(view()); 2580 return *this; 2581 } 2582 2583 typedef ListAssigner<T,rowmajor_iterator> MyListAssigner; 2584 inline MyListAssigner operator<<(const T& x) 2585 { return MyListAssigner(rowmajor_begin(),size()*(size()+1)/2,x); } 2586 2587 // 2588 // Access 2589 // 2590 2591 inline reference operator()(ptrdiff_t i,ptrdiff_t j) 2592 { 2593 TMVAssert(i>=0 && i<size()); 2594 TMVAssert(j>=0 && j<size()); 2595 TMVAssert(i>=j); 2596 return ref(i,j); 2597 } 2598 2599 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 2600 { 2601 TMVAssert(i>=0 && i<size()); 2602 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 2603 TMVAssert(j1==j2 || okij(i,j2-1)); 2604 return vec_type( 2605 ptr()+i*stepi()+j1*stepj(),j2-j1,stepj(),ct() TMV_FIRSTLAST); 2606 } 2607 2608 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 2609 { 2610 TMVAssert(j>=0 && j<size()); 2611 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 2612 TMVAssert(i1==i2 || okij(i1,j)); 2613 return vec_type( 2614 ptr()+i1*stepi()+j*stepj(),i2-i1,stepi(),ct() TMV_FIRSTLAST); 2615 } 2616 2617 inline vec_type diag() 2618 { 2619 TMVAssert(!isunit()); 2620 return vec_type( 2621 ptr(),size(),stepi()+stepj(),ct() TMV_FIRSTLAST); 2622 } 2623 2624 inline vec_type diag(ptrdiff_t i) 2625 { 2626 TMVAssert(i>=-size()); 2627 TMVAssert(isunit() ? i<0 : i<=0); 2628 return vec_type( 2629 ptr()-i*stepi(),size()+i,stepi()+stepj(),ct() TMV_FIRSTLAST); 2630 } 2631 2632 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 2633 { 2634 TMVAssert(i>=-size()); 2635 TMVAssert(isunit() ? i<0 : i<=0); 2636 TMVAssert(j1>=0 && j1 <= j2 && j2 <= size()+i); 2637 const ptrdiff_t ds = stepi()+stepj(); 2638 return vec_type( 2639 ptr()-i*stepi()+j1*ds,j2-j1,ds,ct() TMV_FIRSTLAST); 2640 } 2641 2642 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 2643 { return base::operator()(i,j); } 2644 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 2645 { return base::row(i,j1,j2); } 2646 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 2647 { return base::col(j,i1,i2); } 2648 inline const_vec_type diag() const 2649 { return base::diag(); } 2650 inline const_vec_type diag(ptrdiff_t i) const 2651 { return base::diag(i); } 2652 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 2653 { return base::diag(i,j1,j2); } 2654 2655 // 2656 // Modifying Functions 2657 // 2658 2659 inline type& setZero() 2660 { transpose().setZero(); return *this; } 2661 2662 inline type& setAllTo(const T& x) 2663 { transpose().setAllTo(x); return *this; } 2664 2665 inline type& addToAll(const T& x) 2666 { transpose().addToAll(x); return *this; } 2667 2668 inline type& clip(RT thresh) 2669 { transpose().clip(thresh); return *this; } 2670 2671 inline type& conjugateSelf() 2672 { transpose().conjugateSelf(); return *this; } 2673 2674 inline type& invertSelf() 2675 { transpose().invertSelf(); return *this; } 2676 2677 inline type& setToIdentity(const T& x=T(1)) 2678 { transpose().setToIdentity(x); return *this; } 2679 2680 // 2681 // subMatrix 2682 // 2683 2684 using base::hasSubMatrix; 2685 using base::hasSubVector; 2686 using base::hasSubTriMatrix; 2687 2688 inline rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 2689 { 2690 return rec_type( 2691 ptr()+i1*stepi()+j1*stepj(), 2692 i2-i1, j2-j1, stepi(),stepj(),ct() TMV_FIRSTLAST ); 2693 } 2694 2695 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 2696 { 2697 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 2698 return cSubMatrix(i1,i2,j1,j2); 2699 } 2700 2701 inline rec_type cSubMatrix( 2702 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 2703 { 2704 return rec_type( 2705 ptr()+i1*stepi()+j1*stepj(), 2706 (i2-i1)/istep, (j2-j1)/jstep, istep*stepi(), jstep*stepj(), 2707 ct() TMV_FIRSTLAST ); 2708 } 2709 2710 inline rec_type subMatrix( 2711 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 2712 { 2713 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 2714 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 2715 } 2716 2717 inline vec_type cSubVector( 2718 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 2719 { 2720 return vec_type( 2721 ptr()+i*stepi()+j*stepj(),size, 2722 istep*stepi()+jstep*stepj(),ct() TMV_FIRSTLAST ); 2723 } 2724 2725 inline vec_type subVector( 2726 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 2727 { 2728 TMVAssert(size >= 0); 2729 TMVAssert(hasSubVector(i,j,istep,jstep,size)); 2730 return cSubVector(i,j,istep,jstep,size); 2731 } 2732 2733 inline lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 2734 { 2735 return lowertri_type( 2736 ptr()+i1*(stepi()+stepj()), 2737 i2-i1,stepi(),stepj(),dt(),ct() TMV_FIRSTLAST); 2738 } 2739 2740 inline lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 2741 { 2742 TMVAssert(hasSubTriMatrix(i1,i2,1)); 2743 return cSubTriMatrix(i1,i2); 2744 } 2745 2746 inline lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 2747 { 2748 return lowertri_type( 2749 ptr()+i1*(stepi()+stepj()), 2750 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(),ct() 2751 TMV_FIRSTLAST); 2752 } 2753 2754 inline lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 2755 { 2756 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 2757 return cSubTriMatrix(i1,i2,istep); 2758 } 2759 2760 inline lowertri_type offDiag(ptrdiff_t noff=1) 2761 { 2762 TMVAssert(noff >= 1); 2763 TMVAssert(noff <= size()); 2764 return lowertri_type( 2765 ptr()+noff*stepi(),size()-noff, 2766 stepi(),stepj(),NonUnitDiag,ct() TMV_FIRSTLAST); 2767 } 2768 2769 inline realpart_type realPart() 2770 { 2771 return realpart_type( 2772 reinterpret_cast<RT*>(ptr()), size(), 2773 isReal(T()) ? stepi() : 2*stepi(), 2774 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj 2775 #ifdef TMVFLDEBUG 2776 ,reinterpret_cast<const RT*>(_first) 2777 ,reinterpret_cast<const RT*>(_last) 2778 #endif 2779 ); 2780 } 2781 2782 inline realpart_type imagPart() 2783 { 2784 TMVAssert(isComplex(T())); 2785 TMVAssert(!isunit()); 2786 return realpart_type( 2787 reinterpret_cast<RT*>(ptr())+1, size(), 2788 2*stepi(), 2*stepj(), dt(), NonConj 2789 #ifdef TMVFLDEBUG 2790 ,reinterpret_cast<const RT*>(_first)+1 2791 ,reinterpret_cast<const RT*>(_last)+1 2792 #endif 2793 ); 2794 } 2795 2796 inline lowertri_type view() 2797 { return *this; } 2798 2799 inline lowertri_type viewAsUnitDiag() 2800 { 2801 return lowertri_type( 2802 ptr(),size(), 2803 stepi(),stepj(),UnitDiag,ct() TMV_FIRSTLAST); 2804 } 2805 2806 inline uppertri_type transpose() 2807 { 2808 return uppertri_type( 2809 ptr(),size(), 2810 stepj(),stepi(),dt(),ct() TMV_FIRSTLAST); 2811 } 2812 2813 inline lowertri_type conjugate() 2814 { 2815 return lowertri_type( 2816 ptr(),size(), 2817 stepi(),stepj(),dt(),TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 2818 } 2819 2820 inline uppertri_type adjoint() 2821 { 2822 return uppertri_type( 2823 ptr(),size(), 2824 stepj(),stepi(),dt(),TMV_ConjOf(T,ct()) 2825 TMV_FIRSTLAST); 2826 } 2827 2828 inline const_rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 2829 { return base::cSubMatrix(i1,i2,j1,j2); } 2830 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 2831 { return base::subMatrix(i1,i2,j1,j2); } 2832 inline const_rec_type cSubMatrix( 2833 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 2834 { return base::cSubMatrix(i1,i2,j1,j2,istep,jstep); } 2835 inline const_rec_type subMatrix( 2836 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 2837 { return base::subMatrix(i1,i2,j1,j2,istep,jstep); } 2838 inline const_vec_type cSubVector( 2839 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 2840 { return base::cSubVector(i,j,istep,jstep,size); } 2841 inline const_vec_type subVector( 2842 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 2843 { return base::subVector(i,j,istep,jstep,size); } 2844 inline const_lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 2845 { return base::cSubTriMatrix(i1,i2); } 2846 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 2847 { return base::subTriMatrix(i1,i2); } 2848 inline const_lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 2849 { return base::cSubTriMatrix(i1,i2,istep); } 2850 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 2851 { return base::subTriMatrix(i1,i2,istep); } 2852 inline const_lowertri_type offDiag(ptrdiff_t noff=1) const 2853 { return base::offDiag(noff); } 2854 inline const_realpart_type realPart() const 2855 { return base::realPart(); } 2856 inline const_realpart_type imagPart() const 2857 { return base::imagPart(); } 2858 inline const_lowertri_type view() const 2859 { return base::view(); } 2860 inline const_lowertri_type viewAsUnitDiag() const 2861 { return base::viewAsUnitDiag(); } 2862 inline const_uppertri_type transpose() const 2863 { return base::transpose(); } 2864 inline const_lowertri_type conjugate() const 2865 { return base::conjugate(); } 2866 inline const_uppertri_type adjoint() const 2867 { return base::adjoint(); } 2868 2869 // 2870 // I/O 2871 // 2872 2873 void read(const TMV_Reader& reader); 2874 2875 inline ptrdiff_t size() const { return itss; } 2876 inline const T* cptr() const { return itsm; } 2877 inline T* ptr() { return itsm; } 2878 inline ptrdiff_t stepi() const { return itssi; } 2879 inline ptrdiff_t stepj() const { return itssj; } 2880 using base::isconj; 2881 using base::isrm; 2882 using base::iscm; 2883 using base::isunit; 2884 inline DiagType dt() const { return itsdiag; } 2885 inline ConjType ct() const { return itsct; } 2886 2887 inline reference ref(ptrdiff_t i, ptrdiff_t j) 2888 { 2889 T* mi = ptr() + i*stepi() + j*stepj(); 2890 return reference(isunit() && i==j,*mi,ct()); 2891 } 2892 2893 inline rowmajor_iterator rowmajor_begin() 2894 { return rowmajor_iterator(this,0,0); } 2895 inline rowmajor_iterator rowmajor_end() 2896 { return rowmajor_iterator(this,size(),0); } 2897 2898 inline colmajor_iterator colmajor_begin() 2899 { return colmajor_iterator(this,0,0); } 2900 inline colmajor_iterator colmajor_end() 2901 { return colmajor_iterator(this,size(),size()); } 2902 2903 inline const_rowmajor_iterator rowmajor_begin() const 2904 { return const_rowmajor_iterator(this,0,0); } 2905 inline const_rowmajor_iterator rowmajor_end() const 2906 { return const_rowmajor_iterator(this,size(),0); } 2907 2908 inline const_colmajor_iterator colmajor_begin() const 2909 { return const_colmajor_iterator(this,0,0); } 2910 inline const_colmajor_iterator colmajor_end() const 2911 { return const_colmajor_iterator(this,size(),size()); } 2912 2913 2914 protected : 2915 2916 T*const itsm; 2917 const ptrdiff_t itss; 2918 const ptrdiff_t itssi; 2919 const ptrdiff_t itssj; 2920 DiagType itsdiag; 2921 ConjType itsct; 2922 2923 using base::okij; 2924 2925 #ifdef TMVFLDEBUG 2926 public: 2927 const T* _first; 2928 const T* _last; 2929 #endif 2930 2931 }; // LowerTriMatrixView 2932 2933 template <typename T> 2934 class UpperTriMatrixView<T,FortranStyle> : 2935 public UpperTriMatrixView<T,CStyle> 2936 { 2937 public: 2938 2939 typedef TMV_RealType(T) RT; 2940 typedef TMV_ComplexType(T) CT; 2941 typedef GenUpperTriMatrix<T> base; 2942 typedef UpperTriMatrixView<T,FortranStyle> type; 2943 typedef UpperTriMatrixView<T,CStyle> c_type; 2944 typedef ConstUpperTriMatrixView<T,FortranStyle> const_type; 2945 typedef LowerTriMatrixView<T,FortranStyle> lowertri_type; 2946 typedef UpperTriMatrixView<T,FortranStyle> uppertri_type; 2947 typedef uppertri_type view_type; 2948 typedef lowertri_type transpose_type; 2949 typedef uppertri_type conjugate_type; 2950 typedef lowertri_type adjoint_type; 2951 typedef MatrixView<T,FortranStyle> rec_type; 2952 typedef VectorView<T,FortranStyle> vec_type; 2953 typedef UpperTriMatrixView<RT,FortranStyle> realpart_type; 2954 typedef ConstLowerTriMatrixView<T,FortranStyle> const_lowertri_type; 2955 typedef ConstUpperTriMatrixView<T,FortranStyle> const_uppertri_type; 2956 typedef const_uppertri_type const_view_type; 2957 typedef const_lowertri_type const_transpose_type; 2958 typedef const_uppertri_type const_conjugate_type; 2959 typedef const_lowertri_type const_adjoint_type; 2960 typedef ConstMatrixView<T,FortranStyle> const_rec_type; 2961 typedef ConstVectorView<T,FortranStyle> const_vec_type; 2962 typedef ConstUpperTriMatrixView<RT,FortranStyle> const_realpart_type; 2963 typedef TriRef<T,true> reference; 2964 2965 // 2966 // Constructors 2967 // 2968 2969 inline UpperTriMatrixView(const type& rhs) : c_type(rhs) {} 2970 2971 inline UpperTriMatrixView(const c_type& rhs) : c_type(rhs) {} 2972 2973 inline UpperTriMatrixView( 2974 T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, DiagType _dt, ConjType _ct 2975 TMV_PARAMFIRSTLAST(T) ) : 2976 c_type(_m,_s,_si,_sj,_dt,_ct TMV_FIRSTLAST1(_first,_last) ) 2977 {} 2978 2979 virtual inline ~UpperTriMatrixView() {} 2980 2981 // 2982 // Op= 2983 // 2984 2985 inline type& operator=(const type& m2) 2986 { c_type::operator=(m2); return *this; } 2987 2988 inline type& operator=(const c_type& m2) 2989 { c_type::operator=(m2); return *this; } 2990 2991 inline type& operator=(const GenUpperTriMatrix<RT>& m2) 2992 { c_type::operator=(m2); return *this; } 2993 2994 inline type& operator=(const GenUpperTriMatrix<CT>& m2) 2995 { c_type::operator=(m2); return *this; } 2996 2997 inline type& operator=(const GenDiagMatrix<RT>& m2) 2998 { c_type::operator=(m2); return *this; } 2999 3000 inline type& operator=(const GenDiagMatrix<CT>& m2) 3001 { c_type::operator=(m2); return *this; } 3002 3003 template <typename T2> 3004 inline type& operator=(const GenUpperTriMatrix<T2>& m2) 3005 { c_type::operator=(m2); return *this; } 3006 3007 inline type& operator=(const T& x) 3008 { c_type::operator=(x); return *this; } 3009 3010 inline type& operator=(const AssignableToUpperTriMatrix<RT>& m2) 3011 { c_type::operator=(m2); return *this; } 3012 3013 inline type& operator=(const AssignableToUpperTriMatrix<CT>& m2) 3014 { c_type::operator=(m2); return *this; } 3015 3016 typedef typename c_type::MyListAssigner MyListAssigner; 3017 inline MyListAssigner operator<<(const T& x) 3018 { return c_type::operator<<(x); } 3019 3020 // 3021 // Access 3022 // 3023 3024 inline reference operator()(ptrdiff_t i,ptrdiff_t j) 3025 { 3026 TMVAssert(i>0 && i <= c_type::size()); 3027 TMVAssert(j>0 && j <= c_type::size()); 3028 TMVAssert(i<=j); 3029 return c_type::ref(i-1,j-1); 3030 } 3031 3032 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3033 { 3034 TMVAssert(i>0 && i<=c_type::size()); 3035 TMVAssert(j1>0 && j1<=j2 && j2<=c_type::size()); 3036 return c_type::row(i-1,j1-1,j2); 3037 } 3038 3039 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 3040 { 3041 TMVAssert(j>0 && j<=c_type::size()); 3042 TMVAssert(i1>0 && i1<=i2 && i2<=c_type::size()); 3043 return c_type::col(j-1,i1-1,i2); 3044 } 3045 3046 inline vec_type diag() 3047 { return c_type::diag(); } 3048 3049 inline vec_type diag(ptrdiff_t i) 3050 { return c_type::diag(i); } 3051 3052 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3053 { 3054 TMVAssert(j1>0); 3055 return c_type::diag(i,j1-1,j2); 3056 } 3057 3058 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 3059 { 3060 TMVAssert(i>0 && i <= c_type::size()); 3061 TMVAssert(j>0 && j <= c_type::size()); 3062 TMVAssert(i<=j); 3063 return c_type::cref(i-1,j-1); 3064 } 3065 3066 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3067 { 3068 TMVAssert(i>0 && i<=c_type::size()); 3069 TMVAssert(j1>0 && j1<=j2 && j2<=c_type::size()); 3070 return c_type::row(i-1,j1-1,j2); 3071 } 3072 3073 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 3074 { 3075 TMVAssert(j>0 && j<=c_type::size()); 3076 TMVAssert(i1>0 && i1<=i2 && i2<=c_type::size()); 3077 return c_type::col(j-1,i1-1,i2); 3078 } 3079 3080 inline const_vec_type diag() const 3081 { return c_type::diag(); } 3082 3083 inline const_vec_type diag(ptrdiff_t i) const 3084 { return c_type::diag(i); } 3085 3086 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3087 { 3088 TMVAssert(j1>0); 3089 return c_type::diag(i,j1-1,j2); 3090 } 3091 3092 // 3093 // Modifying Functions 3094 // 3095 3096 inline type& setZero() 3097 { c_type::setZero(); return *this; } 3098 3099 inline type& setAllTo(const T& x) 3100 { c_type::setAllTo(x); return *this; } 3101 3102 inline type& addToAll(const T& x) 3103 { c_type::addToAll(x); return *this; } 3104 3105 inline type& clip(RT thresh) 3106 { c_type::clip(thresh); return *this; } 3107 3108 inline type& conjugateSelf() 3109 { c_type::conjugateSelf(); return *this; } 3110 3111 inline type& invertSelf() 3112 { c_type::invertSelf(); return *this; } 3113 3114 inline type& setToIdentity(const T& x=T(1)) 3115 { c_type::setToIdentity(x); return *this; } 3116 3117 // 3118 // SubMatrix 3119 // 3120 3121 inline bool hasSubMatrix( 3122 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 3123 { return const_type(*this).hasSubMatrix(i1,i2,j1,j2,istep,jstep); } 3124 3125 inline bool hasSubVector( 3126 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t s) const 3127 { return const_type(*this).hasSubVector(i,j,istep,jstep,s); } 3128 3129 inline bool hasSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 3130 { return const_type(*this).hasSubTriMatrix(i1,i2,istep); } 3131 3132 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 3133 { 3134 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 3135 return c_type::cSubMatrix(i1-1,i2,j1-1,j2); 3136 } 3137 3138 inline rec_type subMatrix( 3139 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 3140 { 3141 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 3142 return c_type::cSubMatrix( 3143 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 3144 } 3145 3146 inline vec_type subVector( 3147 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t s) 3148 { 3149 TMVAssert(hasSubVector(i,j,istep,jstep,s)); 3150 return c_type::cSubVector(i-1,j-1,istep,jstep,s); 3151 } 3152 3153 inline uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 3154 { 3155 TMVAssert(hasSubTriMatrix(i1,i2,1)); 3156 return c_type::cSubTriMatrix(i1-1,i2); 3157 } 3158 3159 inline uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 3160 { 3161 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 3162 return c_type::cSubTriMatrix(i1-1,i2-1+istep,istep); 3163 } 3164 3165 inline uppertri_type offDiag(ptrdiff_t noff=1) 3166 { 3167 TMVAssert(noff >= 1); 3168 TMVAssert(noff <= c_type::size()); 3169 return c_type::offDiag(noff); 3170 } 3171 3172 inline realpart_type realPart() 3173 { return c_type::realPart(); } 3174 3175 inline realpart_type imagPart() 3176 { return c_type::imagPart(); } 3177 3178 inline uppertri_type view() 3179 { return *this; } 3180 3181 inline uppertri_type viewAsUnitDiag() 3182 { return c_type::viewAsUnitDiag(); } 3183 3184 inline lowertri_type transpose() 3185 { return c_type::transpose(); } 3186 3187 inline uppertri_type conjugate() 3188 { return c_type::conjugate(); } 3189 3190 inline lowertri_type adjoint() 3191 { return c_type::adjoint(); } 3192 3193 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 3194 { 3195 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 3196 return c_type::cSubMatrix(i1-1,i2,j1-1,j2); 3197 } 3198 3199 inline const_rec_type subMatrix( 3200 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 3201 { 3202 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 3203 return c_type::cSubMatrix( 3204 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 3205 } 3206 3207 inline const_vec_type subVector( 3208 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t s) const 3209 { 3210 TMVAssert(hasSubVector(i,j,istep,jstep,s)); 3211 return c_type::cSubVector(i-1,j-1,istep,jstep,s); 3212 } 3213 3214 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 3215 { 3216 TMVAssert(hasSubTriMatrix(i1,i2,1)); 3217 return c_type::cSubTriMatrix(i1-1,i2); 3218 } 3219 3220 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 3221 { 3222 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 3223 return c_type::cSubTriMatrix(i1-1,i2-1+istep,istep); 3224 } 3225 3226 inline const_uppertri_type offDiag(ptrdiff_t noff=1) const 3227 { 3228 TMVAssert(noff >= 1); 3229 TMVAssert(noff <= c_type::size()); 3230 return c_type::offDiag(noff); 3231 } 3232 3233 inline const_realpart_type realPart() const 3234 { return c_type::realPart(); } 3235 3236 inline const_realpart_type imagPart() const 3237 { return c_type::imagPart(); } 3238 3239 inline const_uppertri_type view() const 3240 { return *this; } 3241 3242 inline const_uppertri_type viewAsUnitDiag() const 3243 { return c_type::viewAsUnitDiag(); } 3244 3245 inline const_lowertri_type transpose() const 3246 { return c_type::transpose(); } 3247 3248 inline const_uppertri_type conjugate() const 3249 { return c_type::conjugate(); } 3250 3251 inline const_lowertri_type adjoint() const 3252 { return c_type::adjoint(); } 3253 3254 }; // FortranStyle UpperTriMatrixView 3255 3256 template <typename T> 3257 class LowerTriMatrixView<T,FortranStyle> : 3258 public LowerTriMatrixView<T,CStyle> 3259 { 3260 public: 3261 3262 typedef TMV_RealType(T) RT; 3263 typedef TMV_ComplexType(T) CT; 3264 typedef GenLowerTriMatrix<T> base; 3265 typedef LowerTriMatrixView<T,FortranStyle> type; 3266 typedef LowerTriMatrixView<T,CStyle> c_type; 3267 typedef ConstLowerTriMatrixView<T,FortranStyle> const_type; 3268 typedef UpperTriMatrixView<T,FortranStyle> uppertri_type; 3269 typedef LowerTriMatrixView<T,FortranStyle> lowertri_type; 3270 typedef lowertri_type view_type; 3271 typedef uppertri_type transpose_type; 3272 typedef lowertri_type conjugate_type; 3273 typedef uppertri_type adjoint_type; 3274 typedef MatrixView<T,FortranStyle> rec_type; 3275 typedef VectorView<T,FortranStyle> vec_type; 3276 typedef LowerTriMatrixView<RT,FortranStyle> realpart_type; 3277 typedef ConstUpperTriMatrixView<T,FortranStyle> const_uppertri_type; 3278 typedef ConstLowerTriMatrixView<T,FortranStyle> const_lowertri_type; 3279 typedef const_lowertri_type const_view_type; 3280 typedef const_uppertri_type const_transpose_type; 3281 typedef const_lowertri_type const_conjugate_type; 3282 typedef const_uppertri_type const_adjoint_type; 3283 typedef ConstMatrixView<T,FortranStyle> const_rec_type; 3284 typedef ConstVectorView<T,FortranStyle> const_vec_type; 3285 typedef ConstLowerTriMatrixView<RT,FortranStyle> const_realpart_type; 3286 typedef TriRef<T,true> reference; 3287 3288 // 3289 // Constructors 3290 // 3291 3292 inline LowerTriMatrixView(const type& rhs) : c_type(rhs) {} 3293 3294 inline LowerTriMatrixView(const c_type& rhs) : c_type(rhs) {} 3295 3296 inline LowerTriMatrixView( 3297 T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, DiagType _dt, ConjType _ct 3298 TMV_PARAMFIRSTLAST(T) ) : 3299 c_type(_m,_s,_si,_sj,_dt,_ct TMV_FIRSTLAST1(_first,_last) ) 3300 {} 3301 3302 virtual inline ~LowerTriMatrixView() {} 3303 3304 // 3305 // Op= 3306 // 3307 3308 inline type& operator=(const type& m2) 3309 { c_type::operator=(m2); return *this; } 3310 3311 inline type& operator=(const c_type& m2) 3312 { c_type::operator=(m2); return *this; } 3313 3314 inline type& operator=(const GenLowerTriMatrix<RT>& m2) 3315 { c_type::operator=(m2); return *this; } 3316 3317 inline type& operator=(const GenLowerTriMatrix<CT>& m2) 3318 { c_type::operator=(m2); return *this; } 3319 3320 inline type& operator=(const GenDiagMatrix<RT>& m2) 3321 { c_type::operator=(m2); return *this; } 3322 3323 inline type& operator=(const GenDiagMatrix<CT>& m2) 3324 { c_type::operator=(m2); return *this; } 3325 3326 template <typename T2> 3327 inline type& operator=(const GenLowerTriMatrix<T2>& m2) 3328 { c_type::operator=(m2); return *this; } 3329 3330 inline type& operator=(const T& x) 3331 { c_type::operator=(x); return *this; } 3332 3333 inline type& operator=(const AssignableToLowerTriMatrix<RT>& m2) 3334 { c_type::operator=(m2); return *this; } 3335 3336 inline type& operator=(const AssignableToLowerTriMatrix<CT>& m2) 3337 { c_type::operator=(m2); return *this; } 3338 3339 typedef typename c_type::MyListAssigner MyListAssigner; 3340 inline MyListAssigner operator<<(const T& x) 3341 { return c_type::operator<<(x); } 3342 3343 // 3344 // Access 3345 // 3346 3347 inline reference operator()(ptrdiff_t i,ptrdiff_t j) 3348 { 3349 TMVAssert(i>0 && i <= c_type::size()); 3350 TMVAssert(j>0 && j <= c_type::size()); 3351 TMVAssert(i>=j); 3352 return c_type::ref(i-1,j-1); 3353 } 3354 3355 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3356 { 3357 TMVAssert(i>0 && i<=c_type::size()); 3358 TMVAssert(j1>0 && j1<=j2 && j2<=c_type::size()); 3359 return c_type::row(i-1,j1-1,j2); 3360 } 3361 3362 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 3363 { 3364 TMVAssert(j>0 && j<=c_type::size()); 3365 TMVAssert(i1>0 && i1<=i2 && i2<=c_type::size()); 3366 return c_type::col(j-1,i1-1,i2); 3367 } 3368 3369 inline vec_type diag() 3370 { return c_type::diag(); } 3371 3372 inline vec_type diag(ptrdiff_t i) 3373 { return c_type::diag(i); } 3374 3375 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3376 { 3377 TMVAssert(j1>0); 3378 return c_type::diag(i,j1-1,j2); 3379 } 3380 3381 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 3382 { 3383 TMVAssert(i>0 && i <= c_type::size()); 3384 TMVAssert(j>0 && j <= c_type::size()); 3385 TMVAssert(i>=j); 3386 return c_type::ref(i-1,j-1); 3387 } 3388 3389 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3390 { 3391 TMVAssert(i>0 && i<=c_type::size()); 3392 TMVAssert(j1>0 && j1<=j2 && j2<=c_type::size()); 3393 return c_type::row(i-1,j1-1,j2); 3394 } 3395 3396 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 3397 { 3398 TMVAssert(j>0 && j<=c_type::size()); 3399 TMVAssert(i1>0 && i1<=i2 && i2<=c_type::size()); 3400 return c_type::col(j-1,i1-1,i2); 3401 } 3402 3403 inline const_vec_type diag() const 3404 { return c_type::diag(); } 3405 3406 inline const_vec_type diag(ptrdiff_t i) const 3407 { return c_type::diag(i); } 3408 3409 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3410 { 3411 TMVAssert(j1>0); 3412 return c_type::diag(i,j1-1,j2); 3413 } 3414 3415 // 3416 // Modifying Functions 3417 // 3418 3419 inline type& setZero() 3420 { c_type::setZero(); return *this; } 3421 3422 inline type& setAllTo(const T& x) 3423 { c_type::setAllTo(x); return *this; } 3424 3425 inline type& addToAll(const T& x) 3426 { c_type::addToAll(x); return *this; } 3427 3428 inline type& clip(RT thresh) 3429 { c_type::clip(thresh); return *this; } 3430 3431 inline type& conjugateSelf() 3432 { c_type::conjugateSelf(); return *this; } 3433 3434 inline type& invertSelf() 3435 { c_type::invertSelf(); return *this; } 3436 3437 inline type& setToIdentity(const T& x=T(1)) 3438 { c_type::setToIdentity(x); return *this; } 3439 3440 // 3441 // subMatrix 3442 // 3443 3444 inline bool hasSubMatrix( 3445 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 3446 { return const_type(*this).hasSubMatrix(i1,i2,j1,j2,istep,jstep); } 3447 3448 inline bool hasSubVector( 3449 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t s) const 3450 { return const_type(*this).hasSubVector(i,j,istep,jstep,s); } 3451 3452 inline bool hasSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 3453 { return const_type(*this).hasSubTriMatrix(i1,i2,istep); } 3454 3455 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 3456 { 3457 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 3458 return c_type::cSubMatrix(i1-1,i2,j1-1,j2); 3459 } 3460 3461 inline rec_type subMatrix( 3462 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 3463 { 3464 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 3465 return c_type::cSubMatrix( 3466 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 3467 } 3468 3469 inline vec_type subVector( 3470 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t s) 3471 { 3472 TMVAssert(hasSubVector(i,j,istep,jstep,s)); 3473 return c_type::cSubVector(i-1,j-1,istep,jstep,s); 3474 } 3475 3476 inline lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 3477 { 3478 TMVAssert(hasSubTriMatrix(i1,i2,1)); 3479 return c_type::cSubTriMatrix(i1-1,i2); 3480 } 3481 3482 inline lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 3483 { 3484 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 3485 return c_type::cSubTriMatrix(i1-1,i2-1+istep,istep); 3486 } 3487 3488 inline lowertri_type offDiag(ptrdiff_t noff=1) 3489 { 3490 TMVAssert(noff >= 1); 3491 TMVAssert(noff <= c_type::size()); 3492 return c_type::offDiag(noff); 3493 } 3494 3495 inline realpart_type realPart() 3496 { return c_type::realPart(); } 3497 3498 inline realpart_type imagPart() 3499 { return c_type::imagPart(); } 3500 3501 inline lowertri_type view() 3502 { return *this; } 3503 3504 inline lowertri_type viewAsUnitDiag() 3505 { return c_type::viewAsUnitDiag(); } 3506 3507 inline uppertri_type transpose() 3508 { return c_type::transpose(); } 3509 3510 inline lowertri_type conjugate() 3511 { return c_type::conjugate(); } 3512 3513 inline uppertri_type adjoint() 3514 { return c_type::adjoint(); } 3515 3516 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 3517 { 3518 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 3519 return c_type::cSubMatrix(i1-1,i2,j1-1,j2); 3520 } 3521 3522 inline const_rec_type subMatrix( 3523 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 3524 { 3525 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 3526 return c_type::cSubMatrix( 3527 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 3528 } 3529 3530 inline const_vec_type subVector( 3531 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t s) const 3532 { 3533 TMVAssert(hasSubVector(i,j,istep,jstep,s)); 3534 return c_type::cSubVector(i-1,j-1,istep,jstep,s); 3535 } 3536 3537 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 3538 { 3539 TMVAssert(hasSubTriMatrix(i1,i2,1)); 3540 return c_type::cSubTriMatrix(i1-1,i2); 3541 } 3542 3543 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 3544 { 3545 TMVAssert(hasSubTriMatrix(i1,i2,istep)); 3546 return c_type::cSubTriMatrix(i1-1,i2-1+istep,istep); 3547 } 3548 3549 inline const_lowertri_type offDiag(ptrdiff_t noff=1) const 3550 { 3551 TMVAssert(noff >= 1); 3552 TMVAssert(noff <= c_type::size()); 3553 return c_type::offDiag(noff); 3554 } 3555 3556 inline const_realpart_type realPart() const 3557 { return c_type::realPart(); } 3558 3559 inline const_realpart_type imagPart() const 3560 { return c_type::imagPart(); } 3561 3562 inline const_lowertri_type view() const 3563 { return *this; } 3564 3565 inline const_lowertri_type viewAsUnitDiag() const 3566 { return c_type::viewAsUnitDiag(); } 3567 3568 inline const_uppertri_type transpose() const 3569 { return c_type::transpose(); } 3570 3571 inline const_lowertri_type conjugate() const 3572 { return c_type::conjugate(); } 3573 3574 inline const_uppertri_type adjoint() const 3575 { return c_type::adjoint(); } 3576 3577 }; // FortranStyle LowerTriMatrixView 3578 3579 3580 template <typename T, int A> 3581 class UpperTriMatrix : public GenUpperTriMatrix<T> 3582 { 3583 public: 3584 3585 enum { D = A & UnitDiag }; 3586 enum { S = A & AllStorageType }; 3587 enum { I = A & FortranStyle }; 3588 3589 typedef TMV_RealType(T) RT; 3590 typedef TMV_ComplexType(T) CT; 3591 typedef GenUpperTriMatrix<T> base; 3592 typedef UpperTriMatrix<T,A> type; 3593 typedef ConstVectorView<T,I> const_vec_type; 3594 typedef ConstMatrixView<T,I> const_rec_type; 3595 typedef ConstUpperTriMatrixView<T,I> const_uppertri_type; 3596 typedef ConstLowerTriMatrixView<T,I> const_lowertri_type; 3597 typedef const_uppertri_type const_view_type; 3598 typedef const_lowertri_type const_transpose_type; 3599 typedef const_uppertri_type const_conjugate_type; 3600 typedef const_lowertri_type const_adjoint_type; 3601 typedef ConstUpperTriMatrixView<RT,I> const_realpart_type; 3602 typedef VectorView<T,I> vec_type; 3603 typedef MatrixView<T,I> rec_type; 3604 typedef UpperTriMatrixView<T,I> uppertri_type; 3605 typedef LowerTriMatrixView<T,I> lowertri_type; 3606 typedef uppertri_type view_type; 3607 typedef lowertri_type transpose_type; 3608 typedef uppertri_type conjugate_type; 3609 typedef lowertri_type adjoint_type; 3610 typedef UpperTriMatrixView<RT,I> realpart_type; 3611 typedef typename TriRefHelper2<T,D>::reference reference; 3612 typedef RMIt<type> rowmajor_iterator; 3613 typedef CRMIt<type> const_rowmajor_iterator; 3614 typedef CMIt<type> colmajor_iterator; 3615 typedef CCMIt<type> const_colmajor_iterator; 3616 3617 // 3618 // Constructors 3619 // 3620 3621 #define NEW_SIZE(s) \ 3622 itslen((s)*(s)), itsm(itslen), itss(s) \ 3623 TMV_DEFFIRSTLAST(itsm.get(),itsm.get()+itslen) 3624 3625 explicit inline UpperTriMatrix(ptrdiff_t _size=0) : NEW_SIZE(_size) 3626 { 3627 TMVAssert(Attrib<A>::trimatrixok); 3628 TMVAssert(_size >= 0); 3629 #ifdef TMV_EXTRA_DEBUG 3630 setAllTo(T(888)); 3631 #endif 3632 } 3633 3634 inline UpperTriMatrix(ptrdiff_t _size, const T& x) : NEW_SIZE(_size) 3635 { 3636 TMVAssert(Attrib<A>::trimatrixok); 3637 TMVAssert(_size >= 0); 3638 setAllTo(x); 3639 } 3640 3641 template <typename T2> 3642 inline UpperTriMatrix(const GenMatrix<T2>& rhs) : 3643 NEW_SIZE(rhs.rowsize()) 3644 { 3645 TMVAssert(Attrib<A>::trimatrixok); 3646 TMVAssert(isReal(T2()) || isComplex(T())); 3647 Copy(rhs.upperTri(dt()),view()); 3648 } 3649 3650 template <typename T2> 3651 inline UpperTriMatrix(const GenUpperTriMatrix<T2>& rhs) : 3652 NEW_SIZE(rhs.size()) 3653 { 3654 TMVAssert(Attrib<A>::trimatrixok); 3655 TMVAssert(isReal(T2()) || isComplex(T())); 3656 if (isunit() && !rhs.isunit()) { 3657 if (rhs.size() > 0) 3658 Copy(rhs.offDiag(),offDiag()); 3659 } else { 3660 Copy(rhs,view()); 3661 } 3662 } 3663 3664 inline UpperTriMatrix(const type& rhs) : 3665 itslen(rhs.itslen), itsm(itslen), itss(rhs.itss) 3666 TMV_DEFFIRSTLAST(itsm.get(),itsm.get()+itslen) 3667 { 3668 TMVAssert(Attrib<A>::trimatrixok); 3669 std::copy(rhs.cptr(),rhs.cptr()+itslen,itsm.get()); 3670 } 3671 3672 inline UpperTriMatrix(const GenMatrix<T>& rhs) : 3673 NEW_SIZE(rhs.rowsize()) 3674 { 3675 TMVAssert(Attrib<A>::trimatrixok); 3676 Copy(rhs.upperTri(dt()),view()); 3677 } 3678 3679 inline UpperTriMatrix(const GenUpperTriMatrix<RT>& rhs) : 3680 NEW_SIZE(rhs.size()) 3681 { 3682 TMVAssert(Attrib<A>::trimatrixok); 3683 if (isunit() && !rhs.isunit()) { 3684 if (rhs.size() > 0) offDiag() = rhs.offDiag(); 3685 } else 3686 rhs.assignToU(view()); 3687 } 3688 3689 inline UpperTriMatrix(const GenUpperTriMatrix<CT>& rhs) : 3690 NEW_SIZE(rhs.size()) 3691 { 3692 TMVAssert(Attrib<A>::trimatrixok); 3693 TMVAssert(isComplex(T())); 3694 if (isunit() && !rhs.isunit()) { 3695 if (rhs.size() > 0) offDiag() = rhs.offDiag(); 3696 } else 3697 rhs.assignToU(view()); 3698 } 3699 3700 inline UpperTriMatrix(const AssignableToUpperTriMatrix<RT>& m2) : 3701 NEW_SIZE(m2.size()) 3702 { 3703 TMVAssert(Attrib<A>::trimatrixok); 3704 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3705 m2.assignToU(view()); 3706 } 3707 3708 inline UpperTriMatrix(const AssignableToUpperTriMatrix<CT>& m2) : 3709 NEW_SIZE(m2.size()) 3710 { 3711 TMVAssert(Attrib<A>::trimatrixok); 3712 TMVAssert(isComplex(T())); 3713 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3714 m2.assignToU(view()); 3715 } 3716 3717 #undef NEW_SIZE 3718 3719 virtual inline ~UpperTriMatrix() 3720 { 3721 TMV_SETFIRSTLAST(0,0); 3722 #ifdef TMV_EXTRA_DEBUG 3723 setAllTo(T(999)); 3724 #endif 3725 } 3726 3727 3728 // 3729 // Op= 3730 // 3731 3732 inline type& operator=(const type& m2) 3733 { 3734 TMVAssert(size() == m2.size()); 3735 if (&m2 != this) 3736 std::copy(m2.cptr(),m2.cptr()+itslen,itsm.get()); 3737 return *this; 3738 } 3739 3740 inline type& operator=(const GenUpperTriMatrix<RT>& m2) 3741 { 3742 TMVAssert(size() == m2.size()); 3743 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3744 m2.assignToU(view()); 3745 return *this; 3746 } 3747 3748 inline type& operator=(const GenUpperTriMatrix<CT>& m2) 3749 { 3750 TMVAssert(size() == m2.size()); 3751 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3752 TMVAssert(isComplex(T())); 3753 m2.assignToU(view()); 3754 return *this; 3755 } 3756 3757 template <typename T2> 3758 inline type& operator=(const GenUpperTriMatrix<T2>& m2) 3759 { 3760 TMVAssert(size() == m2.size()); 3761 TMVAssert(isReal(T2()) || sComplex(T())); 3762 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3763 Copy(m2,view()); 3764 return *this; 3765 } 3766 3767 inline type& operator=(const T& x) 3768 { 3769 TMVAssert(!this->isunit() || x==T(1)); 3770 return setToIdentity(x); 3771 } 3772 3773 inline type& operator=(const AssignableToUpperTriMatrix<RT>& m2) 3774 { 3775 TMVAssert(size() == m2.size()); 3776 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3777 m2.assignToU(view()); 3778 return *this; 3779 } 3780 3781 inline type& operator=(const AssignableToUpperTriMatrix<CT>& m2) 3782 { 3783 TMVAssert(size() == m2.size()); 3784 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 3785 TMVAssert(isComplex(T())); 3786 m2.assignToU(view()); 3787 return *this; 3788 } 3789 3790 typedef ListAssigner<T,rowmajor_iterator> MyListAssigner; 3791 inline MyListAssigner operator<<(const T& x) 3792 { return MyListAssigner(rowmajor_begin(),size()*(size()+1)/2,x); } 3793 3794 // 3795 // Access 3796 // 3797 3798 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 3799 { 3800 if (I == int(CStyle)) { 3801 TMVAssert(i>=0 && i<size()); 3802 TMVAssert(j>=0 && j<size()); 3803 } else { 3804 TMVAssert(i>0 && i<=size()); --i; 3805 TMVAssert(j>0 && j<=size()); --j; 3806 } 3807 return cref(i,j); 3808 } 3809 3810 inline reference operator()(ptrdiff_t i, ptrdiff_t j) 3811 { 3812 if (I == int(CStyle)) { 3813 TMVAssert(i>=0 && i<size()); 3814 TMVAssert(j>=0 && j<size()); 3815 TMVAssert(i<=j); 3816 } else { 3817 TMVAssert(i>0 && i<= size()); --i; 3818 TMVAssert(j>0 && j<= size()); --j; 3819 TMVAssert(i<=j); 3820 } 3821 return ref(i,j); 3822 } 3823 3824 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3825 { 3826 if (I==int(FortranStyle)) { 3827 TMVAssert(i>0 && i<=size()); --i; 3828 TMVAssert(j1>0 && j1<=j2 && j2<=size()); --j1; 3829 } else { 3830 TMVAssert(i>=0 && i<size()); 3831 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 3832 } 3833 TMVAssert(j1==j2 || okij(i,j1)); 3834 return const_vec_type( 3835 itsm.get()+i*stepi()+j1*stepj(),j2-j1,stepj(),NonConj); 3836 } 3837 3838 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 3839 { 3840 if (I==int(FortranStyle)) { 3841 TMVAssert(j>0 && j<=size()); --j; 3842 TMVAssert(i1>0 && i1<=i2 && i2<=size()); --i1; 3843 } else { 3844 TMVAssert(j>=0 && j<size()); 3845 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 3846 } 3847 TMVAssert(i1==i2 || okij(i2-1,j)); 3848 return const_vec_type( 3849 itsm.get()+i1*stepi()+j*stepj(),i2-i1,stepi(),NonConj); 3850 } 3851 3852 inline const_vec_type diag() const 3853 { 3854 TMVAssert(!isunit()); 3855 return const_vec_type(itsm.get(),size(),stepi()+stepj(),NonConj); 3856 } 3857 3858 inline const_vec_type diag(ptrdiff_t i) const 3859 { 3860 TMVAssert(isunit() ? i>0 : i>=0); 3861 TMVAssert(i<=size()); 3862 return const_vec_type( 3863 itsm.get()+i*stepj(),size()-i,stepi()+stepj(),NonConj); 3864 } 3865 3866 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3867 { 3868 TMVAssert(isunit() ? i>0 : i>=0); 3869 TMVAssert(i<=size()); 3870 if (I == int(FortranStyle)) { 3871 TMVAssert(j1 > 0 && j1<=j2 && j2<=size()-i); --j1; 3872 } else { 3873 TMVAssert(j1>=0 && j1<=j2 && j2<=size()-i); 3874 } 3875 const ptrdiff_t ds = stepi()+stepj(); 3876 return const_vec_type( 3877 itsm.get()+i*stepj()+j1*ds,j2-j1,ds,NonConj); 3878 } 3879 3880 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3881 { 3882 if (I==int(FortranStyle)) { 3883 TMVAssert(i>0 && i<=size()); --i; 3884 TMVAssert(j1>0 && j1<=j2 && j2<=size()); --j1; 3885 } else { 3886 TMVAssert(i>=0 && i<size()); 3887 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 3888 } 3889 TMVAssert(j1==j2 || okij(i,j1)); 3890 return vec_type( 3891 itsm.get()+i*stepi()+j1*stepj(), 3892 j2-j1,stepj(),NonConj TMV_FIRSTLAST); 3893 } 3894 3895 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 3896 { 3897 if (I==int(FortranStyle)) { 3898 TMVAssert(j>0 && j<=size()); --j; 3899 TMVAssert(i1>0 && i1<=i2 && i2<=size()); --i1; 3900 } else { 3901 TMVAssert(j>=0 && j<size()); 3902 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 3903 } 3904 TMVAssert(i1==i2 || okij(i2-1,j)); 3905 return vec_type( 3906 itsm.get()+i1*stepi()+j*stepj(), 3907 i2-i1,stepi(),NonConj TMV_FIRSTLAST); 3908 } 3909 3910 inline vec_type diag() 3911 { 3912 TMVAssert(!isunit()); 3913 return vec_type( 3914 itsm.get(),size(),stepi()+stepj(),NonConj TMV_FIRSTLAST); 3915 } 3916 3917 inline vec_type diag(ptrdiff_t i) 3918 { 3919 TMVAssert(isunit() ? i>0 : i>=0); 3920 TMVAssert(i<=size()); 3921 return vec_type( 3922 itsm.get()+i*stepj(),size()-i,stepi()+stepj(),NonConj 3923 TMV_FIRSTLAST); 3924 } 3925 3926 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3927 { 3928 TMVAssert(isunit() ? i>0 : i>=0); 3929 TMVAssert(i<=size()); 3930 if (I == int(FortranStyle)) { 3931 TMVAssert(j1 > 0 && j1<=j2 && j2<=size()-i); --j1; 3932 } else { 3933 TMVAssert(j1>=0 && j1<=j2 && j2<=size()-i); 3934 } 3935 const ptrdiff_t ds = stepi()+stepj(); 3936 return vec_type( 3937 itsm.get()+i*stepj()+j1*ds,j2-j1,ds,NonConj TMV_FIRSTLAST); 3938 } 3939 3940 // 3941 // Modifying Functions 3942 // 3943 3944 inline type& setZero() 3945 { fill_n(itsm.get(),itslen,T(0)); return *this; } 3946 3947 inline type& setAllTo(const T& x) 3948 { VectorViewOf(itsm.get(),itslen).setAllTo(x); return *this; } 3949 3950 inline type& addToAll(const T& x) 3951 { 3952 TMVAssert(!isunit()); 3953 VectorViewOf(itsm.get(),itslen).addToAll(x); 3954 return *this; 3955 } 3956 3957 inline type& clip(RT thresh) 3958 { VectorViewOf(itsm.get(),itslen).clip(thresh); return *this; } 3959 3960 inline type& conjugateSelf() 3961 { VectorViewOf(itsm.get(),itslen).conjugateSelf(); return *this; } 3962 3963 inline type& invertSelf() 3964 { view().invertSelf(); return *this; } 3965 3966 inline type& setToIdentity(const T& x=T(1)) 3967 { 3968 TMVAssert(!isunit() || x==T(1)); 3969 setZero(); if (!isunit()) diag().setAllTo(x); 3970 return *this; 3971 } 3972 3973 // 3974 // subMatrix 3975 // 3976 3977 inline const_rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 3978 { 3979 return const_rec_type( 3980 itsm.get()+i1*stepi()+j1*stepj(), 3981 i2-i1, j2-j1,stepi(),stepj(),NonConj); 3982 } 3983 3984 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 3985 { 3986 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 3987 if (I==int(FortranStyle)) { --i1; --j1; } 3988 return cSubMatrix(i1,i2,j1,j2); 3989 } 3990 3991 inline const_rec_type cSubMatrix( 3992 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 3993 { 3994 return const_rec_type( 3995 itsm.get()+i1*stepi()+j1*stepj(), 3996 (i2-i1)/istep, (j2-j1)/jstep, istep*stepi(), jstep*stepj(), 3997 NonConj); 3998 } 3999 4000 inline const_rec_type subMatrix( 4001 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 4002 { 4003 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 4004 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 4005 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 4006 } 4007 4008 inline const_vec_type cSubVector( 4009 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 4010 { 4011 return const_vec_type( 4012 itsm.get()+i*stepi()+j*stepj(),size, 4013 istep*stepi()+jstep*stepj(),NonConj); 4014 } 4015 4016 inline const_vec_type subVector( 4017 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 4018 { 4019 TMVAssert(size >= 0); 4020 TMVAssert(view().hasSubVector(i,j,istep,jstep,size)); 4021 if (I==int(FortranStyle)) { --i; --j; } 4022 return cSubVector(i,j,istep,jstep,size); 4023 } 4024 4025 inline const_uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 4026 { 4027 return const_uppertri_type( 4028 itsm.get()+i1*(stepi()+stepj()), 4029 i2-i1,stepi(),stepj(),dt(),NonConj); 4030 } 4031 4032 inline const_uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 4033 { 4034 TMVAssert(view().hasSubTriMatrix(i1,i2,1)); 4035 if (I==int(FortranStyle)) { --i1; } 4036 return cSubTriMatrix(i1,i2); 4037 } 4038 4039 inline const_uppertri_type cSubTriMatrix( 4040 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 4041 { 4042 return const_uppertri_type( 4043 itsm.get()+i1*(stepi()+stepj()), 4044 (i2-i1)/istep, istep*stepi(),istep*stepj(), dt(), NonConj); 4045 } 4046 4047 inline const_uppertri_type subTriMatrix( 4048 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 4049 { 4050 TMVAssert(view().hasSubTriMatrix(i1,i2,istep)); 4051 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 4052 return cSubTriMatrix(i1,i2,istep); 4053 } 4054 4055 inline const_uppertri_type offDiag(ptrdiff_t noff=1) const 4056 { 4057 TMVAssert(noff >= 1); 4058 TMVAssert(noff <= size()); 4059 return const_uppertri_type( 4060 itsm.get()+noff*stepj(), 4061 size()-noff,stepi(),stepj(),NonUnitDiag,NonConj); 4062 } 4063 4064 inline const_realpart_type realPart() const 4065 { 4066 return const_realpart_type( 4067 reinterpret_cast<RT*>(itsm.get()), size(), 4068 isReal(T()) ? stepi() : 2*stepi(), 4069 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj); 4070 } 4071 4072 inline const_realpart_type imagPart() const 4073 { 4074 TMVAssert(isComplex(T())); 4075 TMVAssert(!isunit()); 4076 return const_realpart_type( 4077 reinterpret_cast<RT*>(itsm.get())+1, size(), 4078 2*stepi(), 2*stepj(), NonUnitDiag, NonConj); 4079 } 4080 4081 inline rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 4082 { 4083 return rec_type( 4084 itsm.get()+i1*stepi()+j1*stepj(), 4085 i2-i1, j2-j1, stepi(),stepj(),NonConj TMV_FIRSTLAST); 4086 } 4087 4088 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 4089 { 4090 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 4091 if (I==int(FortranStyle)) { --i1; --j1; } 4092 return cSubMatrix(i1,i2,j1,j2); 4093 } 4094 4095 inline rec_type cSubMatrix( 4096 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 4097 { 4098 return rec_type( 4099 itsm.get()+i1*stepi()+j1*stepj(), 4100 (i2-i1)/istep, (j2-j1)/jstep, istep*stepi(), jstep*stepj(), 4101 NonConj TMV_FIRSTLAST); 4102 } 4103 4104 inline rec_type subMatrix( 4105 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 4106 { 4107 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 4108 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 4109 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 4110 } 4111 4112 inline vec_type cSubVector(ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 4113 { 4114 return vec_type( 4115 itsm.get()+i*stepi()+j*stepj(),size, 4116 istep*stepi()+jstep*stepj(),NonConj TMV_FIRSTLAST); 4117 } 4118 4119 inline vec_type subVector(ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 4120 { 4121 TMVAssert(size >= 0); 4122 TMVAssert(view().hasSubVector(i,j,istep,jstep,size)); 4123 if (I==int(FortranStyle)) { --i; --j; } 4124 return cSubVector(i,j,istep,jstep,size); 4125 } 4126 4127 inline uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 4128 { 4129 return uppertri_type( 4130 itsm.get()+i1*(stepi()+stepj()), 4131 i2-i1,stepi(),stepj(),dt(),NonConj TMV_FIRSTLAST); 4132 } 4133 4134 inline uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 4135 { 4136 TMVAssert(view().hasSubTriMatrix(i1,i2,1)); 4137 if (I==int(FortranStyle)) { --i1; } 4138 return cSubTriMatrix(i1,i2); 4139 } 4140 4141 inline uppertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 4142 { 4143 return uppertri_type( 4144 itsm.get()+i1*(stepi()+stepj()), 4145 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(), NonConj 4146 TMV_FIRSTLAST); 4147 } 4148 4149 inline uppertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 4150 { 4151 TMVAssert(view().hasSubTriMatrix(i1,i2,istep)); 4152 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 4153 return cSubTriMatrix(i1,i2,istep); 4154 } 4155 4156 inline uppertri_type offDiag(ptrdiff_t noff=1) 4157 { 4158 TMVAssert(noff >= 1); 4159 TMVAssert(noff <= size()); 4160 return uppertri_type( 4161 itsm.get()+noff*stepj(),size()-noff, 4162 stepi(),stepj(),NonUnitDiag,NonConj TMV_FIRSTLAST); 4163 } 4164 4165 inline realpart_type realPart() 4166 { 4167 return realpart_type( 4168 reinterpret_cast<RT*>(itsm.get()), size(), 4169 isReal(T()) ? stepi() : 2*stepi(), 4170 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj 4171 #ifdef TMVFLDEBUG 4172 ,reinterpret_cast<const RT*>(_first) 4173 ,reinterpret_cast<const RT*>(_last) 4174 #endif 4175 ); 4176 } 4177 4178 inline realpart_type imagPart() 4179 { 4180 TMVAssert(isComplex(T())); 4181 TMVAssert(!isunit()); 4182 return realpart_type( 4183 reinterpret_cast<RT*>(itsm.get())+1, size(), 4184 2*stepi(), 2*stepj(), NonUnitDiag, NonConj 4185 #ifdef TMVFLDEBUG 4186 ,reinterpret_cast<const RT*>(_first)+1 4187 ,reinterpret_cast<const RT*>(_last)+1 4188 #endif 4189 ); 4190 } 4191 4192 inline const_uppertri_type view() const 4193 { 4194 return const_uppertri_type( 4195 itsm.get(),size(),stepi(),stepj(),dt(),NonConj); 4196 } 4197 4198 inline const_uppertri_type viewAsUnitDiag() const 4199 { 4200 return const_uppertri_type( 4201 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj); 4202 } 4203 4204 inline const_lowertri_type transpose() const 4205 { 4206 return const_lowertri_type( 4207 itsm.get(),size(),stepj(),stepi(),dt(),NonConj); 4208 } 4209 4210 inline const_uppertri_type conjugate() const 4211 { 4212 return const_uppertri_type( 4213 itsm.get(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,NonConj)); 4214 } 4215 4216 inline const_lowertri_type adjoint() const 4217 { 4218 return const_lowertri_type( 4219 itsm.get(),size(),stepj(),stepi(),dt(), 4220 TMV_ConjOf(T,NonConj)); 4221 } 4222 4223 inline uppertri_type view() 4224 { 4225 return uppertri_type( 4226 itsm.get(),size(),stepi(),stepj(),dt(),NonConj TMV_FIRSTLAST); 4227 } 4228 4229 inline uppertri_type viewAsUnitDiag() 4230 { 4231 return uppertri_type( 4232 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj 4233 TMV_FIRSTLAST); 4234 } 4235 4236 inline lowertri_type transpose() 4237 { 4238 return lowertri_type( 4239 itsm.get(),size(),stepj(),stepi(),dt(),NonConj TMV_FIRSTLAST); 4240 } 4241 4242 inline uppertri_type conjugate() 4243 { 4244 return uppertri_type( 4245 itsm.get(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,NonConj) 4246 TMV_FIRSTLAST); 4247 } 4248 4249 inline lowertri_type adjoint() 4250 { 4251 return lowertri_type( 4252 itsm.get(),size(),stepj(),stepi(),dt(), 4253 TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 4254 } 4255 4256 // 4257 // I/O 4258 // 4259 4260 void read(const TMV_Reader& reader); 4261 4262 inline ptrdiff_t size() const { return itss; } 4263 inline const T* cptr() const { return itsm.get(); } 4264 inline T* ptr() { return itsm.get(); } 4265 inline ptrdiff_t stepi() const { return S==int(RowMajor) ? itss : 1; } 4266 inline ptrdiff_t stepj() const { return S==int(RowMajor) ? 1 : itss; } 4267 inline DiagType dt() const { return static_cast<DiagType>(D); } 4268 inline ConjType ct() const { return NonConj; } 4269 inline bool isrm() const { return S==int(RowMajor); } 4270 inline bool iscm() const { return S==int(ColMajor); } 4271 inline bool isunit() const { return D == int(UnitDiag); } 4272 inline bool isconj() const { return false; } 4273 4274 inline reference ref(ptrdiff_t i, ptrdiff_t j) 4275 { 4276 return TriRefHelper2<T,D>::makeRef( 4277 i==j, 4278 itsm.get()[S==int(RowMajor) ? i*itss + j : j*itss + i]); 4279 } 4280 4281 inline T cref(ptrdiff_t i, ptrdiff_t j) const 4282 { 4283 return ( 4284 (isunit() && i==j) ? T(1) : 4285 (i>j) ? T(0) : 4286 itsm.get()[S==int(RowMajor) ? i*itss + j : j*itss + i]); 4287 } 4288 4289 inline void resize(ptrdiff_t s) 4290 { 4291 TMVAssert(s >= 0); 4292 itslen = s*s; 4293 itsm.resize(itslen); 4294 itss = s; 4295 #ifdef TMVFLDEBUG 4296 _first = itsm.get(); 4297 _last = _first+itslen; 4298 #endif 4299 #ifdef TMV_EXTRA_DEBUG 4300 setAllTo(T(888)); 4301 #endif 4302 } 4303 4304 inline rowmajor_iterator rowmajor_begin() 4305 { return rowmajor_iterator(this,0,0); } 4306 inline rowmajor_iterator rowmajor_end() 4307 { return rowmajor_iterator(this,size(),size()); } 4308 4309 inline const_rowmajor_iterator rowmajor_begin() const 4310 { return const_rowmajor_iterator(this,0,0); } 4311 inline const_rowmajor_iterator rowmajor_end() const 4312 { return const_rowmajor_iterator(this,size(),size()); } 4313 4314 inline colmajor_iterator colmajor_begin() 4315 { return colmajor_iterator(this,0,0); } 4316 inline colmajor_iterator colmajor_end() 4317 { return colmajor_iterator(this,0,size()); } 4318 4319 inline const_colmajor_iterator colmajor_begin() const 4320 { return const_colmajor_iterator(this,0,0); } 4321 inline const_colmajor_iterator colmajor_end() const 4322 { return const_colmajor_iterator(this,0,size()); } 4323 4324 protected : 4325 4326 ptrdiff_t itslen; 4327 AlignedArray<T> itsm; 4328 ptrdiff_t itss; 4329 4330 #ifdef TMVFLDEBUG 4331 public : 4332 const T* _first; 4333 const T* _last; 4334 protected : 4335 #endif 4336 4337 inline bool okij(ptrdiff_t i, ptrdiff_t j) const 4338 { 4339 TMVAssert(i>=0 && i < size()); 4340 TMVAssert(j>=0 && j < size()); 4341 return isunit() ? i-j<0 : i-j<=0; 4342 } 4343 4344 friend void Swap( 4345 UpperTriMatrix<T,A>& m1, UpperTriMatrix<T,A>& m2) 4346 { 4347 TMVAssert(m1.size() == m2.size()); 4348 m1.itsm.swapWith(m2.itsm); 4349 #ifdef TMVFLDEBUG 4350 TMV_SWAP(m1._first,m2._first); 4351 TMV_SWAP(m1._last,m2._last); 4352 #endif 4353 } 4354 4355 }; // UpperTriMatrix 4356 4357 template <typename T, int A> 4358 class LowerTriMatrix : public GenLowerTriMatrix<T> 4359 { 4360 public: 4361 4362 enum { D = A & UnitDiag }; 4363 enum { S = A & AllStorageType }; 4364 enum { I = A & FortranStyle }; 4365 4366 typedef TMV_RealType(T) RT; 4367 typedef TMV_ComplexType(T) CT; 4368 typedef GenLowerTriMatrix<T> base; 4369 typedef LowerTriMatrix<T,A> type; 4370 typedef ConstVectorView<T,I> const_vec_type; 4371 typedef ConstMatrixView<T,I> const_rec_type; 4372 typedef ConstUpperTriMatrixView<T,I> const_uppertri_type; 4373 typedef ConstLowerTriMatrixView<T,I> const_lowertri_type; 4374 typedef const_lowertri_type const_view_type; 4375 typedef const_uppertri_type const_transpose_type; 4376 typedef const_lowertri_type const_conjugate_type; 4377 typedef const_uppertri_type const_adjoint_type; 4378 typedef ConstLowerTriMatrixView<RT,I> const_realpart_type; 4379 typedef VectorView<T,I> vec_type; 4380 typedef MatrixView<T,I> rec_type; 4381 typedef UpperTriMatrixView<T,I> uppertri_type; 4382 typedef LowerTriMatrixView<T,I> lowertri_type; 4383 typedef lowertri_type view_type; 4384 typedef uppertri_type transpose_type; 4385 typedef lowertri_type conjugate_type; 4386 typedef uppertri_type adjoint_type; 4387 typedef LowerTriMatrixView<RT,I> realpart_type; 4388 typedef typename TriRefHelper2<T,D>::reference reference; 4389 typedef RMIt<type> rowmajor_iterator; 4390 typedef CRMIt<type> const_rowmajor_iterator; 4391 typedef CMIt<type> colmajor_iterator; 4392 typedef CCMIt<type> const_colmajor_iterator; 4393 4394 // 4395 // Constructors 4396 // 4397 4398 #define NEW_SIZE(s) \ 4399 itslen((s)*(s)), itsm(itslen), itss(s) \ 4400 TMV_DEFFIRSTLAST(itsm.get(),itsm.get()+itslen) 4401 4402 explicit inline LowerTriMatrix(ptrdiff_t _size=0) : NEW_SIZE(_size) 4403 { 4404 TMVAssert(Attrib<A>::trimatrixok); 4405 TMVAssert(_size >= 0); 4406 #ifdef TMV_EXTRA_DEBUG 4407 setAllTo(T(888)); 4408 #endif 4409 } 4410 4411 inline LowerTriMatrix(ptrdiff_t _size, const T& x) : NEW_SIZE(_size) 4412 { 4413 TMVAssert(Attrib<A>::trimatrixok); 4414 TMVAssert(_size >= 0); 4415 setAllTo(x); 4416 } 4417 4418 template <typename T2> 4419 inline LowerTriMatrix(const GenMatrix<T2>& rhs) : 4420 NEW_SIZE(rhs.colsize()) 4421 { 4422 TMVAssert(Attrib<A>::trimatrixok); 4423 TMVAssert(isReal(T2()) || isComplex(T())); 4424 Copy(rhs.lowerTri(dt()).transpose(),transpose()); 4425 } 4426 4427 template <typename T2> 4428 inline LowerTriMatrix(const GenLowerTriMatrix<T2>& rhs) : 4429 NEW_SIZE(rhs.size()) 4430 { 4431 TMVAssert(Attrib<A>::trimatrixok); 4432 TMVAssert(isReal(T2()) || isComplex(T())); 4433 if (isunit() && !rhs.isunit()) { 4434 if (rhs.size() > 0) 4435 Copy(rhs.offDiag().transpose(),offDiag().transpose()); 4436 } else { 4437 Copy(rhs.transpose(),transpose()); 4438 } 4439 } 4440 4441 inline LowerTriMatrix(const type& rhs) : 4442 itslen(rhs.itslen), itsm(itslen), itss(rhs.itss) 4443 TMV_DEFFIRSTLAST(itsm.get(),itsm.get()+itslen) 4444 { 4445 TMVAssert(Attrib<A>::trimatrixok); 4446 std::copy(rhs.cptr(),rhs.cptr()+itslen,itsm.get()); 4447 } 4448 4449 inline LowerTriMatrix(const GenMatrix<T>& rhs) : 4450 NEW_SIZE(rhs.rowsize()) 4451 { 4452 TMVAssert(Attrib<A>::trimatrixok); 4453 Copy(rhs.lowerTri(dt()).transpose(),transpose()); 4454 } 4455 4456 inline LowerTriMatrix(const GenLowerTriMatrix<RT>& rhs) : 4457 NEW_SIZE(rhs.size()) 4458 { 4459 TMVAssert(Attrib<A>::trimatrixok); 4460 if (isunit() && !rhs.isunit()) { 4461 if (rhs.size() > 0) offDiag() = rhs.offDiag(); 4462 } else 4463 rhs.assignToL(view()); 4464 } 4465 4466 inline LowerTriMatrix(const GenLowerTriMatrix<CT>& rhs) : 4467 NEW_SIZE(rhs.size()) 4468 { 4469 TMVAssert(Attrib<A>::trimatrixok); 4470 TMVAssert(isComplex(T())); 4471 if (isunit() && !rhs.isunit()) { 4472 if (rhs.size() > 0) offDiag() = rhs.offDiag(); 4473 } else 4474 rhs.assignToL(view()); 4475 } 4476 4477 inline LowerTriMatrix(const AssignableToLowerTriMatrix<RT>& m2) : 4478 NEW_SIZE(m2.size()) 4479 { 4480 TMVAssert(Attrib<A>::trimatrixok); 4481 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4482 m2.assignToL(view()); 4483 } 4484 4485 inline LowerTriMatrix(const AssignableToLowerTriMatrix<CT>& m2) : 4486 NEW_SIZE(m2.size()) 4487 { 4488 TMVAssert(Attrib<A>::trimatrixok); 4489 TMVAssert(isComplex(T())); 4490 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4491 m2.assignToL(view()); 4492 } 4493 4494 #undef NEW_SIZE 4495 4496 virtual inline ~LowerTriMatrix() 4497 { 4498 TMV_SETFIRSTLAST(0,0); 4499 #ifdef TMV_EXTRA_DEBUG 4500 setAllTo(T(999)); 4501 #endif 4502 } 4503 4504 4505 // 4506 // Op= 4507 // 4508 4509 inline type& operator=(const type& m2) 4510 { 4511 TMVAssert(size() == m2.size()); 4512 if (&m2 != this) 4513 std::copy(m2.cptr(),m2.cptr()+itslen,itsm.get()); 4514 return *this; 4515 } 4516 4517 inline type& operator=(const GenLowerTriMatrix<RT>& m2) 4518 { 4519 TMVAssert(size() == m2.size()); 4520 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4521 m2.assignToL(view()); 4522 return *this; 4523 } 4524 4525 inline type& operator=(const GenLowerTriMatrix<CT>& m2) 4526 { 4527 TMVAssert(size() == m2.size()); 4528 TMVAssert(isComplex(T())); 4529 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4530 m2.assignToL(view()); 4531 return *this; 4532 } 4533 4534 template <typename T2> 4535 inline type& operator=(const GenLowerTriMatrix<T2>& m2) 4536 { 4537 TMVAssert(size() == m2.size()); 4538 TMVAssert(isReal(T2()) || isComplex(T())); 4539 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4540 Copy(m2,view()); 4541 return *this; 4542 } 4543 4544 inline type& operator=(const T& x) 4545 { 4546 TMVAssert(!this->isunit() || x==T(1)); 4547 return setToIdentity(x); 4548 } 4549 4550 inline type& operator=(const AssignableToLowerTriMatrix<RT>& m2) 4551 { 4552 TMVAssert(size() == m2.size()); 4553 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4554 m2.assignToL(view()); 4555 return *this; 4556 } 4557 4558 inline type& operator=(const AssignableToLowerTriMatrix<CT>& m2) 4559 { 4560 TMVAssert(size() == m2.size()); 4561 TMVAssert(!(m2.dt()==NonUnitDiag && dt()==UnitDiag)); 4562 TMVAssert(isComplex(T())); 4563 m2.assignToL(view()); 4564 return *this; 4565 } 4566 4567 typedef ListAssigner<T,rowmajor_iterator> MyListAssigner; 4568 inline MyListAssigner operator<<(const T& x) 4569 { return MyListAssigner(rowmajor_begin(),size()*(size()+1)/2,x); } 4570 4571 // 4572 // Access 4573 // 4574 4575 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 4576 { 4577 if (I == int(CStyle)) { 4578 TMVAssert(i>=0 && i<size()); 4579 TMVAssert(j>=0 && j<size()); 4580 } else { 4581 TMVAssert(i>0 && i<=size()); --i; 4582 TMVAssert(j>0 && j<=size()); --j; 4583 } 4584 return cref(i,j); 4585 } 4586 4587 inline reference operator()(ptrdiff_t i, ptrdiff_t j) 4588 { 4589 if (I == int(CStyle)) { 4590 TMVAssert(i>=0 && i<size()); 4591 TMVAssert(j>=0 && j<size()); 4592 TMVAssert(i>=j); 4593 } else { 4594 TMVAssert(i>0 && i<= size()); --i; 4595 TMVAssert(j>0 && j<= size()); --j; 4596 TMVAssert(i>=j); 4597 } 4598 return ref(i,j); 4599 } 4600 4601 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 4602 { 4603 if (I==int(FortranStyle)) { 4604 TMVAssert(i>0 && i<=size()); --i; 4605 TMVAssert(j1>0 && j1<=j2 && j2<=size()); --j1; 4606 } else { 4607 TMVAssert(i>=0 && i<size()); 4608 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 4609 } 4610 TMVAssert(j1==j2 || okij(i,j2-1)); 4611 return const_vec_type( 4612 itsm.get()+i*stepi()+j1*stepj(),j2-j1,stepj(),NonConj); 4613 } 4614 4615 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 4616 { 4617 if (I==int(FortranStyle)) { 4618 TMVAssert(j>0 && j<=size()); --j; 4619 TMVAssert(i1>0 && i1<=i2 && i2<=size()); --i1; 4620 } else { 4621 TMVAssert(j>=0 && j<size()); 4622 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 4623 } 4624 TMVAssert(i1==i2 || okij(i1,j)); 4625 return const_vec_type( 4626 itsm.get()+i1*stepi()+j*stepj(),i2-i1,stepi(),NonConj); 4627 } 4628 4629 inline const_vec_type diag() const 4630 { 4631 TMVAssert(!isunit()); 4632 return const_vec_type(itsm.get(),size(),stepi()+stepj(),NonConj); 4633 } 4634 4635 inline const_vec_type diag(ptrdiff_t i) const 4636 { 4637 TMVAssert(i>=-size()); 4638 TMVAssert(isunit() ? i<0 : i<=0); 4639 return const_vec_type( 4640 itsm.get()-i*stepi(),size()+i,stepi()+stepj(),NonConj); 4641 } 4642 4643 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 4644 { 4645 TMVAssert(i>=-size()); 4646 TMVAssert(isunit() ? i<0 : i<=0); 4647 if (I == int(FortranStyle)) { 4648 TMVAssert(j1>0 && j1<=j2 && j2<=size()+i); --j1; 4649 } else { 4650 TMVAssert(j1>=0 && j1<=j2 && j2<=size()+i); 4651 } 4652 const ptrdiff_t ds = stepi()+stepj(); 4653 return const_vec_type( 4654 itsm.get()-i*stepi()+j1*ds,j2-j1,ds,NonConj); 4655 } 4656 4657 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 4658 { 4659 if (I==int(FortranStyle)) { 4660 TMVAssert(i>0 && i<=size()); --i; 4661 TMVAssert(j1>0 && j1<=j2 && j2<=size()); --j1; 4662 } else { 4663 TMVAssert(i>=0 && i<size()); 4664 TMVAssert(j1>=0 && j1<=j2 && j2<=size()); 4665 } 4666 TMVAssert(j1==j2 || okij(i,j2-1)); 4667 return vec_type( 4668 itsm.get()+i*stepi()+j1*stepj(), 4669 j2-j1,stepj(),NonConj TMV_FIRSTLAST); 4670 } 4671 4672 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 4673 { 4674 if (I==int(FortranStyle)) { 4675 TMVAssert(j>0 && j<=size()); --j; 4676 TMVAssert(i1>0 && i1<=i2 && i2<=size()); --i1; 4677 } else { 4678 TMVAssert(j>=0 && j<size()); 4679 TMVAssert(i1>=0 && i1<=i2 && i2<=size()); 4680 } 4681 TMVAssert(i1==i2 || okij(i1,j)); 4682 return vec_type( 4683 itsm.get()+i1*stepi()+j*stepj(),i2-i1,stepi(),NonConj 4684 TMV_FIRSTLAST); 4685 } 4686 4687 inline vec_type diag() 4688 { 4689 TMVAssert(!isunit()); 4690 return vec_type( 4691 itsm.get(),size(),stepi()+stepj(),NonConj TMV_FIRSTLAST); 4692 } 4693 4694 inline vec_type diag(ptrdiff_t i) 4695 { 4696 TMVAssert(i>=-size()); 4697 TMVAssert(isunit() ? i<0 : i<=0); 4698 return vec_type( 4699 itsm.get()-i*stepi(),size()+i,stepi()+stepj(),NonConj 4700 TMV_FIRSTLAST); 4701 } 4702 4703 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 4704 { 4705 TMVAssert(i>=-size()); 4706 TMVAssert(isunit() ? i<0 : i<=0); 4707 if (I == int(FortranStyle)) { 4708 TMVAssert(j1>0 && j1<=j2 && j2<=size()+i); --j1; 4709 } else { 4710 TMVAssert(j1>=0 && j1<=j2 && j2<=size()+i); 4711 } 4712 const ptrdiff_t ds = stepi()+stepj(); 4713 return vec_type( 4714 itsm.get()-i*stepi()+j1*ds,j2-j1,ds,NonConj TMV_FIRSTLAST); 4715 } 4716 4717 // 4718 // Modifying Functions 4719 // 4720 4721 inline type& setZero() 4722 { fill_n(itsm.get(),itslen,T(0)); return *this; } 4723 4724 inline type& setAllTo(const T& x) 4725 { VectorViewOf(itsm.get(),itslen).setAllTo(x); return *this; } 4726 4727 inline type& addToAll(const T& x) 4728 { 4729 TMVAssert(!isunit()); 4730 VectorViewOf(itsm.get(),itslen).addToAll(x); 4731 return *this; 4732 } 4733 4734 inline type& clip(RT thresh) 4735 { VectorViewOf(itsm.get(),itslen).clip(thresh); return *this; } 4736 4737 inline type& conjugateSelf() 4738 { VectorViewOf(itsm.get(),itslen).conjugateSelf(); return *this; } 4739 4740 inline type& invertSelf() 4741 { view().invertSelf(); return *this; } 4742 4743 inline type& setToIdentity(const T& x=T(1)) 4744 { 4745 TMVAssert(!isunit() || x == T(1)); 4746 setZero(); if (!isunit()) diag().setAllTo(x); 4747 return *this; 4748 } 4749 4750 // 4751 // subMatrix 4752 // 4753 4754 inline const_rec_type cSubMatrix( 4755 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 4756 { 4757 return const_rec_type( 4758 itsm.get()+i1*stepi()+j1*stepj(), 4759 i2-i1, j2-j1,stepi(),stepj(),NonConj); 4760 } 4761 4762 inline const_rec_type subMatrix( 4763 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 4764 { 4765 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 4766 if (I==int(FortranStyle)) { --i1; --j1; } 4767 return cSubMatrix(i1,i2,j1,j2); 4768 } 4769 4770 inline const_rec_type cSubMatrix( 4771 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 4772 { 4773 return const_rec_type( 4774 itsm.get()+i1*stepi()+j1*stepj(), 4775 (i2-i1)/istep, (j2-j1)/jstep, istep*stepi(), jstep*stepj(), 4776 NonConj); 4777 } 4778 4779 inline const_rec_type subMatrix( 4780 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 4781 { 4782 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 4783 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 4784 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 4785 } 4786 4787 inline const_vec_type cSubVector( 4788 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 4789 { 4790 TMVAssert(size >= 0); 4791 return const_vec_type( 4792 itsm.get()+i*stepi()+j*stepj(),size, 4793 istep*stepi()+jstep*stepj(),NonConj); 4794 } 4795 4796 inline const_vec_type subVector( 4797 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 4798 { 4799 TMVAssert(size >= 0); 4800 TMVAssert(view().hasSubVector(i,j,istep,jstep,size)); 4801 if (I==int(FortranStyle)) { --i; --j; } 4802 return cSubVector(i,j,istep,jstep,size); 4803 } 4804 4805 inline const_lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 4806 { 4807 return const_lowertri_type( 4808 itsm.get()+i1*(stepi()+stepj()), 4809 i2-i1,stepi(),stepj(),dt(),NonConj); 4810 } 4811 4812 inline const_lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) const 4813 { 4814 TMVAssert(view().hasSubTriMatrix(i1,i2,1)); 4815 if (I==int(FortranStyle)) { --i1; } 4816 return cSubTriMatrix(i1,i2); 4817 } 4818 4819 inline const_lowertri_type cSubTriMatrix( 4820 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 4821 { 4822 return const_lowertri_type( 4823 itsm.get()+i1*(stepi()+stepj()), 4824 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(), NonConj); 4825 } 4826 4827 inline const_lowertri_type subTriMatrix( 4828 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 4829 { 4830 TMVAssert(view().hasSubTriMatrix(i1,i2,istep)); 4831 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 4832 return cSubTriMatrix(i1,i2,istep); 4833 } 4834 4835 inline const_lowertri_type offDiag(ptrdiff_t noff=1) const 4836 { 4837 TMVAssert(noff >= 1); 4838 TMVAssert(noff <= size()); 4839 return const_lowertri_type( 4840 itsm.get()+noff*stepi(),size()-noff, 4841 stepi(),stepj(),NonUnitDiag,NonConj); 4842 } 4843 4844 inline const_realpart_type realPart() const 4845 { 4846 return const_realpart_type( 4847 reinterpret_cast<RT*>(itsm.get()), size(), 4848 isReal(T()) ? stepi() : 2*stepi(), 4849 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj 4850 ); 4851 } 4852 4853 inline const_realpart_type imagPart() const 4854 { 4855 TMVAssert(isComplex(T())); 4856 TMVAssert(!isunit()); 4857 return const_realpart_type( 4858 reinterpret_cast<RT*>(itsm.get())+1, size(), 4859 2*stepi(), 2*stepj(), NonUnitDiag, NonConj 4860 ); 4861 } 4862 4863 inline rec_type cSubMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 4864 { 4865 return rec_type( 4866 itsm.get()+i1*stepi()+j1*stepj(), 4867 i2-i1, j2-j1, stepi(),stepj(),NonConj TMV_FIRSTLAST ); 4868 } 4869 4870 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 4871 { 4872 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 4873 if (I==int(FortranStyle)) { --i1; --j1; } 4874 return cSubMatrix(i1,i2,j1,j2); 4875 } 4876 4877 inline rec_type cSubMatrix( 4878 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 4879 { 4880 return rec_type( 4881 itsm.get()+i1*stepi()+j1*stepj(), 4882 (i2-i1)/istep, (j2-j1)/jstep, istep*stepi(), jstep*stepj(), 4883 NonConj TMV_FIRSTLAST ); 4884 } 4885 4886 inline rec_type subMatrix( 4887 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 4888 { 4889 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 4890 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 4891 return cSubMatrix(i1,i2,j1,j2,istep,jstep); 4892 } 4893 4894 inline vec_type cSubVector( 4895 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 4896 { 4897 TMVAssert(size >= 0); 4898 return vec_type( 4899 itsm.get()+i*stepi()+j*stepj(),size, 4900 istep*stepi()+jstep*stepj(),NonConj TMV_FIRSTLAST ); 4901 } 4902 4903 inline vec_type subVector( 4904 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) 4905 { 4906 TMVAssert(size >= 0); 4907 TMVAssert(view().hasSubVector(i,j,istep,jstep,size)); 4908 if (I==int(FortranStyle)) { --i; --j; } 4909 return cSubVector(i,j,istep,jstep,size); 4910 } 4911 4912 inline lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 4913 { 4914 return lowertri_type( 4915 itsm.get()+i1*(stepi()+stepj()), 4916 i2-i1,stepi(),stepj(),dt(),NonConj TMV_FIRSTLAST); 4917 } 4918 4919 inline lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2) 4920 { 4921 TMVAssert(view().hasSubTriMatrix(i1,i2,1)); 4922 if (I==int(FortranStyle)) { --i1; } 4923 return cSubTriMatrix(i1,i2); 4924 } 4925 4926 inline lowertri_type cSubTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 4927 { 4928 return lowertri_type( 4929 itsm.get()+i1*(stepi()+stepj()), 4930 (i2-i1)/istep,istep*stepi(),istep*stepj(),dt(), NonConj 4931 TMV_FIRSTLAST); 4932 } 4933 4934 inline lowertri_type subTriMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 4935 { 4936 TMVAssert(view().hasSubTriMatrix(i1,i2,istep)); 4937 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 4938 return cSubTriMatrix(i1,i2,istep); 4939 } 4940 4941 inline lowertri_type offDiag(ptrdiff_t noff=1) 4942 { 4943 TMVAssert(noff >= 1); 4944 TMVAssert(noff <= size()); 4945 return lowertri_type( 4946 itsm.get()+noff*stepi(),size()-noff, 4947 stepi(),stepj(),NonUnitDiag,NonConj TMV_FIRSTLAST); 4948 } 4949 4950 inline realpart_type realPart() 4951 { 4952 return realpart_type( 4953 reinterpret_cast<RT*>(itsm.get()), size(), 4954 isReal(T()) ? stepi() : 2*stepi(), 4955 isReal(T()) ? stepj() : 2*stepj(), dt(), NonConj 4956 #ifdef TMVFLDEBUG 4957 ,reinterpret_cast<const RT*>(_first) 4958 ,reinterpret_cast<const RT*>(_last) 4959 #endif 4960 ); 4961 } 4962 4963 inline realpart_type imagPart() 4964 { 4965 TMVAssert(isComplex(T())); 4966 TMVAssert(!isunit()); 4967 return realpart_type( 4968 reinterpret_cast<RT*>(itsm.get())+1, 4969 size(), 2*stepi(), 2*stepj(), NonUnitDiag, NonConj 4970 #ifdef TMVFLDEBUG 4971 ,reinterpret_cast<const RT*>(_first)+1 4972 ,reinterpret_cast<const RT*>(_last)+1 4973 #endif 4974 ); 4975 } 4976 4977 inline const_lowertri_type view() const 4978 { 4979 return const_lowertri_type( 4980 itsm.get(),size(),stepi(),stepj(),dt(),NonConj); 4981 } 4982 4983 inline const_lowertri_type viewAsUnitDiag() const 4984 { 4985 return const_lowertri_type( 4986 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj); 4987 } 4988 4989 inline const_uppertri_type transpose() const 4990 { 4991 return const_uppertri_type( 4992 itsm.get(),size(),stepj(),stepi(),dt(),NonConj); 4993 } 4994 4995 inline const_lowertri_type conjugate() const 4996 { 4997 return const_lowertri_type( 4998 itsm.get(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,NonConj)); 4999 } 5000 5001 inline const_uppertri_type adjoint() const 5002 { 5003 return const_uppertri_type( 5004 itsm.get(),size(),stepj(),stepi(),dt(), 5005 TMV_ConjOf(T,NonConj)); 5006 } 5007 5008 inline lowertri_type view() 5009 { 5010 return lowertri_type( 5011 itsm.get(),size(),stepi(),stepj(),dt(),NonConj TMV_FIRSTLAST); 5012 } 5013 5014 inline lowertri_type viewAsUnitDiag() 5015 { 5016 return lowertri_type( 5017 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj 5018 TMV_FIRSTLAST); 5019 } 5020 5021 inline uppertri_type transpose() 5022 { 5023 return uppertri_type( 5024 itsm.get(),size(),stepj(),stepi(),dt(), 5025 NonConj TMV_FIRSTLAST); 5026 } 5027 5028 inline lowertri_type conjugate() 5029 { 5030 return lowertri_type( 5031 itsm.get(),size(),stepi(),stepj(),dt(),TMV_ConjOf(T,NonConj) 5032 TMV_FIRSTLAST); 5033 } 5034 5035 inline uppertri_type adjoint() 5036 { 5037 return uppertri_type( 5038 itsm.get(),size(),stepj(),stepi(),dt(), 5039 TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 5040 } 5041 5042 // 5043 // I/O 5044 // 5045 5046 void read(const TMV_Reader& reader); 5047 5048 inline ptrdiff_t size() const { return itss; } 5049 inline const T* cptr() const { return itsm.get(); } 5050 inline T* ptr() { return itsm.get(); } 5051 inline ptrdiff_t stepi() const { return S==int(RowMajor) ? itss : 1; } 5052 inline ptrdiff_t stepj() const { return S==int(RowMajor) ? 1 : itss; } 5053 inline DiagType dt() const { return static_cast<DiagType>(D); } 5054 inline ConjType ct() const { return NonConj; } 5055 inline bool isrm() const { return S==int(RowMajor); } 5056 inline bool iscm() const { return S==int(ColMajor); } 5057 inline bool isunit() const { return D == int(UnitDiag); } 5058 inline bool isconj() const { return false; } 5059 5060 inline reference ref(ptrdiff_t i, ptrdiff_t j) 5061 { 5062 return TriRefHelper2<T,D>::makeRef( 5063 i==j, 5064 itsm.get()[S==int(RowMajor) ? i*itss + j : j*itss + i]); 5065 } 5066 5067 inline T cref(ptrdiff_t i, ptrdiff_t j) const 5068 { 5069 return ( 5070 (isunit() && i==j) ? T(1) : 5071 (i<j) ? T(0) : 5072 itsm.get()[S==int(RowMajor) ? i*itss + j : j*itss + i]); 5073 } 5074 5075 inline void resize(ptrdiff_t s) 5076 { 5077 TMVAssert(s >= 0); 5078 itslen = s*s; 5079 itsm.resize(itslen); 5080 itss = s; 5081 #ifdef TMVFLDEBUG 5082 _first = itsm.get(); 5083 _last = _first+itslen; 5084 #endif 5085 #ifdef TMV_EXTRA_DEBUG 5086 setAllTo(T(888)); 5087 #endif 5088 } 5089 5090 inline rowmajor_iterator rowmajor_begin() 5091 { return rowmajor_iterator(this,0,0); } 5092 inline rowmajor_iterator rowmajor_end() 5093 { return rowmajor_iterator(this,size(),0); } 5094 5095 inline const_rowmajor_iterator rowmajor_begin() const 5096 { return const_rowmajor_iterator(this,0,0); } 5097 inline const_rowmajor_iterator rowmajor_end() const 5098 { return const_rowmajor_iterator(this,size(),0); } 5099 5100 inline colmajor_iterator colmajor_begin() 5101 { return colmajor_iterator(this,0,0); } 5102 inline colmajor_iterator colmajor_end() 5103 { return colmajor_iterator(this,size(),size()); } 5104 5105 inline const_colmajor_iterator colmajor_begin() const 5106 { return const_colmajor_iterator(this,0,0); } 5107 inline const_colmajor_iterator colmajor_end() const 5108 { return const_colmajor_iterator(this,size(),size()); } 5109 5110 protected : 5111 5112 ptrdiff_t itslen; 5113 AlignedArray<T> itsm; 5114 ptrdiff_t itss; 5115 5116 #ifdef TMVFLDEBUG 5117 public: 5118 const T* _first; 5119 const T* _last; 5120 protected : 5121 #endif 5122 5123 inline bool okij(ptrdiff_t i, ptrdiff_t j) const 5124 { 5125 TMVAssert(i>=0 && i < size()); 5126 TMVAssert(j>=0 && j < size()); 5127 return isunit() ? i-j>0 : i-j>=0; 5128 } 5129 5130 friend void Swap( 5131 LowerTriMatrix<T,A>& m1, LowerTriMatrix<T,A>& m2) 5132 { 5133 TMVAssert(m1.size() == m2.size()); 5134 m1.itsm.swapWith(m2.itsm); 5135 #ifdef TMVFLDEBUG 5136 TMV_SWAP(m1._first,m2._first); 5137 TMV_SWAP(m1._last,m2._last); 5138 #endif 5139 } 5140 5141 }; // LowerTriMatrix 5142 5143 //------------------------------------------------------------------------- 5144 5145 // 5146 // Special Creators: 5147 // UpperTriMatrixViewOf(T* m, n, S, D) 5148 // UnitUpperTriMatrixViewOf(T* m, n, S) 5149 // UpperTriMatrixViewOf(T* m, n, Si, Sj, D) 5150 // UnitUpperTriMatrixViewOf(T* m, n, Si, Sj) 5151 // 5152 5153 template <typename T> 5154 inline ConstUpperTriMatrixView<T> UpperTriMatrixViewOf( 5155 const T* m, ptrdiff_t size, StorageType stor, DiagType dt=NonUnitDiag) 5156 { 5157 TMVAssert(size >= 0); 5158 TMVAssert(stor == RowMajor || stor == ColMajor); 5159 if (stor == RowMajor) 5160 return ConstUpperTriMatrixView<T>(m,size,size,1,dt,NonConj); 5161 else 5162 return ConstUpperTriMatrixView<T>(m,size,1,size,dt,NonConj); 5163 } 5164 5165 template <typename T> 5166 inline UpperTriMatrixView<T> UpperTriMatrixViewOf( 5167 T* m, ptrdiff_t size, StorageType stor, DiagType dt=NonUnitDiag) 5168 { 5169 TMVAssert(size >= 0); 5170 TMVAssert(stor == RowMajor || stor == ColMajor); 5171 if (stor == RowMajor) 5172 return UpperTriMatrixView<T>( 5173 m,size,size,1,dt,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5174 else 5175 return UpperTriMatrixView<T>( 5176 m,size,1,size,dt,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5177 } 5178 5179 template <typename T> 5180 inline ConstLowerTriMatrixView<T> LowerTriMatrixViewOf( 5181 const T* m, ptrdiff_t size, StorageType stor, DiagType dt=NonUnitDiag) 5182 { 5183 TMVAssert(size >= 0); 5184 TMVAssert(stor == RowMajor || stor == ColMajor); 5185 if (stor == RowMajor) 5186 return ConstLowerTriMatrixView<T>(m,size,size,1,dt,NonConj); 5187 else 5188 return ConstLowerTriMatrixView<T>(m,size,1,size,dt,NonConj); 5189 } 5190 5191 template <typename T> 5192 inline LowerTriMatrixView<T> LowerTriMatrixViewOf( 5193 T* m, ptrdiff_t size, StorageType stor, DiagType dt=NonUnitDiag) 5194 { 5195 TMVAssert(size >= 0); 5196 TMVAssert(stor == RowMajor || stor == ColMajor); 5197 if (stor == RowMajor) 5198 return LowerTriMatrixView<T>( 5199 m,size,size,1,dt,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5200 else 5201 return LowerTriMatrixView<T>( 5202 m,size,1,size,dt,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5203 } 5204 5205 template <typename T> 5206 inline ConstUpperTriMatrixView<T> UnitUpperTriMatrixViewOf( 5207 const T* m, ptrdiff_t size, StorageType stor) 5208 { 5209 TMVAssert(size >= 0); 5210 TMVAssert(stor == RowMajor || stor == ColMajor); 5211 if (stor == RowMajor) 5212 return ConstUpperTriMatrixView<T>(m,size,size,1,UnitDiag,NonConj); 5213 else 5214 return ConstUpperTriMatrixView<T>(m,size,1,size,UnitDiag,NonConj); 5215 } 5216 5217 template <typename T> 5218 inline UpperTriMatrixView<T> UnitUpperTriMatrixViewOf( 5219 T* m, ptrdiff_t size, StorageType stor) 5220 { 5221 TMVAssert(size >= 0); 5222 TMVAssert(stor == RowMajor || stor == ColMajor); 5223 if (stor == RowMajor) 5224 return UpperTriMatrixView<T>( 5225 m,size,size,1,UnitDiag,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5226 else 5227 return UpperTriMatrixView<T>( 5228 m,size,1,size,UnitDiag,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5229 } 5230 5231 template <typename T> 5232 inline ConstLowerTriMatrixView<T> UnitLowerTriMatrixViewOf( 5233 const T* m, ptrdiff_t size, StorageType stor) 5234 { 5235 TMVAssert(size >= 0); 5236 TMVAssert(stor == RowMajor || stor == ColMajor); 5237 if (stor == RowMajor) 5238 return ConstLowerTriMatrixView<T>(m,size,size,1,UnitDiag,NonConj); 5239 else 5240 return ConstLowerTriMatrixView<T>(m,size,1,size,UnitDiag,NonConj); 5241 } 5242 5243 template <typename T> 5244 inline LowerTriMatrixView<T> UnitLowerTriMatrixViewOf( 5245 T* m, ptrdiff_t size, StorageType stor) 5246 { 5247 TMVAssert(size >= 0); 5248 TMVAssert(stor == RowMajor || stor == ColMajor); 5249 if (stor == RowMajor) 5250 return LowerTriMatrixView<T>( 5251 m,size,size,1,UnitDiag,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5252 else 5253 return LowerTriMatrixView<T>( 5254 m,size,1,size,UnitDiag,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5255 } 5256 5257 5258 template <typename T> 5259 inline ConstUpperTriMatrixView<T> UpperTriMatrixViewOf( 5260 const T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj, DiagType dt=NonUnitDiag) 5261 { 5262 TMVAssert(size >= 0); 5263 return ConstUpperTriMatrixView<T>(m,size,stepi,stepj,dt,NonConj); 5264 } 5265 5266 template <typename T> 5267 inline UpperTriMatrixView<T> UpperTriMatrixViewOf( 5268 T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj, DiagType dt=NonUnitDiag) 5269 { 5270 TMVAssert(size >= 0); 5271 return UpperTriMatrixView<T>( 5272 m,size,stepi,stepj,dt,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5273 } 5274 5275 template <typename T> 5276 inline ConstLowerTriMatrixView<T> LowerTriMatrixViewOf( 5277 const T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj, DiagType dt=NonUnitDiag) 5278 { 5279 TMVAssert(size >= 0); 5280 return ConstLowerTriMatrixView<T>(m,size,stepi,stepj,dt,NonConj); 5281 } 5282 5283 template <typename T> 5284 inline LowerTriMatrixView<T> LowerTriMatrixViewOf( 5285 T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj, DiagType dt=NonUnitDiag) 5286 { 5287 TMVAssert(size >= 0); 5288 return LowerTriMatrixView<T>( 5289 m,size,stepi,stepj,dt,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5290 } 5291 5292 template <typename T> 5293 inline ConstUpperTriMatrixView<T> UnitUpperTriMatrixViewOf( 5294 const T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj) 5295 { 5296 TMVAssert(size >= 0); 5297 return ConstUpperTriMatrixView<T>(m,size,stepi,stepj,UnitDiag,NonConj); 5298 } 5299 5300 template <typename T> 5301 inline UpperTriMatrixView<T> UnitUpperTriMatrixViewOf( 5302 T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj) 5303 { 5304 TMVAssert(size >= 0); 5305 return UpperTriMatrixView<T>( 5306 m,size,size,1,UnitDiag,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5307 } 5308 5309 template <typename T> 5310 inline ConstLowerTriMatrixView<T> UnitLowerTriMatrixViewOf( 5311 const T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj) 5312 { 5313 TMVAssert(size >= 0); 5314 return ConstLowerTriMatrixView<T>(m,size,size,1,UnitDiag,NonConj); 5315 } 5316 5317 template <typename T> 5318 inline LowerTriMatrixView<T> UnitLowerTriMatrixViewOf( 5319 T* m, ptrdiff_t size, ptrdiff_t stepi, ptrdiff_t stepj) 5320 { 5321 TMVAssert(size >= 0); 5322 return LowerTriMatrixView<T>( 5323 m,size,size,1,UnitDiag,NonConj TMV_FIRSTLAST1(m,m+size*size)); 5324 } 5325 5326 // 5327 // Copy 5328 // 5329 5330 template <typename T1, typename T2> 5331 inline void nonUnitDiagCopy( 5332 const GenUpperTriMatrix<T1>& m1, UpperTriMatrixView<T2> m2) 5333 { 5334 TMVAssert(isReal(T1()) || isComplex(T2())); 5335 TMVAssert(m1.size() == m2.size()); 5336 TMVAssert(m1.dt() == NonUnitDiag); 5337 TMVAssert(m2.dt() == NonUnitDiag); 5338 const ptrdiff_t N = m1.size(); 5339 5340 if (!m1.isSameAs(m2) && m1.size() > 0) { 5341 if (m1.iscm() && m2.iscm()) 5342 for(ptrdiff_t j=0;j<N;++j) m2.col(j,0,j+1) = m1.col(j,0,j+1); 5343 else 5344 for(ptrdiff_t i=0;i<N;++i) m2.row(i,i,N) = m1.row(i,i,N); 5345 } 5346 } 5347 5348 template <typename T1, typename T2> 5349 inline void Copy( 5350 const GenUpperTriMatrix<T1>& m1, UpperTriMatrixView<T2> m2) 5351 { 5352 TMVAssert(isReal(T1()) || isComplex(T2())); 5353 TMVAssert(m1.size() == m2.size()); 5354 TMVAssert(m1.isunit() || !m2.isunit()); 5355 5356 if (m1.isunit()) { 5357 if (m1.size() > 0) 5358 nonUnitDiagCopy(m1.offDiag(),m2.offDiag()); 5359 if (!m2.isunit()) 5360 m2.diag().setAllTo(T2(1)); 5361 } else { 5362 nonUnitDiagCopy(m1,m2); 5363 } 5364 } 5365 5366 5367 // 5368 // Swap Matrices 5369 // 5370 5371 template <typename T> 5372 void Swap( 5373 UpperTriMatrixView<T> m1, UpperTriMatrixView<T> m2); 5374 5375 template <typename T, int A> 5376 inline void Swap( 5377 UpperTriMatrixView<T> m1, UpperTriMatrix<T,A>& m2) 5378 { 5379 TMVAssert(m1.isunit() == m2.isunit()); 5380 Swap(m1,m2.view()); 5381 } 5382 5383 template <typename T, int A> 5384 inline void Swap( 5385 UpperTriMatrix<T,A>& m1, UpperTriMatrixView<T> m2) 5386 { 5387 TMVAssert(m1.isunit() == m2.isunit()); 5388 Swap(m1.view(),m2); 5389 } 5390 5391 template <typename T, int A1, int A2> 5392 inline void Swap( 5393 UpperTriMatrix<T,A1>& m1, UpperTriMatrix<T,A2>& m2) 5394 { 5395 TMVAssert(m1.isunit() == m2.isunit()); 5396 Swap(m1.view(),m2.view()); 5397 } 5398 5399 template <typename T> 5400 inline void Swap( 5401 LowerTriMatrixView<T> m1, LowerTriMatrixView<T> m2) 5402 { Swap(m1.transpose(),m2.transpose()); } 5403 5404 template <typename T, int A> 5405 inline void Swap( 5406 LowerTriMatrixView<T> m1, LowerTriMatrix<T,A>& m2) 5407 { Swap(m1.transpose(),m2.transpose()); } 5408 5409 template <typename T, int A> 5410 inline void Swap( 5411 LowerTriMatrix<T,A>& m1, LowerTriMatrixView<T> m2) 5412 { Swap(m1.transpose(),m2.transpose()); } 5413 5414 template <typename T, int A1, int A2> 5415 inline void Swap( 5416 LowerTriMatrix<T,A1>& m1, LowerTriMatrix<T,A2>& m2) 5417 { Swap(m1.transpose(),m2.transpose()); } 5418 5419 5420 // 5421 // Views: 5422 // 5423 5424 template <typename T> 5425 inline ConstLowerTriMatrixView<T> Transpose(const GenUpperTriMatrix<T>& m) 5426 { return m.transpose(); } 5427 template <typename T> 5428 inline ConstUpperTriMatrixView<T> Transpose(const GenLowerTriMatrix<T>& m) 5429 { return m.transpose(); } 5430 5431 template <typename T, int A> 5432 inline ConstLowerTriMatrixView<T,A> Transpose( 5433 const ConstUpperTriMatrixView<T,A>& m) 5434 { return m.transpose(); } 5435 template <typename T, int A> 5436 inline ConstUpperTriMatrixView<T,A> Transpose( 5437 const ConstLowerTriMatrixView<T,A>& m) 5438 { return m.transpose(); } 5439 5440 template <typename T, int A> 5441 inline ConstLowerTriMatrixView<T,A&FortranStyle> Transpose( 5442 const UpperTriMatrix<T,A>& m) 5443 { return m.transpose(); } 5444 template <typename T, int A> 5445 inline ConstUpperTriMatrixView<T,A&FortranStyle> Transpose( 5446 const LowerTriMatrix<T,A>& m) 5447 { return m.transpose(); } 5448 5449 template <typename T, int A> 5450 inline LowerTriMatrixView<T,A> Transpose(UpperTriMatrixView<T,A> m) 5451 { return m.transpose(); } 5452 template <typename T, int A> 5453 inline UpperTriMatrixView<T,A> Transpose(LowerTriMatrixView<T,A> m) 5454 { return m.transpose(); } 5455 5456 template <typename T, int A> 5457 inline LowerTriMatrixView<T,A&FortranStyle> Transpose( 5458 UpperTriMatrix<T,A>& m) 5459 { return m.transpose(); } 5460 template <typename T, int A> 5461 inline UpperTriMatrixView<T,A&FortranStyle> Transpose( 5462 LowerTriMatrix<T,A>& m) 5463 { return m.transpose(); } 5464 5465 template <typename T> 5466 inline ConstUpperTriMatrixView<T> Conjugate(const GenUpperTriMatrix<T>& m) 5467 { return m.conjugate(); } 5468 template <typename T> 5469 inline ConstLowerTriMatrixView<T> Conjugate(const GenLowerTriMatrix<T>& m) 5470 { return m.conjugate(); } 5471 5472 template <typename T, int A> 5473 inline ConstUpperTriMatrixView<T,A> Conjugate( 5474 const ConstUpperTriMatrixView<T,A>& m) 5475 { return m.conjugate(); } 5476 template <typename T, int A> 5477 inline ConstLowerTriMatrixView<T,A> Conjugate( 5478 const ConstLowerTriMatrixView<T,A>& m) 5479 { return m.conjugate(); } 5480 5481 template <typename T, int A> 5482 inline ConstUpperTriMatrixView<T,A&FortranStyle> Conjugate( 5483 const UpperTriMatrix<T,A>& m) 5484 { return m.conjugate(); } 5485 template <typename T, int A> 5486 inline ConstLowerTriMatrixView<T,A&FortranStyle> Conjugate( 5487 const LowerTriMatrix<T,A>& m) 5488 { return m.conjugate(); } 5489 5490 template <typename T, int A> 5491 inline UpperTriMatrixView<T,A> Conjugate(UpperTriMatrixView<T,A> m) 5492 { return m.conjugate(); } 5493 template <typename T, int A> 5494 inline LowerTriMatrixView<T,A> Conjugate(LowerTriMatrixView<T,A> m) 5495 { return m.conjugate(); } 5496 5497 template <typename T, int A> 5498 inline UpperTriMatrixView<T,A&FortranStyle> Conjugate( 5499 UpperTriMatrix<T,A>& m) 5500 { return m.conjugate(); } 5501 template <typename T, int A> 5502 inline LowerTriMatrixView<T,A&FortranStyle> Conjugate( 5503 LowerTriMatrix<T,A>& m) 5504 { return m.conjugate(); } 5505 5506 template <typename T> 5507 inline ConstLowerTriMatrixView<T> Adjoint(const GenUpperTriMatrix<T>& m) 5508 { return m.adjoint(); } 5509 template <typename T> 5510 inline ConstUpperTriMatrixView<T> Adjoint(const GenLowerTriMatrix<T>& m) 5511 { return m.adjoint(); } 5512 5513 template <typename T, int A> 5514 inline ConstLowerTriMatrixView<T,A> Adjoint( 5515 const ConstUpperTriMatrixView<T,A>& m) 5516 { return m.adjoint(); } 5517 template <typename T, int A> 5518 inline ConstUpperTriMatrixView<T,A> Adjoint( 5519 const ConstLowerTriMatrixView<T,A>& m) 5520 { return m.adjoint(); } 5521 5522 template <typename T, int A> 5523 inline ConstLowerTriMatrixView<T,A&FortranStyle> Adjoint( 5524 const UpperTriMatrix<T,A>& m) 5525 { return m.adjoint(); } 5526 template <typename T, int A> 5527 inline ConstUpperTriMatrixView<T,A&FortranStyle> Adjoint( 5528 const LowerTriMatrix<T,A>& m) 5529 { return m.adjoint(); } 5530 5531 template <typename T, int A> 5532 inline LowerTriMatrixView<T,A> Adjoint(UpperTriMatrixView<T,A> m) 5533 { return m.adjoint(); } 5534 template <typename T, int A> 5535 inline UpperTriMatrixView<T,A> Adjoint(LowerTriMatrixView<T,A> m) 5536 { return m.adjoint(); } 5537 5538 template <typename T, int A> 5539 inline LowerTriMatrixView<T,A&FortranStyle> Adjoint(UpperTriMatrix<T,A>& m) 5540 { return m.adjoint(); } 5541 template <typename T, int A> 5542 inline UpperTriMatrixView<T,A&FortranStyle> Adjoint(LowerTriMatrix<T,A>& m) 5543 { return m.adjoint(); } 5544 5545 template <typename T> 5546 inline QuotXU<T,T> Inverse(const GenUpperTriMatrix<T>& m) 5547 { return m.inverse(); } 5548 template <typename T> 5549 inline QuotXL<T,T> Inverse(const GenLowerTriMatrix<T>& m) 5550 { return m.inverse(); } 5551 5552 5553 // 5554 // TriMatrix ==, != TriMatrix 5555 // 5556 5557 template <typename T1, typename T2> 5558 bool operator==( 5559 const GenUpperTriMatrix<T1>& m1, const GenUpperTriMatrix<T2>& m2); 5560 5561 template <typename T1, typename T2> 5562 inline bool operator==( 5563 const GenLowerTriMatrix<T1>& m1, const GenLowerTriMatrix<T2>& m2) 5564 { return m1.transpose() == m2.transpose(); } 5565 5566 template <typename T1, typename T2> 5567 inline bool operator!=( 5568 const GenUpperTriMatrix<T1>& m1, const GenUpperTriMatrix<T2>& m2) 5569 { return !(m1 == m2); } 5570 5571 template <typename T1, typename T2> 5572 inline bool operator!=( 5573 const GenLowerTriMatrix<T1>& m1, const GenLowerTriMatrix<T2>& m2) 5574 { return !(m1 == m2); } 5575 5576 template <typename T1, typename T2> 5577 inline bool operator==( 5578 const GenUpperTriMatrix<T1>& m1, const GenMatrix<T2>& m2) 5579 { 5580 return 5581 m1 == m2.upperTri() && 5582 m2.lowerTri().offDiag().maxAbs2Element() == T2(0); 5583 } 5584 5585 template <typename T1, typename T2> 5586 inline bool operator==( 5587 const GenLowerTriMatrix<T1>& m1, const GenMatrix<T2>& m2) 5588 { 5589 return 5590 m1 == m2.lowerTri() && 5591 m2.upperTri().offDiag().maxAbs2Element() == T2(0); 5592 } 5593 5594 template <typename T1, typename T2> 5595 inline bool operator==( 5596 const GenMatrix<T1>& m1, const GenUpperTriMatrix<T2>& m2) 5597 { return m2 == m1; } 5598 5599 template <typename T1, typename T2> 5600 inline bool operator==( 5601 const GenMatrix<T1>& m1, const GenLowerTriMatrix<T2>& m2) 5602 { return m2 == m1; } 5603 5604 template <typename T1, typename T2> 5605 inline bool operator!=( 5606 const GenUpperTriMatrix<T1>& m1, const GenMatrix<T2>& m2) 5607 { return !(m1 == m2); } 5608 5609 template <typename T1, typename T2> 5610 inline bool operator!=( 5611 const GenMatrix<T1>& m1, const GenUpperTriMatrix<T2>& m2) 5612 { return !(m1 == m2); } 5613 5614 template <typename T1, typename T2> 5615 inline bool operator!=( 5616 const GenLowerTriMatrix<T1>& m1, const GenMatrix<T2>& m2) 5617 { return !(m1 == m2); } 5618 5619 template <typename T1, typename T2> 5620 inline bool operator!=( 5621 const GenMatrix<T1>& m1, const GenLowerTriMatrix<T2>& m2) 5622 { return !(m1 == m2); } 5623 5624 5625 // 5626 // I/O 5627 // 5628 5629 template <typename T> 5630 inline std::istream& operator>>( 5631 std::istream& is, UpperTriMatrixView<T> m) 5632 { return is >> IOStyle() >> m; } 5633 5634 template <typename T, int A> 5635 inline std::istream& operator>>(std::istream& is, UpperTriMatrix<T,A>& m) 5636 { return is >> IOStyle() >> m; } 5637 5638 template <typename T> 5639 inline std::istream& operator>>( 5640 std::istream& is, LowerTriMatrixView<T> m) 5641 { return is >> IOStyle() >> m; } 5642 5643 template <typename T, int A> 5644 inline std::istream& operator>>(std::istream& is, LowerTriMatrix<T,A>& m) 5645 { return is >> IOStyle() >> m; } 5646 5647 template <typename T> 5648 inline std::istream& operator>>( 5649 const TMV_Reader& reader, UpperTriMatrixView<T> m) 5650 { m.read(reader); return reader.getis(); } 5651 5652 template <typename T, int A> 5653 inline std::istream& operator>>( 5654 const TMV_Reader& reader, UpperTriMatrix<T,A>& m) 5655 { m.read(reader); return reader.getis(); } 5656 5657 template <typename T> 5658 inline std::istream& operator>>( 5659 const TMV_Reader& reader, LowerTriMatrixView<T> m) 5660 { m.read(reader); return reader.getis(); } 5661 5662 template <typename T, int A> 5663 inline std::istream& operator>>( 5664 const TMV_Reader& reader, LowerTriMatrix<T,A>& m) 5665 { m.read(reader); return reader.getis(); } 5666 5667 5668 template <typename T, bool C> 5669 inline std::string TMV_Text(const TriRef<T,C>& ref) 5670 { 5671 return std::string("TriRef<") + TMV_Text(T()) + "," + 5672 TMV_Text(C ? Conj : NonConj) + ">"; 5673 } 5674 5675 } // namespace tmv 5676 5677 #endif 5678