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 SymMatrix and HermMatrix classes. 26 // 27 // SymMatrix is used for symmetric matrices, and 28 // HermMatrix is used for Hermitian matrices. 29 // For real matrices, these are the same thing: 30 // A = A.transpose(). 31 // But for complex, they are different: 32 // A_sym = A_sym.transpose() 33 // A_herm = A_herm.adjoint() 34 // 35 // For these notes, I will always write SymMatrix, but (except where 36 // otherwise indicated) everything applies the same for Sym and Herm. 37 // Also, the Views keep track of sym/herm difference with a parameter, 38 // so it is always a GenSymMatrix, ConstSymMatrixView, or 39 // SymMatrixView - never Herm in any of these. 40 // 41 // Caveat: Complex Hermitian matrices are such that A = At, which 42 // implies that their diagonal elements are real. Many routines 43 // involving HermMatrixes assume the reality of the diagonal. 44 // However, it is possible to assign a non-real value to a diagonal 45 // element. If the user does this, the results are undefined. 46 // 47 // As usual, the first template parameter is the type of the data, 48 // and the optional second template parameter specifies the known 49 // attributes. The valid attributs for a SymMatrix are: 50 // - ColMajor or RoMajor 51 // - CStyle or FortranStyle 52 // - Lower or Upper 53 // The defaults are (ColMajor,CStyle,Lower) if you do not specify 54 // otherwise. 55 // 56 // The storage and index style follow the same meaning as for regular 57 // Matrices. 58 // Lower or Upper refers to which triangular half stores the actual data. 59 // 60 // Constructors: 61 // 62 // SymMatrix<T,A>(int n) 63 // Makes a Symmetric Matrix with column size = row size = n 64 // with _uninitialized_ values. 65 // 66 // SymMatrix<T,A>(int n, T x) 67 // Makes a Symmetric Matrix with values all initialized to x 68 // For Hermitian matrixces, x must be real. 69 // 70 // SymMatrix<T,A>(const Matrix<T>& m) 71 // Makes a SymMatrix which copies the corresponding elements of m. 72 // 73 // ConstSymMatrixView<T> SymMatrixViewOf(const Matrix<T>& m, uplo) 74 // ConstSymMatrixView<T> HermMatrixViewOf(const Matrix<T>& m, uplo) 75 // Makes a constant SymMatrix view of the corresponding part of m. 76 // While this view cannot be modified, changing the original m 77 // will cause corresponding changes in this view of m. 78 // 79 // SymMatrixView<T> SymMatrixViewOf(Matrix<T>& m, uplo) 80 // SymMatrixView<T> HermMatrixViewOf(Matrix<T>& m, uplo) 81 // Makes a modifiable SymMatrix view of the corresponding part of m. 82 // Modifying this matrix will change the corresponding elements in 83 // the original Matrix. 84 // 85 // ConstSymMatrixView<T> SymMatrixViewOf(const T* m, size, uplo, stor) 86 // ConstSymMatrixView<T> HermMatrixViewOf(const T* m, size, uplo, stor) 87 // SymMatrixView<T> SymMatrixViewOf(T* m, size, uplo, stor) 88 // SymMatrixView<T> HermMatrixViewOf(T* m, size, uplo, stor) 89 // View the actual memory pointed to by m as a SymMatrix/HermMatrix 90 // with the given size and storage. 91 // 92 // Access Functions 93 // 94 // int colsize() const 95 // int rowsize() const 96 // int size() const 97 // Return the dimensions of the SymMatrix 98 // 99 // T& operator()(int i, int j) 100 // T operator()(int i, int j) const 101 // Return the (i,j) element of the SymMatrix 102 // 103 // Vector& row(int i, int j1, int j2) 104 // Return a portion of the ith row 105 // This range must be a valid range for the requested row 106 // and be entirely in the upper or lower triangle. 107 // 108 // Vector& col(int j, int i1, int i2) 109 // Return a portion of the jth column 110 // This range must be a valid range for the requested column 111 // and be entirely in the upper or lower triangle. 112 // 113 // Vector& diag() 114 // Return the main diagonal 115 // 116 // Vector& diag(int i, int j1, int j2) 117 // Vector& diag(int i) 118 // Return the super- or sub-diagonal i 119 // If i > 0 return the super diagonal starting at m_0i 120 // If i < 0 return the sub diagonal starting at m_|i|0 121 // If j1,j2 are given, it returns the diagonal Vector 122 // either from m_j1,i+j1 to m_j2,i+j2 (for i>0) 123 // or from m_|i|+j1,j1 to m_|i|+j2,j2 (for i<0) 124 // 125 // Modifying Functions: 126 // 127 // setZero() 128 // setAllTo(T x) 129 // addToAll(T x) 130 // For HermMatrix, x must be real in both setAllTo and addToAll. 131 // clip(RT thresh) 132 // conjugateSelf() 133 // transposeSelf() 134 // setToIdentity(x = 1) 135 // swapRowsCols(i1,i2) 136 // Equivalent to swapping rows i1,i2 then swapping cols i1,i2. 137 // permuteRowsCols(const int* p) 138 // reversePermuteRowsCols(const int* p) 139 // Perform a series of row/col swaps. 140 // void Swap(SymMatrix& m1, SymMatrix& m2) 141 // The SymMatrices must be the same size and Hermitianity. 142 // 143 // Views of a SymMatrix: 144 // 145 // subMatrix(int i1, int i2, int j1, int j2, int istep=1, int jstep=1) 146 // This member function will return a submatrix using rows i1 to i2 147 // and columns j1 to j2 which refers 148 // to the same physical elements as the original. 149 // The submatrix must be completely contained within the upper 150 // or lower triangle. 151 // 152 // subVector(int i, int j, int istep, int jstep, int size) 153 // Returns a VectorView which starts at position (i,j) in the 154 // matrix, moves in the directions (istep,jstep) and has a length 155 // of size. The subvector must be completely with the upper or 156 // lower triangle. 157 // 158 // subSymMatrix(int i1, int i2, int istep) 159 // Returns the SymMatrix which runs from i1 to i2 along the diagonal 160 // (not including i2) with an optional step, and includes the 161 // off diagonal in the same rows/cols. 162 // 163 // For example, with a SymMatrix of size 10, the x's below 164 // are the original data, the O's are the subSymMatrix returned 165 // with the command subSymMatrix(3,11,2), and the #'s are the 166 // subSymMatrix returned with subSymMatrix(0,3) 167 // 168 // ###xxxxxxx 169 // ###xxxxxxx 170 // ###xxxxxxx 171 // xxxOxOxOxO 172 // xxxxxxxxxx 173 // xxxOxOxOxO 174 // xxxxxxxxxx 175 // xxxOxOxOxO 176 // xxxxxxxxxx 177 // xxxOxOxOxO 178 // 179 // transpose(m) 180 // adjoint(m) 181 // conjugate(m) 182 // 183 // 184 // Functions of Matrixs: 185 // 186 // m.det() or Det(m) 187 // m.logDet() or m.logDet(T* sign) or LogDet(m) 188 // m.trace() or Trace(m) 189 // m.sumElements() or SumElements(m) 190 // m.sumAbsElements() or SumAbsElements(m) 191 // m.sumAbs2Elements() or SumAbs2Elements(m) 192 // m.norm() or m.normF() or Norm(m) or NormF(m) 193 // m.normSq() or NormSq(m) 194 // m.norm1() or Norm1(m) 195 // m.norm2() or Norm2(m) 196 // m.normInf() or NormInf(m) 197 // m.maxAbsElement() or MaxAbsElements(m) 198 // m.maxAbs2Element() or MaxAbs2Elements(m) 199 // 200 // m.inverse() or Inverse(m) 201 // m.makeInverse(minv) // Takes either a SymMatrix or Matrix argument 202 // m.makeInverseATA(invata) 203 // 204 // 205 // I/O: 206 // 207 // os << m 208 // Writes m to ostream os in the usual Matrix format 209 // 210 // os << CompactIO() << m 211 // Writes m to ostream os in the following compact format: 212 // size 213 // ( m(0,0) ) 214 // ( m(1,0) m(1,1) ) 215 // ... 216 // ( m(size-1,0) ... m(size-1,size-1) ) 217 // 218 // is >> m 219 // is >> CompactIO() >> m 220 // Reads m from istream is in either format 221 // 222 // 223 // Division Control Functions: 224 // 225 // m.divideUsing(dt) 226 // where dt is LU, CH, or SV 227 // 228 // m.lud(), m.chd(), m.svd(), and m.symsvd() return the 229 // corresponding Divider classes. 230 // 231 // For SymMatrixes, LU actually does an LDLT decomposition. 232 // ie. L(Lower Triangle) * Diagonal * L.transpose() 233 // For non-positive definite matrices, D may have 2x2 blocks along 234 // the diagonal, which can then be further decomposed via normal LU 235 // decomposition. 236 // 237 // The option unique to hermitian matrixes is CH - Cholskey decomposition. 238 // This is only appropriate if you know that your HermMatrix is 239 // positive definite (ie. all eigenvalues are positive). 240 // (This is guaranteed, for example, if all the square 241 // submatrices have positive determinant.) 242 // In this case, the SymMatrix can be decomposed into L*L.adjoint(). 243 // 244 // Finally, a difference for the SV Decomposition for Hermitian matrices 245 // is that the decomposition can be done with V = Ut. In this case, 246 // the "singular values" are really the eigenvalues of A. 247 // The only caveat is that they may be negative, whereas the usual 248 // definition of singular values is that they are all positive. 249 // Requiring positivity would destroy V = Ut, so it seemed more 250 // useful to leave them as the actual eigenvalues. Just keep that 251 // in mind if you use the singular values for anything that expects 252 // them to be positive. 253 // 254 // If m is complex, symmetric (i.e. not hermitian), then the SVD 255 // does not result in an eigenvalue decomposition, and we actually 256 // just do a regular SVD. In this case, the access is through 257 // symsvd(), rather than the usual svd() method. 258 259 260 #ifndef TMV_SymMatrix_H 261 #define TMV_SymMatrix_H 262 263 #include "tmv/TMV_BaseSymMatrix.h" 264 #include "tmv/TMV_BaseTriMatrix.h" 265 #include "tmv/TMV_BaseDiagMatrix.h" 266 #include "tmv/TMV_SymMatrixArithFunc.h" 267 #include "tmv/TMV_TriMatrix.h" 268 #include "tmv/TMV_Matrix.h" 269 #include "tmv/TMV_DiagMatrix.h" 270 #include "tmv/TMV_Array.h" 271 #include <vector> 272 273 namespace tmv { 274 275 template <typename T, typename T1, typename T2> 276 class OProdVV; 277 template <typename T, typename T1, typename T2> 278 class ProdMM; 279 template <typename T, typename T1, typename T2> 280 class ProdUL; 281 template <typename T, typename T1, typename T2> 282 class ProdLU; 283 284 template <typename T> 285 struct SymCopyHelper // real 286 { typedef SymMatrix<T> type; }; 287 template <typename T> 288 struct SymCopyHelper<std::complex<T> > // complex 289 // Have to copy to matrix, since don't know whether herm or sym. 290 { typedef Matrix<std::complex<T> > type; }; 291 292 template <typename T> 293 class GenSymMatrix : 294 virtual public AssignableToSymMatrix<T>, 295 public BaseMatrix<T>, 296 public DivHelper<T> 297 { 298 public: 299 300 typedef TMV_RealType(T) RT; 301 typedef TMV_ComplexType(T) CT; 302 typedef T value_type; 303 typedef RT real_type; 304 typedef CT complex_type; 305 typedef GenSymMatrix<T> type; 306 typedef typename SymCopyHelper<T>::type copy_type; 307 typedef ConstSymMatrixView<T> const_view_type; 308 typedef const_view_type const_transpose_type; 309 typedef const_view_type const_conjugate_type; 310 typedef const_view_type const_adjoint_type; 311 typedef ConstVectorView<T> const_vec_type; 312 typedef ConstMatrixView<T> const_rec_type; 313 typedef ConstUpperTriMatrixView<T> const_uppertri_type; 314 typedef ConstLowerTriMatrixView<T> const_lowertri_type; 315 typedef ConstSymMatrixView<RT> const_realpart_type; 316 typedef SymMatrixView<T> nonconst_type; 317 318 // 319 // Constructors 320 // 321 322 inline GenSymMatrix() {} 323 inline GenSymMatrix(const type&) {} 324 virtual inline ~GenSymMatrix() {} 325 326 // 327 // Access Functions 328 // 329 330 using AssignableToSymMatrix<T>::size; 331 inline ptrdiff_t colsize() const { return size(); } 332 inline ptrdiff_t rowsize() const { return size(); } 333 using AssignableToSymMatrix<T>::sym; 334 335 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 336 { 337 TMVAssert(i>=0 && i<size()); 338 TMVAssert(j>=0 && j<size()); 339 return cref(i,j); 340 } 341 342 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 343 { 344 TMVAssert(i>=0 && i<size()); 345 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()); 346 if ((i-j1<=0 && uplo()==Upper) || (j2-i<=1 && uplo()==Lower)) 347 return const_vec_type( 348 cptr()+i*stepi()+j1*stepj(),j2-j1,stepj(), 349 ct()); 350 else 351 return const_vec_type( 352 cptr()+i*stepj()+j1*stepi(),j2-j1,stepi(), 353 issym()?ct():TMV_ConjOf(T,ct())); 354 } 355 356 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 357 { 358 TMVAssert(j>=0 && j<size()); 359 TMVAssert(i1>=0 && i1-i2<=0 && i2<=size()); 360 if ((i2-j<=1 && uplo()==Upper) || (j-i1<=0 && uplo()==Lower)) 361 return const_vec_type( 362 cptr()+i1*stepi()+j*stepj(),i2-i1,stepi(), 363 ct()); 364 else 365 return const_vec_type( 366 cptr()+i1*stepj()+j*stepi(),i2-i1,stepj(), 367 issym()?ct():TMV_ConjOf(T,ct())); 368 } 369 370 inline const_vec_type diag() const 371 { return const_vec_type( 372 cptr(),size(),stepi()+stepj(),ct()); } 373 374 inline const_vec_type diag(ptrdiff_t i) const 375 { 376 TMVAssert(i>=-size() && i<=size()); 377 if (i>=0) 378 if (uplo()==Upper) 379 return const_vec_type( 380 cptr()+i*stepj(),size()-i, 381 stepi()+stepj(),ct()); 382 else 383 return const_vec_type( 384 cptr()+i*stepi(),size()-i, 385 stepi()+stepj(),issym()?ct():TMV_ConjOf(T,ct())); 386 else 387 if (uplo()==Upper) 388 return const_vec_type( 389 cptr()-i*stepj(),size()+i, 390 stepi()+stepj(),issym()?ct():TMV_ConjOf(T,ct())); 391 else 392 return const_vec_type( 393 cptr()-i*stepi(),size()+i, 394 stepi()+stepj(),ct()); 395 } 396 397 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 398 { 399 TMVAssert(i>=-size() && i<=size()); 400 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()-std::abs(i)); 401 const ptrdiff_t ds = stepi()+stepj(); 402 if (i>=0) 403 if (uplo()==Upper) 404 return const_vec_type( 405 cptr()+i*stepj()+j1*ds,j2-j1,ds,ct()); 406 else 407 return const_vec_type( 408 cptr()+i*stepi()+j1*ds,j2-j1,ds, 409 issym()?ct():TMV_ConjOf(T,ct())); 410 else 411 if (uplo()==Upper) 412 return const_vec_type( 413 cptr()-i*stepj()+j1*ds,j2-j1,ds, 414 issym()?ct():TMV_ConjOf(T,ct())); 415 else 416 return const_vec_type( 417 cptr()-i*stepi()+j1*ds,j2-j1,ds,ct()); 418 } 419 420 template <typename T2> 421 inline bool isSameAs(const BaseMatrix<T2>& ) const 422 { return false; } 423 424 inline bool isSameAs(const type& m2) const 425 { 426 if (this == &m2) return true; 427 else if (cptr()==m2.cptr() && size()==m2.size() && 428 (isReal(T()) || sym() == m2.sym())) { 429 if (uplo() == m2.uplo()) 430 return (stepi() == m2.stepi() && stepj() == m2.stepj() 431 && ct() == m2.ct()); 432 else 433 return (stepi() == m2.stepj() && stepj() == m2.stepi() 434 && issym() == (ct()==m2.ct())); 435 } else return false; 436 } 437 438 inline void assignToM(MatrixView<RT> m2) const 439 { 440 TMVAssert(isReal(T())); 441 TMVAssert(m2.colsize() == size()); 442 TMVAssert(m2.rowsize() == size()); 443 assignToS(SymMatrixViewOf(m2,Upper)); 444 if (size() > 0) 445 m2.lowerTri().offDiag() = 446 m2.upperTri().offDiag().transpose(); 447 } 448 449 inline void assignToM(MatrixView<CT> m2) const 450 { 451 TMVAssert(m2.colsize() == size()); 452 TMVAssert(m2.rowsize() == size()); 453 if (issym()) { 454 assignToS(SymMatrixViewOf(m2,Upper)); 455 if (size() > 0) 456 m2.lowerTri().offDiag() = 457 m2.upperTri().offDiag().transpose(); 458 } else { 459 m2.diag().imagPart().setZero(); 460 assignToS(HermMatrixViewOf(m2,Upper)); 461 if (size() > 0) 462 m2.lowerTri().offDiag() = 463 m2.upperTri().offDiag().adjoint(); 464 } 465 } 466 467 inline void assignToS(SymMatrixView<RT> m2) const 468 { 469 TMVAssert(isReal(T())); 470 TMVAssert(m2.size() == size()); 471 if (!isSameAs(m2)) m2.upperTri() = upperTri(); 472 } 473 474 inline void assignToS(SymMatrixView<CT> m2) const 475 { 476 TMVAssert(m2.size() == size()); 477 TMVAssert(isReal(T()) || m2.sym() == sym()); 478 if (!isSameAs(m2)) m2.upperTri() = upperTri(); 479 } 480 481 // 482 // subMatrix 483 // 484 485 bool hasSubMatrix( 486 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const; 487 488 inline const_rec_type subMatrix( 489 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 490 { 491 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 492 if ( (uplo()==Upper && i2-j1<=1) || (uplo()==Lower && j2-i1<=1) ) 493 return const_rec_type( 494 cptr()+i1*stepi()+j1*stepj(), 495 i2-i1,j2-j1,stepi(),stepj(),ct()); 496 else 497 return const_rec_type( 498 cptr()+i1*stepj()+j1*stepi(), 499 i2-i1,j2-j1,stepj(),stepi(), 500 issym()?ct():TMV_ConjOf(T,ct())); 501 } 502 503 inline const_rec_type subMatrix( 504 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 505 { 506 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 507 if ( (uplo()==Upper && i2-j1<=istep) || 508 (uplo()==Lower && j2-i1<=istep) ) 509 return const_rec_type( 510 cptr()+i1*stepi()+j1*stepj(), 511 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(), 512 ct()); 513 else 514 return const_rec_type( 515 cptr()+i1*stepj()+j1*stepi(), 516 (i2-i1)/istep,(j2-j1)/jstep,istep*stepj(),jstep*stepi(), 517 issym()?ct():TMV_ConjOf(T,ct())); 518 } 519 520 bool hasSubVector( 521 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const; 522 523 inline const_vec_type subVector( 524 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 525 { 526 TMVAssert(hasSubVector(i,j,istep,jstep,n)); 527 if ((i-j<=0 && uplo()==Upper) || (j-i<=0 && uplo()==Lower)) 528 return const_vec_type( 529 cptr()+i*stepi()+j*stepj(),n, 530 istep*stepi()+jstep*stepj(),ct()); 531 else 532 return const_vec_type( 533 cptr()+i*stepj()+j*stepi(),n,istep*stepj()+jstep*stepi(), 534 issym()?ct():TMV_ConjOf(T,ct())); 535 } 536 537 bool hasSubSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const; 538 539 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) const 540 { 541 TMVAssert(hasSubSymMatrix(i1,i2,1)); 542 return const_view_type( 543 cptr()+i1*(stepi()+stepj()),i2-i1, 544 stepi(),stepj(),sym(),uplo(),ct()); 545 } 546 547 inline const_view_type subSymMatrix( 548 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 549 { 550 TMVAssert(hasSubSymMatrix(i1,i2,istep)); 551 return const_view_type( 552 cptr()+i1*(stepi()+stepj()), 553 (i2-i1)/istep,istep*stepi(),istep*stepj(),sym(),uplo(),ct()); 554 } 555 556 inline const_uppertri_type upperTri(DiagType dt = NonUnitDiag) const 557 { 558 if (uplo() == Upper) 559 return const_uppertri_type( 560 cptr(),size(),stepi(),stepj(),dt,ct()); 561 else 562 return const_uppertri_type( 563 cptr(),size(),stepj(),stepi(),dt, 564 issym()?ct():TMV_ConjOf(T,ct())); 565 } 566 567 inline const_uppertri_type unitUpperTri() const 568 { 569 if (uplo() == Upper) 570 return const_uppertri_type( 571 cptr(),size(),stepi(),stepj(),UnitDiag,ct()); 572 else 573 return const_uppertri_type( 574 cptr(),size(),stepj(),stepi(),UnitDiag, 575 issym()?ct():TMV_ConjOf(T,ct())); 576 } 577 578 inline const_lowertri_type lowerTri(DiagType dt = NonUnitDiag) const 579 { 580 if (uplo() == Lower) 581 return const_lowertri_type( 582 cptr(),size(),stepi(),stepj(),dt,ct()); 583 else 584 return const_lowertri_type( 585 cptr(),size(),stepj(),stepi(),dt, 586 issym()?ct():TMV_ConjOf(T,ct())); 587 } 588 589 inline const_lowertri_type unitLowerTri() const 590 { 591 if (uplo() == Lower) 592 return const_lowertri_type( 593 cptr(),size(), 594 stepi(),stepj(),UnitDiag,ct()); 595 else 596 return const_lowertri_type( 597 cptr(),size(),stepj(),stepi(),UnitDiag, 598 issym()?ct():TMV_ConjOf(T,ct())); 599 } 600 601 inline const_realpart_type realPart() const 602 { 603 return const_realpart_type( 604 reinterpret_cast<const RT*>(cptr()),size(), 605 isReal(T()) ? stepi() : 2*stepi(), 606 isReal(T()) ? stepj() : 2*stepj(), Sym, uplo(), NonConj); 607 } 608 609 inline const_realpart_type imagPart() const 610 { 611 TMVAssert(isComplex(T())); 612 TMVAssert(issym()); 613 // The imaginary part of a Hermitian matrix is anti-symmetric 614 return const_realpart_type( 615 reinterpret_cast<const RT*>(cptr())+1, 616 size(),2*stepi(),2*stepj(), Sym,uplo(),NonConj); 617 } 618 619 // 620 // Views 621 // 622 623 inline const_view_type view() const 624 { 625 return const_view_type( 626 cptr(),size(),stepi(),stepj(),sym(),uplo(),ct()); 627 } 628 629 inline const_view_type transpose() const 630 { 631 return const_view_type( 632 cptr(),size(),stepj(),stepi(),sym(),TMV_UTransOf(uplo()),ct()); 633 } 634 635 inline const_view_type conjugate() const 636 { 637 return const_view_type( 638 cptr(),size(),stepi(),stepj(),sym(),uplo(),TMV_ConjOf(T,ct())); 639 } 640 641 inline const_view_type adjoint() const 642 { 643 return const_view_type( 644 cptr(),size(),stepj(),stepi(),sym(),TMV_UTransOf(uplo()), 645 TMV_ConjOf(T,ct())); 646 } 647 648 inline nonconst_type nonConst() const 649 { 650 return SymMatrixView<T>( 651 const_cast<T*>(cptr()),size(),stepi(),stepj(),sym(),uplo(),ct() 652 TMV_FIRSTLAST1( 653 cptr(),row(colsize()-1,0,colsize()).end().getP())); 654 } 655 656 // 657 // Functions of Matrix 658 // 659 660 T det() const; 661 662 RT logDet(T* sign=0) const; 663 664 bool isSingular() const; 665 666 inline T trace() const 667 { return diag().sumElements(); } 668 669 T sumElements() const; 670 671 RT sumAbsElements() const; 672 673 RT sumAbs2Elements() const; 674 675 inline RT norm() const 676 { return normF(); } 677 678 RT normF() const; 679 680 RT normSq(const RT scale = RT(1)) const; 681 682 RT norm1() const; 683 684 inline RT normInf() const 685 { return norm1(); } 686 687 RT maxAbsElement() const; 688 RT maxAbs2Element() const; 689 690 RT doNorm2() const; 691 inline RT norm2() const 692 { 693 if (this->divIsSet() && this->getDivType() == SV) 694 return DivHelper<T>::norm2(); 695 else return doNorm2(); 696 } 697 698 RT doCondition() const; 699 inline RT condition() const 700 { 701 if (this->divIsSet() && this->getDivType() == SV) 702 return DivHelper<T>::condition(); 703 else return doCondition(); 704 } 705 706 template <typename T1> 707 void doMakeInverse(SymMatrixView<T1> sinv) const; 708 709 template <typename T1> 710 inline void makeInverse(SymMatrixView<T1> minv) const 711 { 712 TMVAssert(minv.size() == size()); 713 TMVAssert(isherm() == minv.isherm()); 714 TMVAssert(issym() == minv.issym()); 715 doMakeInverse(minv); 716 } 717 718 inline void makeInverse(MatrixView<T> minv) const 719 { DivHelper<T>::makeInverse(minv); } 720 721 template <typename T1> 722 inline void makeInverse(MatrixView<T1> minv) const 723 { DivHelper<T>::makeInverse(minv); } 724 725 template <typename T1, int A> 726 inline void makeInverse(Matrix<T1,A>& minv) const 727 { DivHelper<T>::makeInverse(minv); } 728 729 template <typename T1, int A> 730 inline void makeInverse(SymMatrix<T1,A>& sinv) const 731 { 732 TMVAssert(issym()); 733 makeInverse(sinv.view()); 734 } 735 736 template <typename T1, int A> 737 inline void makeInverse(HermMatrix<T1,A>& sinv) const 738 { 739 TMVAssert(isherm()); 740 makeInverse(sinv.view()); 741 } 742 743 QuotXS<T,T> QInverse() const; 744 inline QuotXS<T,T> inverse() const 745 { return QInverse(); } 746 747 // 748 // Division Control 749 // 750 751 void setDiv() const; 752 753 inline void divideUsing(DivType dt) const 754 { 755 TMVAssert(dt == CH || dt == LU || dt == SV); 756 TMVAssert(isherm() || dt != CH); 757 DivHelper<T>::divideUsing(dt); 758 } 759 760 inline const SymLDLDiv<T>& lud() const 761 { 762 divideUsing(LU); 763 setDiv(); 764 TMVAssert(this->getDiv()); 765 TMVAssert(divIsLUDiv()); 766 return static_cast<const SymLDLDiv<T>&>(*this->getDiv()); 767 } 768 769 inline const HermCHDiv<T>& chd() const 770 { 771 TMVAssert(isherm()); 772 divideUsing(CH); 773 setDiv(); 774 TMVAssert(this->getDiv()); 775 TMVAssert(divIsCHDiv()); 776 return static_cast<const HermCHDiv<T>&>(*this->getDiv()); 777 } 778 779 inline const HermSVDiv<T>& svd() const 780 { 781 TMVAssert(isherm()); 782 divideUsing(SV); 783 setDiv(); 784 TMVAssert(this->getDiv()); 785 TMVAssert(divIsHermSVDiv()); 786 return static_cast<const HermSVDiv<T>&>(*this->getDiv()); 787 } 788 789 inline const SymSVDiv<T>& symsvd() const 790 { 791 TMVAssert(isComplex(T())); 792 TMVAssert(issym()); 793 divideUsing(SV); 794 setDiv(); 795 TMVAssert(this->getDiv()); 796 TMVAssert(divIsSymSVDiv()); 797 return static_cast<const SymSVDiv<T>&>(*this->getDiv()); 798 } 799 800 801 // 802 // I/O 803 // 804 805 void write(const TMV_Writer& writer) const; 806 807 virtual const T* cptr() const = 0; 808 virtual ptrdiff_t stepi() const = 0; 809 virtual ptrdiff_t stepj() const = 0; 810 virtual UpLoType uplo() const = 0; 811 virtual ConjType ct() const = 0; 812 virtual inline bool isrm() const { return stepj() == 1; } 813 virtual inline bool iscm() const { return stepi() == 1; } 814 inline bool isconj() const 815 { 816 TMVAssert(isComplex(T()) || ct()==NonConj); 817 return isComplex(T()) && ct() == Conj; 818 } 819 inline bool isherm() const { return isReal(T()) || sym() == Herm; } 820 inline bool issym() const { return isReal(T()) || sym() == Sym; } 821 inline bool isupper() const { return uplo() == Upper; } 822 823 inline bool isHermOK() const 824 { 825 if (issym()) return true; 826 else return diag().imagPart().normInf() == RT(0); 827 } 828 829 virtual T cref(ptrdiff_t i, ptrdiff_t j) const; 830 831 protected : 832 833 inline const BaseMatrix<T>& getMatrix() const { return *this; } 834 835 private : 836 837 type& operator=(const type&); 838 839 bool divIsLUDiv() const; 840 bool divIsCHDiv() const; 841 bool divIsHermSVDiv() const; 842 bool divIsSymSVDiv() const; 843 844 }; // GenSymMatrix 845 846 template <typename T, int A> 847 class ConstSymMatrixView : public GenSymMatrix<T> 848 { 849 public : 850 851 typedef ConstSymMatrixView<T,A> type; 852 typedef GenSymMatrix<T> base; 853 854 inline ConstSymMatrixView(const type& rhs) : 855 itsm(rhs.itsm), itss(rhs.itss), itssi(rhs.itssi), itssj(rhs.itssj), 856 itssym(rhs.itssym), itsuplo(rhs.itsuplo), 857 itsct(rhs.itsct) 858 { TMVAssert(Attrib<A>::viewok); } 859 860 inline ConstSymMatrixView(const base& rhs) : 861 itsm(rhs.cptr()), itss(rhs.size()), 862 itssi(rhs.stepi()), itssj(rhs.stepj()), 863 itssym(rhs.sym()), itsuplo(rhs.uplo()), 864 itsct(rhs.ct()) 865 { TMVAssert(Attrib<A>::viewok); } 866 867 inline ConstSymMatrixView( 868 const T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 869 SymType _sym, UpLoType _uplo, ConjType _ct) : 870 itsm(_m), itss(_s), itssi(_si), itssj(_sj), 871 itssym(_sym), itsuplo(_uplo), itsct(_ct) 872 { TMVAssert(Attrib<A>::viewok); } 873 874 virtual inline ~ConstSymMatrixView() 875 { 876 #ifdef TMV_EXTRA_DEBUG 877 const_cast<const T*&>(itsm) = 0; 878 #endif 879 } 880 881 inline ptrdiff_t size() const { return itss; } 882 inline const T* cptr() const { return itsm; } 883 inline ptrdiff_t stepi() const { return itssi; } 884 inline ptrdiff_t stepj() const { return itssj; } 885 inline SymType sym() const { return itssym; } 886 inline UpLoType uplo() const { return itsuplo; } 887 inline ConjType ct() const { return itsct; } 888 889 protected : 890 891 const T*const itsm; 892 const ptrdiff_t itss; 893 const ptrdiff_t itssi; 894 const ptrdiff_t itssj; 895 896 const SymType itssym; 897 const UpLoType itsuplo; 898 const ConjType itsct; 899 900 private : 901 902 type& operator=(const type&); 903 904 }; // ConstSymMatrixView 905 906 template <typename T> 907 class ConstSymMatrixView<T,FortranStyle> : 908 public ConstSymMatrixView<T,CStyle> 909 { 910 public : 911 912 typedef TMV_RealType(T) RT; 913 typedef ConstSymMatrixView<T,FortranStyle> type; 914 typedef ConstSymMatrixView<T,CStyle> c_type; 915 typedef GenSymMatrix<T> base; 916 typedef ConstSymMatrixView<T,FortranStyle> const_view_type; 917 typedef const_view_type const_transpose_type; 918 typedef const_view_type const_conjugate_type; 919 typedef const_view_type const_adjoint_type; 920 typedef ConstVectorView<T,FortranStyle> const_vec_type; 921 typedef ConstMatrixView<T,FortranStyle> const_rec_type; 922 typedef ConstUpperTriMatrixView<T,FortranStyle> const_uppertri_type; 923 typedef ConstLowerTriMatrixView<T,FortranStyle> const_lowertri_type; 924 typedef ConstSymMatrixView<RT,FortranStyle> const_realpart_type; 925 typedef SymMatrixView<T,FortranStyle> nonconst_type; 926 927 inline ConstSymMatrixView(const type& rhs) : c_type(rhs) {} 928 929 inline ConstSymMatrixView(const c_type& rhs) : c_type(rhs) {} 930 931 inline ConstSymMatrixView(const base& rhs) : c_type(rhs) {} 932 933 inline ConstSymMatrixView( 934 const T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 935 SymType _sym, UpLoType _uplo, ConjType _ct) : 936 c_type(_m,_s,_si,_sj,_sym,_uplo,_ct) {} 937 938 virtual inline ~ConstSymMatrixView() {} 939 940 // 941 // Access Functions 942 // 943 944 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 945 { 946 TMVAssert(i>0 && i<=this->size()); 947 TMVAssert(j>0 && j<=this->size()); 948 return base::cref(i-1,j-1); 949 } 950 951 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 952 { 953 TMVAssert(i>0 && i<=this->size()); 954 TMVAssert(j1>0 && j1-j2<=0 && j2<=this->size()); 955 return base::row(i-1,j1-1,j2); 956 } 957 958 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 959 { 960 TMVAssert(j>0 && j<=this->size()); 961 TMVAssert(i1>0 && i1-i2<=0 && i2<=this->size()); 962 return base::col(j-1,i1-1,i2); 963 } 964 965 inline const_vec_type diag() const 966 { return base::diag(); } 967 968 inline const_vec_type diag(ptrdiff_t i) const 969 { return base::diag(i); } 970 971 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 972 { 973 TMVAssert(j1>0); 974 return base::diag(i,j1-1,j2); 975 } 976 977 // 978 // SubMatrix 979 // 980 981 bool hasSubMatrix( 982 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const; 983 984 bool hasSubVector(ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const; 985 986 bool hasSubSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const; 987 988 inline const_rec_type subMatrix( 989 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 990 { 991 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 992 return base::subMatrix(i1-1,i2,j1-1,j2); 993 } 994 995 inline const_rec_type subMatrix( 996 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 997 { 998 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 999 return base::subMatrix(i1-1,i2-1+istep,j1-1,j2-1+jstep, 1000 istep,jstep); 1001 } 1002 1003 inline const_vec_type subVector( 1004 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 1005 { 1006 TMVAssert(hasSubVector(i,j,istep,jstep,n)); 1007 return base::subVector(i-1,j-1,istep,jstep,n); 1008 } 1009 1010 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1011 { 1012 TMVAssert(hasSubSymMatrix(i1,i2,1)); 1013 return base::subSymMatrix(i1-1,i2); 1014 } 1015 1016 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1017 { 1018 TMVAssert(hasSubSymMatrix(i1,i2,istep)); 1019 return base::subSymMatrix(i1-1,i2-1+istep,istep); 1020 } 1021 1022 inline const_uppertri_type upperTri(DiagType dt = NonUnitDiag) const 1023 { return base::upperTri(dt); } 1024 1025 inline const_uppertri_type unitUpperTri() const 1026 { return base::unitUpperTri(); } 1027 1028 inline const_lowertri_type lowerTri(DiagType dt = NonUnitDiag) const 1029 { return base::lowerTri(dt); } 1030 1031 inline const_lowertri_type unitLowerTri() const 1032 { return base::unitLowerTri(); } 1033 1034 inline const_realpart_type realPart() const 1035 { return base::realPart(); } 1036 1037 inline const_realpart_type imagPart() const 1038 { return base::imagPart(); } 1039 1040 // 1041 // Views 1042 // 1043 1044 inline const_view_type view() const 1045 { return base::view(); } 1046 1047 inline const_view_type transpose() const 1048 { return base::transpose(); } 1049 1050 inline const_view_type conjugate() const 1051 { return base::conjugate(); } 1052 1053 inline const_view_type adjoint() const 1054 { return base::adjoint(); } 1055 1056 inline nonconst_type nonConst() const 1057 { return base::nonConst(); } 1058 1059 private : 1060 1061 type& operator=(const type&); 1062 1063 }; // FortranStyle ConstSymMatrixView 1064 1065 template <typename T, int A> 1066 class SymMatrixView : public GenSymMatrix<T> 1067 { 1068 public: 1069 1070 typedef TMV_RealType(T) RT; 1071 typedef TMV_ComplexType(T) CT; 1072 typedef SymMatrixView<T,A> type; 1073 typedef GenSymMatrix<T> base; 1074 typedef SymMatrixView<T,A> view_type; 1075 typedef view_type transpose_type; 1076 typedef view_type conjugate_type; 1077 typedef view_type adjoint_type; 1078 typedef VectorView<T,A> vec_type; 1079 typedef MatrixView<T,A> rec_type; 1080 typedef UpperTriMatrixView<T,A> uppertri_type; 1081 typedef LowerTriMatrixView<T,A> lowertri_type; 1082 typedef SymMatrixView<RT,A> realpart_type; 1083 typedef ConstSymMatrixView<T,A> const_view_type; 1084 typedef const_view_type const_transpose_type; 1085 typedef const_view_type const_conjugate_type; 1086 typedef const_view_type const_adjoint_type; 1087 typedef ConstVectorView<T,A> const_vec_type; 1088 typedef ConstMatrixView<T,A> const_rec_type; 1089 typedef ConstUpperTriMatrixView<T,A> const_uppertri_type; 1090 typedef ConstLowerTriMatrixView<T,A> const_lowertri_type; 1091 typedef ConstSymMatrixView<RT,A> const_realpart_type; 1092 typedef typename RefHelper<T>::reference reference; 1093 1094 // 1095 // Constructors 1096 // 1097 1098 inline SymMatrixView(const type& rhs) : 1099 itsm(rhs.itsm), itss(rhs.itss), itssi(rhs.itssi), itssj(rhs.itssj), 1100 itssym(rhs.sym()), itsuplo(rhs.uplo()), 1101 itsct(rhs.ct()) TMV_DEFFIRSTLAST(rhs._first,rhs._last) 1102 { TMVAssert(Attrib<A>::viewok); } 1103 1104 inline SymMatrixView( 1105 T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 1106 SymType _sym, UpLoType _uplo, ConjType _ct 1107 TMV_PARAMFIRSTLAST(T) ) : 1108 itsm(_m), itss(_s), itssi(_si), itssj(_sj), 1109 itssym(_sym), itsuplo(_uplo), itsct(_ct) 1110 TMV_DEFFIRSTLAST(_first,_last) 1111 { TMVAssert(Attrib<A>::viewok); } 1112 1113 virtual inline ~SymMatrixView() 1114 { 1115 TMV_SETFIRSTLAST(0,0); 1116 #ifdef TMV_EXTRA_DEBUG 1117 const_cast<const T*&>(itsm) = 0; 1118 #endif 1119 } 1120 1121 // 1122 // Op= 1123 // 1124 1125 inline type& operator=(const type& m2) 1126 { 1127 TMVAssert(size() == m2.size()); 1128 TMVAssert(isReal(T()) || m2.sym() == sym()); 1129 if (!this->isSameAs(m2)) upperTri() = m2.upperTri(); 1130 return *this; 1131 } 1132 1133 inline type& operator=(const GenSymMatrix<RT>& m2) 1134 { 1135 TMVAssert(size() == m2.size()); 1136 m2.assignToS(*this); 1137 return *this; 1138 } 1139 1140 inline type& operator=(const GenSymMatrix<CT>& m2) 1141 { 1142 TMVAssert(isComplex(T())); 1143 TMVAssert(size() == m2.size()); 1144 TMVAssert(m2.sym() == sym()); 1145 m2.assignToS(*this); 1146 return *this; 1147 } 1148 1149 template <typename T2> 1150 inline type& operator=(const GenSymMatrix<T2>& m2) 1151 { 1152 TMVAssert(size() == m2.size()); 1153 TMVAssert(isComplex(T()) || isReal(T2())); 1154 TMVAssert(isReal(T2()) || m2.sym() == sym()); 1155 if (!this->isSameAs(m2)) upperTri() = m2.upperTri(); 1156 return *this; 1157 } 1158 1159 inline type& operator=(const T& x) 1160 { 1161 TMVAssert(issym() || TMV_IMAG(x) == RT(0)); 1162 return setToIdentity(x); 1163 } 1164 1165 inline type& operator=(const AssignableToSymMatrix<RT>& m2) 1166 { 1167 TMVAssert(size() == m2.size()); 1168 m2.assignToS(view()); 1169 return *this; 1170 } 1171 1172 inline type& operator=(const AssignableToSymMatrix<CT>& m2) 1173 { 1174 TMVAssert(isComplex(T())); 1175 TMVAssert(size() == m2.size()); 1176 TMVAssert(sym() == m2.sym()); 1177 m2.assignToS(view()); 1178 return *this; 1179 } 1180 1181 inline type& operator=(const GenDiagMatrix<RT>& m2) 1182 { 1183 TMVAssert(size() == m2.size()); 1184 m2.assignToD(DiagMatrixViewOf(diag())); 1185 upperTri().offDiag().setZero(); 1186 return *this; 1187 } 1188 1189 inline type& operator=(const GenDiagMatrix<CT>& m2) 1190 { 1191 TMVAssert(isComplex(T())); 1192 TMVAssert(size() == m2.size()); 1193 TMVAssert(issym()); 1194 m2.assignToD(DiagMatrixViewOf(diag())); 1195 upperTri().offDiag().setZero(); 1196 TMVAssert(this->isHermOK()); 1197 return *this; 1198 } 1199 1200 template <typename T2> 1201 inline type& operator=(const OProdVV<RT,T2,T2>& opvv) 1202 { 1203 TMVAssert(size() == opvv.colsize()); 1204 TMVAssert(size() == opvv.rowsize()); 1205 TMVAssert(opvv.getV1().isSameAs( 1206 issym() ? opvv.getV2().view() : opvv.getV2().conjugate())); 1207 TMVAssert(issym() || TMV_IMAG(opvv.getX()) == RT(0)); 1208 Rank1Update<false>(opvv.getX(),opvv.getV1(),*this); 1209 return *this; 1210 } 1211 1212 template <typename T2> 1213 inline type& operator=(const OProdVV<CT,T2,T2>& opvv) 1214 { 1215 TMVAssert(isComplex(T())); 1216 TMVAssert(size() == opvv.colsize()); 1217 TMVAssert(size() == opvv.rowsize()); 1218 TMVAssert(opvv.getV1().isSameAs( 1219 issym() ? opvv.getV2().view() : opvv.getV2().conjugate())); 1220 TMVAssert(issym() || TMV_IMAG(opvv.getX()) == RT(0)); 1221 Rank1Update<false>(opvv.getX(),opvv.getV1(),*this); 1222 return *this; 1223 } 1224 1225 template <typename T2> 1226 inline type& operator=(const ProdMM<RT,T2,T2>& pmm) 1227 { 1228 TMVAssert(size() == pmm.colsize()); 1229 TMVAssert(size() == pmm.rowsize()); 1230 TMVAssert(pmm.getM1().isSameAs(issym() ? pmm.getM2().transpose() : 1231 pmm.getM2().adjoint())); 1232 TMVAssert(issym() || TMV_IMAG(pmm.getX()) == RT(0)); 1233 RankKUpdate<false>(pmm.getX(),pmm.getM1(),*this); 1234 return *this; 1235 } 1236 1237 template <typename T2> 1238 inline type& operator=(const ProdMM<CT,T2,T2>& pmm) 1239 { 1240 TMVAssert(isComplex(T())); 1241 TMVAssert(size() == pmm.colsize()); 1242 TMVAssert(size() == pmm.rowsize()); 1243 TMVAssert(pmm.getM1().isSameAs( 1244 issym() ? pmm.getM2().transpose() : pmm.getM2().adjoint())); 1245 TMVAssert(issym() || TMV_IMAG(pmm.getX()) == RT(0)); 1246 RankKUpdate<false>(pmm.getX(),pmm.getM1(),*this); 1247 return *this; 1248 } 1249 1250 template <typename T2> 1251 inline type& operator=(const ProdUL<RT,T2,T2>& pmm) 1252 { 1253 TMVAssert(size() == pmm.colsize()); 1254 TMVAssert(size() == pmm.rowsize()); 1255 TMVAssert(pmm.getM1().isSameAs( 1256 issym() ? pmm.getM2().transpose() : pmm.getM2().adjoint())); 1257 TMVAssert(issym() || TMV_IMAG(pmm.getX()) == RT(0)); 1258 RankKUpdate<false>(pmm.getX(),pmm.getM1(),*this); 1259 return *this; 1260 } 1261 1262 template <typename T2> 1263 inline type& operator=(const ProdUL<CT,T2,T2>& pmm) 1264 { 1265 TMVAssert(isComplex(T())); 1266 TMVAssert(size() == pmm.colsize()); 1267 TMVAssert(size() == pmm.rowsize()); 1268 TMVAssert(pmm.getM1().isSameAs( 1269 issym() ? pmm.getM2().transpose() : pmm.getM2().adjoint())); 1270 TMVAssert(issym() || TMV_IMAG(pmm.getX()) == RT(0)); 1271 RankKUpdate<false>(pmm.getX(),pmm.getM1(),*this); 1272 return *this; 1273 } 1274 1275 template <typename T2> 1276 inline type& operator=(const ProdLU<RT,T2,T2>& pmm) 1277 { 1278 TMVAssert(size() == pmm.colsize()); 1279 TMVAssert(size() == pmm.rowsize()); 1280 TMVAssert(pmm.getM1().isSameAs( 1281 issym() ? pmm.getM2().transpose() : pmm.getM2().adjoint())); 1282 TMVAssert(issym() || TMV_IMAG(pmm.getX()) == RT(0)); 1283 RankKUpdate<false>(pmm.getX(),pmm.getM1(),*this); 1284 return *this; 1285 } 1286 1287 template <typename T2> 1288 inline type& operator=(const ProdLU<CT,T2,T2>& pmm) 1289 { 1290 TMVAssert(isComplex(T())); 1291 TMVAssert(size() == pmm.colsize()); 1292 TMVAssert(size() == pmm.rowsize()); 1293 TMVAssert(pmm.getM1().isSameAs( 1294 issym() ? pmm.getM2().transpose() : pmm.getM2().adjoint())); 1295 TMVAssert(issym() || TMV_IMAG(pmm.getX()) == RT(0)); 1296 RankKUpdate<false>(pmm.getX(),pmm.getM1(),*this); 1297 return *this; 1298 } 1299 1300 1301 // 1302 // Access 1303 // 1304 1305 inline reference operator()(ptrdiff_t i,ptrdiff_t j) 1306 { 1307 TMVAssert(i>=0 && i<size()); 1308 TMVAssert(j>=0 && j<size()); 1309 return ref(i,j); 1310 } 1311 1312 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 1313 { 1314 TMVAssert(i>=0 && i<size()); 1315 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()); 1316 if ((i-j1<=0 && uplo()==Upper) || (j2-i<=1 && uplo()==Lower)) 1317 return vec_type( 1318 ptr()+i*stepi()+j1*stepj(),j2-j1,stepj(), 1319 ct() TMV_FIRSTLAST); 1320 else 1321 return vec_type( 1322 ptr()+i*stepj()+j1*stepi(),j2-j1,stepi(), 1323 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1324 } 1325 1326 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 1327 { 1328 TMVAssert(j>=0 && j<size()); 1329 TMVAssert(i1>=0 && i1-i2<=0 && i2<=size()); 1330 if ((i2-j<=1 && uplo()==Upper) || (j-i1<=0 && uplo()==Lower)) 1331 return vec_type( 1332 ptr()+i1*stepi()+j*stepj(),i2-i1,stepi(), 1333 ct() TMV_FIRSTLAST); 1334 else 1335 return vec_type( 1336 ptr()+i1*stepj()+j*stepi(),i2-i1,stepj(), 1337 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1338 } 1339 1340 inline vec_type diag() 1341 { return vec_type(ptr(),size(),stepi()+stepj(),ct() TMV_FIRSTLAST); } 1342 1343 inline vec_type diag(ptrdiff_t i) 1344 { 1345 TMVAssert(i>=-size() && i<=size()); 1346 if (i>=0) 1347 if (uplo()==Upper) 1348 return vec_type( 1349 ptr()+i*stepj(),size()-i, 1350 stepi()+stepj(),ct() TMV_FIRSTLAST); 1351 else 1352 return vec_type( 1353 ptr()+i*stepi(),size()-i, 1354 stepi()+stepj(),this->issym()?ct():TMV_ConjOf(T,ct()) 1355 TMV_FIRSTLAST); 1356 else 1357 if (uplo()==Upper) 1358 return vec_type( 1359 ptr()-i*stepj(),size()+i, 1360 stepi()+stepj(),this->issym()?ct():TMV_ConjOf(T,ct()) 1361 TMV_FIRSTLAST); 1362 else 1363 return vec_type( 1364 ptr()-i*stepi(),size()+i, 1365 stepi()+stepj(),ct() TMV_FIRSTLAST); 1366 } 1367 1368 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 1369 { 1370 TMVAssert(i>=-size() && i<=size()); 1371 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()-std::abs(i)); 1372 const ptrdiff_t ds = stepi()+stepj(); 1373 if (i>=0) 1374 if (uplo()==Upper) 1375 return vec_type( 1376 ptr()+i*stepj()+j1*ds,j2-j1,ds,ct() 1377 TMV_FIRSTLAST); 1378 else 1379 return vec_type( 1380 ptr()+i*stepi()+j1*ds,j2-j1,ds, 1381 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1382 else 1383 if (uplo()==Upper) 1384 return vec_type( 1385 ptr()-i*stepj()+j1*ds,j2-j1,ds, 1386 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1387 else 1388 return vec_type( 1389 ptr()-i*stepi()+j1*ds,j2-j1,ds,ct() 1390 TMV_FIRSTLAST); 1391 } 1392 1393 1394 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 1395 { return base::operator()(i,j); } 1396 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1397 { return base::row(i,j1,j2); } 1398 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 1399 { return base::col(j,i1,i2); } 1400 inline const_vec_type diag() const 1401 { return base::diag(); } 1402 inline const_vec_type diag(ptrdiff_t i) const 1403 { return base::diag(i); } 1404 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1405 { return base::diag(i,j1,j2); } 1406 1407 // 1408 // Modifying Functions 1409 // 1410 1411 inline type& setZero() 1412 { upperTri().setZero(); return *this; } 1413 1414 inline type& setAllTo(const T& x) 1415 { 1416 TMVAssert(TMV_IMAG(x)==RT(0) || this->issym()); 1417 upperTri().setAllTo(x); return *this; 1418 } 1419 1420 inline type& addToAll(const T& x) 1421 { 1422 TMVAssert(TMV_IMAG(x)==RT(0) || this->issym()); 1423 upperTri().addToAll(x); return *this; 1424 } 1425 1426 inline type& clip(RT thresh) 1427 { upperTri().clip(thresh); return *this; } 1428 1429 inline type& transposeSelf() 1430 { if (!this->issym()) upperTri().conjugateSelf(); return *this; } 1431 1432 inline type& conjugateSelf() 1433 { if (isComplex(T())) upperTri().conjugateSelf(); return *this; } 1434 1435 inline type& setToIdentity(const T& x=T(1)) 1436 { 1437 TMVAssert(TMV_IMAG(x)==RT(0) || this->issym()); 1438 setZero(); diag().setAllTo(x); return *this; 1439 } 1440 1441 type& swapRowsCols(ptrdiff_t i1, ptrdiff_t i2); 1442 1443 type& permuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2); 1444 1445 type& reversePermuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2); 1446 1447 inline type& permuteRowsCols(const ptrdiff_t* p) 1448 { return permuteRowsCols(p,0,size()); } 1449 1450 inline type& reversePermuteRowsCols(const ptrdiff_t* p) 1451 { return reversePermuteRowsCols(p,0,size()); } 1452 1453 // 1454 // subMatrix 1455 // 1456 1457 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 1458 { 1459 TMVAssert(base::hasSubMatrix(i1,i2,j1,j2,1,1)); 1460 if ( (uplo()==Upper && i2-j1<=1) || 1461 (uplo()==Lower && j2-i1<=1) ) 1462 return rec_type( 1463 ptr()+i1*stepi()+j1*stepj(), 1464 i2-i1,j2-j1,stepi(),stepj(),ct() TMV_FIRSTLAST); 1465 else 1466 return rec_type( 1467 ptr()+i1*stepj()+j1*stepi(),i2-i1,j2-j1,stepj(),stepi(), 1468 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1469 } 1470 1471 inline rec_type subMatrix( 1472 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 1473 { 1474 TMVAssert(base::hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 1475 if ( (uplo()==Upper && i2-j1<=istep) || 1476 (uplo()==Lower && j2-i1<=istep) ) 1477 return rec_type( 1478 ptr()+i1*stepi()+j1*stepj(), 1479 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(), 1480 ct() TMV_FIRSTLAST); 1481 else 1482 return rec_type( 1483 ptr()+i1*stepj()+j1*stepi(),(i2-i1)/istep,(j2-j1)/jstep, 1484 istep*stepj(),jstep*stepi(), 1485 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1486 } 1487 1488 inline vec_type subVector(ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) 1489 { 1490 TMVAssert(base::hasSubVector(i,j,istep,jstep,n)); 1491 if ((i-j<=0 && uplo()==Upper) || (j-i<=0 && uplo()==Lower)) 1492 return vec_type( 1493 ptr()+i*stepi()+j*stepj(),n, 1494 istep*stepi()+jstep*stepj(),ct() TMV_FIRSTLAST); 1495 else 1496 return vec_type( 1497 ptr()+i*stepj()+j*stepi(),n,istep*stepj()+jstep*stepi(), 1498 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1499 } 1500 1501 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) 1502 { 1503 TMVAssert(base::hasSubSymMatrix(i1,i2,1)); 1504 return view_type( 1505 ptr()+i1*(stepi()+stepj()),i2-i1, 1506 stepi(),stepj(),sym(),uplo(),ct() TMV_FIRSTLAST); 1507 } 1508 1509 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 1510 { 1511 TMVAssert(base::hasSubSymMatrix(i1,i2,istep)); 1512 return view_type( 1513 ptr()+i1*(stepi()+stepj()),(i2-i1)/istep, 1514 istep*stepi(),istep*stepj(),sym(),uplo(), ct() TMV_FIRSTLAST); 1515 } 1516 1517 inline uppertri_type upperTri(DiagType dt = NonUnitDiag) 1518 { 1519 if (uplo() == Upper) 1520 return uppertri_type( 1521 ptr(),size(),stepi(),stepj(),dt,ct() TMV_FIRSTLAST); 1522 else 1523 return uppertri_type( 1524 ptr(),size(),stepj(),stepi(),dt, 1525 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1526 } 1527 1528 inline uppertri_type unitUpperTri() 1529 { 1530 if (uplo() == Upper) 1531 return uppertri_type( 1532 ptr(),size(),stepi(),stepj(),UnitDiag,ct() TMV_FIRSTLAST); 1533 else 1534 return uppertri_type( 1535 ptr(),size(),stepj(),stepi(),UnitDiag, 1536 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1537 } 1538 1539 inline lowertri_type lowerTri(DiagType dt = NonUnitDiag) 1540 { 1541 if (uplo() == Lower) 1542 return lowertri_type( 1543 ptr(),size(),stepi(),stepj(),dt,ct() TMV_FIRSTLAST); 1544 else 1545 return lowertri_type( 1546 ptr(),size(),stepj(),stepi(),dt, 1547 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1548 } 1549 1550 inline lowertri_type unitLowerTri() 1551 { 1552 if (uplo() == Lower) 1553 return lowertri_type( 1554 ptr(),size(),stepi(),stepj(),UnitDiag,ct() TMV_FIRSTLAST); 1555 else 1556 return lowertri_type( 1557 ptr(),size(),stepj(),stepi(),UnitDiag, 1558 this->issym()?ct():TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1559 } 1560 1561 inline realpart_type realPart() 1562 { 1563 return realpart_type( 1564 reinterpret_cast<RT*>(ptr()),size(), 1565 isReal(T()) ? stepi() : 2*stepi(), 1566 isReal(T()) ? stepj() : 2*stepj(), sym(),uplo(),NonConj 1567 #ifdef TMVFLDEBUG 1568 ,reinterpret_cast<const RT*>(_first) 1569 ,reinterpret_cast<const RT*>(_last) 1570 #endif 1571 ); 1572 } 1573 1574 inline realpart_type imagPart() 1575 { 1576 TMVAssert(isComplex(T())); 1577 TMVAssert(this->issym()); 1578 return realpart_type( 1579 reinterpret_cast<RT*>(ptr())+1, 1580 size(),2*stepi(),2*stepj(),sym(),uplo(),NonConj 1581 #ifdef TMVFLDEBUG 1582 ,reinterpret_cast<const RT*>(_first)+1 1583 ,reinterpret_cast<const RT*>(_last)+1 1584 #endif 1585 ); 1586 } 1587 1588 inline view_type view() 1589 { return *this; } 1590 1591 inline view_type transpose() 1592 { 1593 return view_type( 1594 ptr(),size(),stepj(),stepi(),sym(),TMV_UTransOf(uplo()),ct() 1595 TMV_FIRSTLAST); 1596 } 1597 1598 inline view_type conjugate() 1599 { 1600 return view_type( 1601 ptr(),size(),stepi(),stepj(),sym(),uplo(),TMV_ConjOf(T,ct()) 1602 TMV_FIRSTLAST); 1603 } 1604 1605 inline view_type adjoint() 1606 { 1607 return view_type( 1608 ptr(),size(),stepj(),stepi(),sym(),TMV_UTransOf(uplo()), 1609 TMV_ConjOf(T,ct()) TMV_FIRSTLAST); 1610 } 1611 1612 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 1613 { return base::subMatrix(i1,i2,j1,j2); } 1614 inline const_rec_type subMatrix( 1615 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1616 { return base::subMatrix(i1,i2,j1,j2,istep,jstep); } 1617 inline const_vec_type subVector( 1618 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 1619 { return base::subVector(i,j,istep,jstep,n); } 1620 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) const 1621 { return base::subSymMatrix(i1,i2); } 1622 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1623 { return base::subSymMatrix(i1,i2,istep); } 1624 inline const_uppertri_type upperTri(DiagType dt = NonUnitDiag) const 1625 { return base::upperTri(dt); } 1626 inline const_uppertri_type unitUpperTri() const 1627 { return base::unitUpperTri(); } 1628 inline const_lowertri_type lowerTri(DiagType dt = NonUnitDiag) const 1629 { return base::lowerTri(dt); } 1630 inline const_lowertri_type unitLowerTri() const 1631 { return base::unitLowerTri(); } 1632 inline const_realpart_type realPart() const 1633 { return base::realPart(); } 1634 inline const_realpart_type imagPart() const 1635 { return base::imagPart(); } 1636 inline const_view_type view() const 1637 { return base::view(); } 1638 inline const_view_type transpose() const 1639 { return base::transpose(); } 1640 inline const_view_type conjugate() const 1641 { return base::conjugate(); } 1642 inline const_view_type adjoint() const 1643 { return base::adjoint(); } 1644 1645 // 1646 // I/O 1647 // 1648 1649 void read(const TMV_Reader& reader); 1650 1651 inline ptrdiff_t size() const { return itss; } 1652 inline const T* cptr() const { return itsm; } 1653 inline T* ptr() { return itsm; } 1654 inline ptrdiff_t stepi() const { return itssi; } 1655 inline ptrdiff_t stepj() const { return itssj; } 1656 inline SymType sym() const { return itssym; } 1657 inline UpLoType uplo() const { return itsuplo; } 1658 inline ConjType ct() const { return itsct; } 1659 using base::issym; 1660 using base::isherm; 1661 1662 reference ref(ptrdiff_t i, ptrdiff_t j); 1663 1664 protected : 1665 1666 T*const itsm; 1667 const ptrdiff_t itss; 1668 const ptrdiff_t itssi; 1669 const ptrdiff_t itssj; 1670 1671 const SymType itssym; 1672 const UpLoType itsuplo; 1673 const ConjType itsct; 1674 1675 #ifdef TMVFLDEBUG 1676 public: 1677 const T* _first; 1678 const T* _last; 1679 protected: 1680 #endif 1681 1682 }; // SymMatrixView 1683 1684 template <typename T> 1685 class SymMatrixView<T,FortranStyle> : public SymMatrixView<T,CStyle> 1686 { 1687 public: 1688 1689 typedef TMV_RealType(T) RT; 1690 typedef TMV_ComplexType(T) CT; 1691 typedef SymMatrixView<T,FortranStyle> type; 1692 typedef SymMatrixView<T,CStyle> c_type; 1693 typedef GenSymMatrix<T> base; 1694 typedef SymMatrixView<T,FortranStyle> view_type; 1695 typedef view_type transpose_type; 1696 typedef view_type conjugate_type; 1697 typedef view_type adjoint_type; 1698 typedef VectorView<T,FortranStyle> vec_type; 1699 typedef MatrixView<T,FortranStyle> rec_type; 1700 typedef UpperTriMatrixView<T,FortranStyle> uppertri_type; 1701 typedef LowerTriMatrixView<T,FortranStyle> lowertri_type; 1702 typedef SymMatrixView<RT,FortranStyle> realpart_type; 1703 typedef ConstSymMatrixView<T,FortranStyle> const_view_type; 1704 typedef const_view_type const_transpose_type; 1705 typedef const_view_type const_conjugate_type; 1706 typedef const_view_type const_adjoint_type; 1707 typedef ConstVectorView<T,FortranStyle> const_vec_type; 1708 typedef ConstMatrixView<T,FortranStyle> const_rec_type; 1709 typedef ConstUpperTriMatrixView<T,FortranStyle> const_uppertri_type; 1710 typedef ConstLowerTriMatrixView<T,FortranStyle> const_lowertri_type; 1711 typedef ConstSymMatrixView<RT,FortranStyle> const_realpart_type; 1712 typedef ConstSymMatrixView<T,FortranStyle> const_type; 1713 typedef typename RefHelper<T>::reference reference; 1714 1715 // 1716 // Constructors 1717 // 1718 1719 inline SymMatrixView(const type& rhs) : c_type(rhs) {} 1720 1721 inline SymMatrixView(const c_type& rhs) : c_type(rhs) {} 1722 1723 inline SymMatrixView( 1724 T* _m, ptrdiff_t _s, ptrdiff_t _si, ptrdiff_t _sj, 1725 SymType _sym, UpLoType _uplo, ConjType _ct 1726 TMV_PARAMFIRSTLAST(T) ) : 1727 c_type(_m,_s,_si,_sj,_sym,_uplo,_ct 1728 TMV_FIRSTLAST1(_first,_last) ) {} 1729 1730 virtual inline ~SymMatrixView() {} 1731 1732 // 1733 // Op= 1734 // 1735 1736 inline type& operator=(const type& m2) 1737 { c_type::operator=(m2); return *this; } 1738 1739 inline type& operator=(const c_type& m2) 1740 { c_type::operator=(m2); return *this; } 1741 1742 inline type& operator=(const GenSymMatrix<RT>& m2) 1743 { c_type::operator=(m2); return *this; } 1744 1745 inline type& operator=(const GenSymMatrix<CT>& m2) 1746 { c_type::operator=(m2); return *this; } 1747 1748 template <typename T2> 1749 inline type& operator=(const GenSymMatrix<T2>& m2) 1750 { c_type::operator=(m2); return *this; } 1751 1752 inline type& operator=(const T& x) 1753 { c_type::operator=(x); return *this; } 1754 1755 inline type& operator=(const AssignableToSymMatrix<RT>& m2) 1756 { c_type::operator=(m2); return *this; } 1757 1758 inline type& operator=(const AssignableToSymMatrix<CT>& m2) 1759 { c_type::operator=(m2); return *this; } 1760 1761 inline type& operator=(const GenDiagMatrix<RT>& m2) 1762 { c_type::operator=(m2); return *this; } 1763 1764 inline type& operator=(const GenDiagMatrix<CT>& m2) 1765 { c_type::operator=(m2); return *this; } 1766 1767 template <typename T2> 1768 inline type& operator=(const OProdVV<T,T2,T2>& opvv) 1769 { c_type::operator=(opvv); return *this; } 1770 1771 template <typename T2> 1772 inline type& operator=(const ProdMM<T,T2,T2>& pmm) 1773 { c_type::operator=(pmm); return *this; } 1774 1775 template <typename T2> 1776 inline type& operator=(const ProdLU<T,T2,T2>& pmm) 1777 { c_type::operator=(pmm); return *this; } 1778 1779 template <typename T2> 1780 inline type& operator=(const ProdUL<T,T2,T2>& pmm) 1781 { c_type::operator=(pmm); return *this; } 1782 1783 1784 // 1785 // Access 1786 // 1787 1788 inline reference operator()(ptrdiff_t i,ptrdiff_t j) 1789 { 1790 TMVAssert(i>0 && i<=this->size()); 1791 TMVAssert(j>0 && j<=this->size()); 1792 return c_type::ref(i-1,j-1); 1793 } 1794 1795 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 1796 { 1797 TMVAssert(i>0 && i<=this->size()); 1798 TMVAssert(j1>0 && j1-j2<=0 && j2<=this->size()); 1799 return c_type::row(i-1,j1-1,j2); 1800 } 1801 1802 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 1803 { 1804 TMVAssert(j>0 && j<=this->size()); 1805 TMVAssert(i1>0 && i1-i2<=0 && i2<=this->size()); 1806 return c_type::col(j-1,i1-1,i2); 1807 } 1808 1809 inline vec_type diag() 1810 { return c_type::diag(); } 1811 1812 inline vec_type diag(ptrdiff_t i) 1813 { return c_type::diag(i); } 1814 1815 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 1816 { 1817 TMVAssert(j1>0); 1818 return c_type::diag(i,j1-1,j2); 1819 } 1820 1821 1822 inline T operator()(ptrdiff_t i,ptrdiff_t j) const 1823 { 1824 TMVAssert(i>0 && i<=this->size()); 1825 TMVAssert(j>0 && j<=this->size()); 1826 return c_type::cref(i-1,j-1); 1827 } 1828 1829 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1830 { 1831 TMVAssert(i>0 && i<=this->size()); 1832 TMVAssert(j1>0 && j1-j2<=0 && j2<=this->size()); 1833 return c_type::row(i-1,j1-1,j2); 1834 } 1835 1836 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 1837 { 1838 TMVAssert(j>0 && j<=this->size()); 1839 TMVAssert(i1>0 && i1-i2<=0 && i2<=this->size()); 1840 return c_type::col(j-1,i1-1,i2); 1841 } 1842 1843 inline const_vec_type diag() const 1844 { return c_type::diag(); } 1845 1846 inline const_vec_type diag(ptrdiff_t i) const 1847 { return c_type::diag(i); } 1848 1849 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 1850 { 1851 TMVAssert(j1>0); 1852 return c_type::diag(i,j1-1,j2); 1853 } 1854 1855 // 1856 // Modifying Functions 1857 // 1858 1859 inline type& setZero() 1860 { c_type::setZero(); return *this; } 1861 1862 inline type& setAllTo(const T& x) 1863 { c_type::setAllTo(x); return *this; } 1864 1865 inline type& addToAll(const T& x) 1866 { c_type::addToAll(x); return *this; } 1867 1868 inline type& clip(RT thresh) 1869 { c_type::clip(thresh); return *this; } 1870 1871 inline type& conjugateSelf() 1872 { c_type::conjugateSelf(); return *this; } 1873 1874 inline type& transposeSelf() 1875 { c_type::transposeSelf(); return *this; } 1876 1877 inline type& setToIdentity(const T& x=T(1)) 1878 { c_type::setToIdentity(x); return *this; } 1879 1880 inline type& swapRowsCols(ptrdiff_t i1, ptrdiff_t i2) 1881 { 1882 TMVAssert(i1>0 && i1<=this->size()); 1883 TMVAssert(i2>0 && i2<=this->size()); 1884 c_type::swapRowsCols(i1-1,i2-1); 1885 return *this; 1886 } 1887 1888 inline type& permuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 1889 { 1890 TMVAssert(i1>0 && i1-i2<=0 && i2<=this->size()); 1891 c_type::permuteRowsCols(p,i1-1,i2); 1892 return *this; 1893 } 1894 1895 inline type& reversePermuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 1896 { 1897 TMVAssert(i1>0 && i1-i2<=0 && i2<=this->size()); 1898 c_type::reversePermuteRowsCols(p,i1-1,i2); 1899 return *this; 1900 } 1901 1902 inline type& permuteRowsCols(const ptrdiff_t* p) 1903 { c_type::permuteRowsCols(p); return *this; } 1904 1905 inline type& reversePermuteRowsCols(const ptrdiff_t* p) 1906 { c_type::reversePermuteRowsCols(p); return *this; } 1907 1908 // 1909 // subMatrix 1910 // 1911 1912 inline bool hasSubMatrix( 1913 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1914 { return const_type(*this).hasSubMatrix(i1,i2,j1,j2,istep,jstep); } 1915 1916 inline bool hasSubVector( 1917 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 1918 { return const_type(*this).hasSubVector(i,j,istep,jstep,n); } 1919 1920 1921 inline bool hasSubSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 1922 { return const_type(*this).hasSubSymMatrix(i1,i2,istep); } 1923 1924 inline rec_type subMatrix( 1925 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 1926 { 1927 TMVAssert(hasSubMatrix(i1,i2,j1,j2,1,1)); 1928 return c_type::subMatrix(i1-1,i2,j1-1,j2); 1929 } 1930 1931 inline rec_type subMatrix( 1932 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 1933 { 1934 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 1935 return c_type::subMatrix( 1936 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 1937 } 1938 1939 inline vec_type subVector(ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) 1940 { 1941 TMVAssert(hasSubVector(i,j,istep,jstep,n)); 1942 return c_type::subVector(i-1,j-1,istep,jstep,n); 1943 } 1944 1945 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) 1946 { 1947 TMVAssert(hasSubSymMatrix(i1,i2,1)); 1948 return c_type::subSymMatrix(i1-1,i2); 1949 } 1950 1951 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 1952 { 1953 TMVAssert(hasSubSymMatrix(i1,i2,istep)); 1954 return c_type::subSymMatrix(i1-1,i2-1+istep,istep); 1955 } 1956 1957 inline uppertri_type upperTri(DiagType dt = NonUnitDiag) 1958 { return c_type::upperTri(dt); } 1959 1960 inline uppertri_type unitUpperTri() 1961 { return c_type::unitUpperTri(); } 1962 1963 inline lowertri_type lowerTri(DiagType dt = NonUnitDiag) 1964 { return c_type::lowerTri(dt); } 1965 1966 inline lowertri_type unitLowerTri() 1967 { return c_type::unitLowerTri(); } 1968 1969 inline realpart_type realPart() 1970 { return c_type::realPart(); } 1971 1972 inline realpart_type imagPart() 1973 { return c_type::imagPart(); } 1974 1975 inline view_type view() 1976 { return c_type::view(); } 1977 1978 inline view_type transpose() 1979 { return c_type::transpose(); } 1980 1981 inline view_type conjugate() 1982 { return c_type::conjugate(); } 1983 1984 inline view_type adjoint() 1985 { return c_type::adjoint(); } 1986 1987 inline const_rec_type subMatrix( 1988 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 1989 { 1990 TMVAssert(hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 1991 return c_type::subMatrix( 1992 i1-1,i2-1+istep,j1-1,j2-1+jstep,istep,jstep); 1993 } 1994 1995 inline const_vec_type subVector( 1996 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 1997 { 1998 TMVAssert(hasSubVector(i,j,istep,jstep,n)); 1999 return c_type::subVector(i-1,j-1,istep,jstep,n); 2000 } 2001 2002 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) const 2003 { 2004 TMVAssert(hasSubSymMatrix(i1,i2,1)); 2005 return c_type::subSymMatrix(i1-1,i2); 2006 } 2007 2008 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 2009 { 2010 TMVAssert(hasSubSymMatrix(i1,i2,istep)); 2011 return c_type::subSymMatrix(i1-1,i2-1+istep,istep); 2012 } 2013 2014 inline const_uppertri_type upperTri(DiagType dt = NonUnitDiag) const 2015 { return c_type::upperTri(dt); } 2016 2017 inline const_uppertri_type unitUpperTri() const 2018 { return c_type::unitUpperTri(); } 2019 2020 inline const_lowertri_type lowerTri(DiagType dt = NonUnitDiag) const 2021 { return c_type::lowerTri(dt); } 2022 2023 inline const_lowertri_type unitLowerTri() const 2024 { return c_type::unitLowerTri(); } 2025 2026 inline const_realpart_type realPart() const 2027 { return c_type::realPart(); } 2028 2029 inline const_realpart_type imagPart() const 2030 { return c_type::imagPart(); } 2031 2032 inline const_view_type view() const 2033 { return c_type::view(); } 2034 2035 inline const_view_type transpose() const 2036 { return c_type::transpose(); } 2037 2038 inline const_view_type conjugate() const 2039 { return c_type::conjugate(); } 2040 2041 inline const_view_type adjoint() const 2042 { return c_type::adjoint(); } 2043 2044 }; // FortranStyle SymMatrixView 2045 2046 template <typename T, int A> 2047 class SymMatrix : public GenSymMatrix<T> 2048 { 2049 public: 2050 2051 enum { S = A & AllStorageType }; 2052 enum { U = A & Upper }; 2053 enum { I = A & FortranStyle }; 2054 2055 typedef TMV_RealType(T) RT; 2056 typedef TMV_ComplexType(T) CT; 2057 typedef SymMatrix<T,A> type; 2058 typedef type copy_type; 2059 typedef GenSymMatrix<T> base; 2060 typedef ConstSymMatrixView<T,I> const_view_type; 2061 typedef const_view_type const_transpose_type; 2062 typedef const_view_type const_conjugate_type; 2063 typedef const_view_type const_adjoint_type; 2064 typedef ConstVectorView<T,I> const_vec_type; 2065 typedef ConstMatrixView<T,I> const_rec_type; 2066 typedef ConstUpperTriMatrixView<T,I> const_uppertri_type; 2067 typedef ConstLowerTriMatrixView<T,I> const_lowertri_type; 2068 typedef ConstSymMatrixView<RT,I> const_realpart_type; 2069 typedef SymMatrixView<T,I> view_type; 2070 typedef view_type transpose_type; 2071 typedef view_type conjugate_type; 2072 typedef view_type adjoint_type; 2073 typedef VectorView<T,I> vec_type; 2074 typedef MatrixView<T,I> rec_type; 2075 typedef UpperTriMatrixView<T,I> uppertri_type; 2076 typedef LowerTriMatrixView<T,I> lowertri_type; 2077 typedef SymMatrixView<RT,I> realpart_type; 2078 typedef T& reference; 2079 2080 // 2081 // Constructors 2082 // 2083 2084 #define NEW_SIZE(s) itslen((s)*(s)), itsm(itslen), itss(s) \ 2085 TMV_DEFFIRSTLAST(itsm.get(),itsm.get()+itslen) 2086 2087 explicit inline SymMatrix(ptrdiff_t _size=0) : NEW_SIZE(_size) 2088 { 2089 TMVAssert(Attrib<A>::symmatrixok); 2090 TMVAssert(_size >= 0); 2091 #ifdef TMV_EXTRA_DEBUG 2092 setAllTo(T(888)); 2093 #endif 2094 } 2095 2096 inline SymMatrix(ptrdiff_t _size, const T& x) : NEW_SIZE(_size) 2097 { 2098 TMVAssert(Attrib<A>::symmatrixok); 2099 TMVAssert(_size >= 0); 2100 setAllTo(x); 2101 } 2102 2103 inline SymMatrix(const type& rhs) : NEW_SIZE(rhs.size()) 2104 { 2105 TMVAssert(Attrib<A>::symmatrixok); 2106 std::copy(rhs.cptr(),rhs.cptr()+itslen,itsm.get()); 2107 } 2108 2109 inline SymMatrix(const GenSymMatrix<RT>& rhs) : NEW_SIZE(rhs.size()) 2110 { 2111 TMVAssert(Attrib<A>::symmatrixok); 2112 if (rhs.issym()) 2113 rhs.assignToS(view()); 2114 else { 2115 if (uplo() == Upper) 2116 upperTri() = rhs.upperTri(); 2117 else 2118 lowerTri() = rhs.lowerTri(); 2119 } 2120 } 2121 2122 inline SymMatrix(const GenSymMatrix<CT>& rhs) : NEW_SIZE(rhs.size()) 2123 { 2124 TMVAssert(Attrib<A>::symmatrixok); 2125 TMVAssert(isComplex(T())); 2126 if (rhs.issym()) 2127 rhs.assignToS(view()); 2128 else { 2129 if (uplo() == Upper) 2130 upperTri() = rhs.upperTri(); 2131 else 2132 lowerTri() = rhs.lowerTri(); 2133 } 2134 } 2135 2136 template <typename T2> 2137 inline SymMatrix(const GenSymMatrix<T2>& rhs) : NEW_SIZE(rhs.size()) 2138 { 2139 TMVAssert(Attrib<A>::symmatrixok); 2140 if (uplo() == Upper) 2141 upperTri() = rhs.upperTri(); 2142 else 2143 lowerTri() = rhs.lowerTri(); 2144 } 2145 2146 template <typename T2> 2147 inline SymMatrix(const GenMatrix<T2>& rhs) : 2148 NEW_SIZE(rhs.rowsize()) 2149 { 2150 TMVAssert(Attrib<A>::symmatrixok); 2151 TMVAssert(rhs.isSquare()); 2152 if (uplo() == Upper) 2153 upperTri() = rhs.upperTri(); 2154 else 2155 lowerTri() = rhs.lowerTri(); 2156 } 2157 2158 inline SymMatrix(const AssignableToSymMatrix<RT>& m2) : 2159 NEW_SIZE(m2.size()) 2160 { 2161 TMVAssert(Attrib<A>::symmatrixok); 2162 TMVAssert(m2.issym()); 2163 m2.assignToS(view()); 2164 } 2165 2166 inline SymMatrix(const AssignableToSymMatrix<CT>& m2) : 2167 NEW_SIZE(m2.size()) 2168 { 2169 TMVAssert(Attrib<A>::symmatrixok); 2170 TMVAssert(isComplex(T())); 2171 TMVAssert(m2.issym()); 2172 m2.assignToS(view()); 2173 } 2174 2175 inline explicit SymMatrix(const GenDiagMatrix<RT>& m2) : 2176 NEW_SIZE(m2.size()) 2177 { 2178 TMVAssert(Attrib<A>::symmatrixok); 2179 TMVAssert(size() == m2.size()); 2180 setZero(); 2181 m2.assignToD(DiagMatrixViewOf(diag())); 2182 } 2183 2184 inline explicit SymMatrix(const GenDiagMatrix<CT>& m2) : 2185 NEW_SIZE(m2.size()) 2186 { 2187 TMVAssert(Attrib<A>::symmatrixok); 2188 TMVAssert(isComplex(T())); 2189 TMVAssert(size() == m2.size()); 2190 setZero(); 2191 m2.assignToD(DiagMatrixViewOf(diag())); 2192 } 2193 2194 template <typename T2> 2195 inline SymMatrix(const OProdVV<T,T2,T2>& opvv) : 2196 NEW_SIZE(opvv.colsize()) 2197 { 2198 TMVAssert(Attrib<A>::symmatrixok); 2199 TMVAssert(opvv.colsize() == opvv.rowsize()); 2200 TMVAssert(opvv.getV1().isSameAs(opvv.getV2().view())); 2201 Rank1Update<false>(opvv.getX(),opvv.getV1(),view()); 2202 } 2203 2204 template <typename T2> 2205 inline SymMatrix(const ProdMM<T,T2,T2>& pmm) : 2206 NEW_SIZE(pmm.colsize()) 2207 { 2208 TMVAssert(Attrib<A>::symmatrixok); 2209 TMVAssert(pmm.colsize() == pmm.rowsize()); 2210 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().transpose())); 2211 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 2212 } 2213 2214 template <typename T2> 2215 inline SymMatrix(const ProdUL<T,T2,T2>& pmm) : 2216 NEW_SIZE(pmm.colsize()) 2217 { 2218 TMVAssert(Attrib<A>::symmatrixok); 2219 TMVAssert(pmm.colsize() == pmm.rowsize()); 2220 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().transpose())); 2221 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 2222 } 2223 2224 template <typename T2> 2225 inline SymMatrix(const ProdLU<T,T2,T2>& pmm) : 2226 NEW_SIZE(pmm.colsize()) 2227 { 2228 TMVAssert(Attrib<A>::symmatrixok); 2229 TMVAssert(pmm.colsize() == pmm.rowsize()); 2230 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().transpose())); 2231 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 2232 } 2233 2234 #undef NEW_SIZE 2235 2236 virtual inline ~SymMatrix() 2237 { 2238 TMV_SETFIRSTLAST(0,0); 2239 #ifdef TMV_EXTRA_DEBUG 2240 setAllTo(T(999)); 2241 #endif 2242 } 2243 2244 // 2245 // Op= 2246 // 2247 2248 inline type& operator=(const type& m2) 2249 { 2250 TMVAssert(size() == m2.size()); 2251 if (&m2 != this) 2252 std::copy(m2.cptr(),m2.cptr()+itslen,itsm.get()); 2253 return *this; 2254 } 2255 2256 inline type& operator=(const GenSymMatrix<RT>& m2) 2257 { 2258 TMVAssert(size() == m2.size()); 2259 TMVAssert(m2.issym()); 2260 m2.assignToS(view()); 2261 return *this; 2262 } 2263 2264 inline type& operator=(const GenSymMatrix<CT>& m2) 2265 { 2266 TMVAssert(isComplex(T())); 2267 TMVAssert(size() == m2.size()); 2268 TMVAssert(m2.issym()); 2269 m2.assignToS(view()); 2270 return *this; 2271 } 2272 2273 template <typename T2> 2274 inline type& operator=(const GenSymMatrix<T2>& m2) 2275 { 2276 TMVAssert(size() == m2.size()); 2277 TMVAssert(isReal(T2()) || m2.issym()); 2278 upperTri() = m2.upperTri(); 2279 return *this; 2280 } 2281 2282 inline type& operator=(const T& x) 2283 { return setToIdentity(x); } 2284 2285 inline type& operator=(const AssignableToSymMatrix<RT>& m2) 2286 { 2287 TMVAssert(size() == m2.size()); 2288 TMVAssert(m2.issym()); 2289 m2.assignToS(view()); 2290 return *this; 2291 } 2292 2293 inline type& operator=(const AssignableToSymMatrix<CT>& m2) 2294 { 2295 TMVAssert(isComplex(T())); 2296 TMVAssert(size() == m2.size()); 2297 TMVAssert(m2.issym()); 2298 m2.assignToS(view()); 2299 return *this; 2300 } 2301 2302 inline type& operator=(const GenDiagMatrix<RT>& m2) 2303 { 2304 TMVAssert(size() == m2.size()); 2305 view() = m2; 2306 return *this; 2307 } 2308 2309 inline type& operator=(const GenDiagMatrix<CT>& m2) 2310 { 2311 TMVAssert(isComplex(T())); 2312 TMVAssert(size() == m2.size()); 2313 view() = m2; 2314 return *this; 2315 } 2316 2317 template <typename T2> 2318 inline type& operator=(const OProdVV<T,T2,T2>& opvv) 2319 { 2320 TMVAssert(size() == opvv.colsize()); 2321 TMVAssert(size() == opvv.rowsize()); 2322 TMVAssert(opvv.getV1().isSameAs(opvv.getV2().view())); 2323 Rank1Update<false>(opvv.getX(),opvv.getV1(),view()); 2324 return *this; 2325 } 2326 2327 template <typename T2> 2328 inline type& operator=(const ProdMM<T,T2,T2>& pmm) 2329 { 2330 TMVAssert(size() == pmm.colsize()); 2331 TMVAssert(size() == pmm.rowsize()); 2332 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().transpose())); 2333 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 2334 return *this; 2335 } 2336 2337 template <typename T2> 2338 inline type& operator=(const ProdUL<T,T2,T2>& pmm) 2339 { 2340 TMVAssert(size() == pmm.colsize()); 2341 TMVAssert(size() == pmm.rowsize()); 2342 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().transpose())); 2343 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 2344 return *this; 2345 } 2346 2347 template <typename T2> 2348 inline type& operator=(const ProdLU<T,T2,T2>& pmm) 2349 { 2350 TMVAssert(size() == pmm.colsize()); 2351 TMVAssert(size() == pmm.rowsize()); 2352 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().transpose())); 2353 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 2354 return *this; 2355 } 2356 2357 2358 // 2359 // Access 2360 // 2361 2362 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 2363 { 2364 if (I==int(CStyle)) { 2365 TMVAssert(i>=0 && i<size()); 2366 TMVAssert(j>=0 && j<size()); 2367 return cref(i,j); 2368 } else { 2369 TMVAssert(i>0 && i<=size()); 2370 TMVAssert(j>0 && j<=size()); 2371 return cref(i-1,j-1); 2372 } 2373 } 2374 2375 inline T& operator()(ptrdiff_t i, ptrdiff_t j) 2376 { 2377 if (I==int(CStyle)) { 2378 TMVAssert(i>=0 && i<size()); 2379 TMVAssert(j>=0 && j<size()); 2380 return ref(i,j); 2381 } else { 2382 TMVAssert(i>0 && i<=size()); 2383 TMVAssert(j>0 && j<=size()); 2384 return ref(i-1,j-1); 2385 } 2386 } 2387 2388 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 2389 { 2390 if (I==int(FortranStyle)) { 2391 TMVAssert(i>0 && i<=size()); 2392 --i; 2393 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()); 2394 --j1; 2395 } else { 2396 TMVAssert(i>=0 && i<size()); 2397 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()); 2398 } 2399 if ((uplo()==Upper && i-j1<=0) || (uplo()==Lower && j2-i<=1)) 2400 return const_vec_type( 2401 itsm.get()+i*stepi()+j1*stepj(), 2402 j2-j1,stepj(),NonConj); 2403 else 2404 return const_vec_type( 2405 itsm.get()+i*stepj()+j1*stepi(), 2406 j2-j1,stepi(),NonConj); 2407 } 2408 2409 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 2410 { 2411 if (I==int(FortranStyle)) { 2412 TMVAssert(j>0 && j<=size()); 2413 --j; 2414 TMVAssert(i1>0 && i1-i2<=0 && i2<=size()); 2415 --i1; 2416 } else { 2417 TMVAssert(j>=0 && j<size()); 2418 TMVAssert(i1>=0 && i1-i2<=0 && i2<=size()); 2419 } 2420 if ((uplo()==Upper && i2-j<=1) || (uplo()==Lower && j-i1<=0)) 2421 return const_vec_type( 2422 itsm.get()+i1*stepi()+j*stepj(), 2423 i2-i1,stepi(),NonConj); 2424 else 2425 return const_vec_type( 2426 itsm.get()+i1*stepj()+j*stepi(), 2427 i2-i1,stepj(),NonConj); 2428 } 2429 2430 inline const_vec_type diag() const 2431 { return const_vec_type( 2432 itsm.get(),size(),stepi()+stepj(),NonConj); } 2433 2434 inline const_vec_type diag(ptrdiff_t i) const 2435 { 2436 TMVAssert(i>=-size() && i<=size()); 2437 i = std::abs(i); 2438 return const_vec_type( 2439 itsm.get()+i*(uplo()==Upper?stepj():stepi()), 2440 size()-i,stepi()+stepj(),NonConj); 2441 } 2442 2443 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 2444 { 2445 TMVAssert(i>=-size() && i<=size()); 2446 if (I==int(FortranStyle)) { 2447 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()-std::abs(i)); 2448 --j1; 2449 } else { 2450 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()-std::abs(i)); 2451 } 2452 i = std::abs(i); 2453 const ptrdiff_t ds = stepi()+stepj(); 2454 return const_vec_type( 2455 itsm.get()+i*(uplo()==Upper?stepj():stepi())+j1*ds, 2456 j2-j1,ds,NonConj); 2457 } 2458 2459 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 2460 { 2461 if (I==int(FortranStyle)) { 2462 TMVAssert(i>0 && i<=size()); 2463 --i; 2464 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()); 2465 --j1; 2466 } else { 2467 TMVAssert(i>=0 && i<size()); 2468 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()); 2469 } 2470 if ((uplo()==Upper && i-j1<=0) || (uplo()==Lower && j2-i<=1)) 2471 return vec_type( 2472 itsm.get()+i*stepi()+j1*stepj(), 2473 j2-j1,stepj(),NonConj TMV_FIRSTLAST); 2474 else 2475 return vec_type( 2476 itsm.get()+i*stepj()+j1*stepi(), 2477 j2-j1,stepi(),NonConj TMV_FIRSTLAST); 2478 } 2479 2480 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 2481 { 2482 if (I==int(FortranStyle)) { 2483 TMVAssert(j>0 && j<=size()); 2484 --j; 2485 TMVAssert(i1>0 && i1-i2<=0 && i2<=size()); 2486 --i1; 2487 } else { 2488 TMVAssert(j>=0 && j<size()); 2489 TMVAssert(i1>=0 && i1-i2<=0 && i2<=size()); 2490 } 2491 if ((uplo()==Upper && i2-j<=1) || (uplo()==Lower && j-i1<=0)) 2492 return vec_type( 2493 itsm.get()+i1*stepi()+j*stepj(), 2494 i2-i1,stepi(),NonConj TMV_FIRSTLAST); 2495 else 2496 return vec_type( 2497 itsm.get()+i1*stepj()+j*stepi(), 2498 i2-i1,stepj(),NonConj TMV_FIRSTLAST); 2499 } 2500 2501 inline vec_type diag() 2502 { 2503 return vec_type( 2504 itsm.get(),size(),stepi()+stepj(),NonConj 2505 TMV_FIRSTLAST); 2506 } 2507 2508 inline vec_type diag(ptrdiff_t i) 2509 { 2510 TMVAssert(i>=-size() && i<=size()); 2511 i = std::abs(i); 2512 TMVAssert(i<=size()); 2513 return vec_type( 2514 itsm.get()+i*(uplo()==Upper?stepj():stepi()), 2515 size()-i,stepi()+stepj(),NonConj TMV_FIRSTLAST); 2516 } 2517 2518 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 2519 { 2520 TMVAssert(i>=-size() && i<=size()); 2521 if (I==int(FortranStyle)) { 2522 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()-std::abs(i)); 2523 --j1; 2524 } else { 2525 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()-std::abs(i)); 2526 } 2527 i = std::abs(i); 2528 const ptrdiff_t ds = stepi()+stepj(); 2529 return vec_type( 2530 itsm.get()+i*(uplo()==Upper?stepj():stepi())+j1*ds, 2531 j2-j1,ds,NonConj TMV_FIRSTLAST); 2532 } 2533 2534 2535 // 2536 // Modifying Functions 2537 // 2538 2539 inline type& setZero() 2540 { std::fill_n(itsm.get(),itslen,T(0)); return *this; } 2541 2542 inline type& setAllTo(const T& x) 2543 { upperTri().setAllTo(x); return *this; } 2544 2545 inline type& addToAll(const T& x) 2546 { upperTri().addToAll(x); return *this; } 2547 2548 inline type& clip(RT thresh) 2549 { upperTri().clip(thresh); return *this; } 2550 2551 inline type& conjugateSelf() 2552 { if (isComplex(T())) upperTri().conjugateSelf(); return *this; } 2553 2554 inline type& transposeSelf() 2555 { return *this; } 2556 2557 inline type& setToIdentity(const T& x=T(1)) 2558 { setZero(); diag().setAllTo(x); return *this; } 2559 2560 inline type& swapRowsCols(ptrdiff_t i1, ptrdiff_t i2) 2561 { view().swapRowsCols(i1,i2); return *this; } 2562 2563 inline type& permuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 2564 { view().permuteRowsCols(p,i1,i2); return *this; } 2565 2566 inline type& reversePermuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 2567 { view().reversePermuteRowsCols(p,i1,i2); return *this; } 2568 2569 inline type& permuteRowsCols(const ptrdiff_t* p) 2570 { view().permuteRowsCols(p); return *this; } 2571 2572 inline type& reversePermuteRowsCols(const ptrdiff_t* p) 2573 { view().reversePermuteRowsCols(p); return *this; } 2574 2575 // 2576 // subMatrix 2577 // 2578 2579 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 2580 { 2581 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 2582 if (I==int(FortranStyle)) { --i1; --j1; } 2583 if ( (uplo()==Upper && i2-j1<=1) || (uplo()==Lower && j2-i1<=1) ) 2584 return const_rec_type( 2585 itsm.get()+i1*stepi()+j1*stepj(), 2586 i2-i1,j2-j1,stepi(),stepj(),NonConj); 2587 else 2588 return const_rec_type( 2589 itsm.get()+i1*stepj()+j1*stepi(), 2590 i2-i1,j2-j1,stepj(),stepi(),NonConj); 2591 } 2592 2593 inline const_rec_type subMatrix( 2594 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 2595 { 2596 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 2597 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 2598 if ( (uplo()==Upper && i2-j1<=istep) || 2599 (uplo()==Lower && j2-i1<=jstep) ) 2600 return const_rec_type( 2601 itsm.get()+i1*stepi()+j1*stepj(), 2602 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(), 2603 NonConj); 2604 else 2605 return const_rec_type( 2606 itsm.get()+i1*stepj()+j1*stepi(), 2607 (i2-i1)/istep,(j2-j1)/jstep,istep*stepj(),jstep*stepi(), 2608 NonConj); 2609 } 2610 2611 inline const_vec_type subVector( 2612 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 2613 { 2614 TMVAssert(view().hasSubVector(i,j,istep,jstep,n)); 2615 if (I==int(FortranStyle)) { --i; --j; } 2616 if ((uplo()==Upper && i-j<=0) || (uplo()==Lower && j-i<=0)) 2617 return const_vec_type( 2618 itsm.get()+i*stepi()+j*stepj(),n, 2619 istep*stepi()+jstep*stepj(),NonConj); 2620 else 2621 return const_vec_type( 2622 itsm.get()+i*stepj()+j*stepi(),n, 2623 istep*stepj()+jstep*stepi(),NonConj); 2624 } 2625 2626 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) const 2627 { 2628 TMVAssert(view().hasSubSymMatrix(i1,i2,1)); 2629 if (I==int(FortranStyle)) { --i1; } 2630 return const_view_type( 2631 itsm.get()+i1*(stepi()+stepj()),i2-i1, 2632 stepi(),stepj(),Sym,uplo(),NonConj); 2633 } 2634 2635 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 2636 { 2637 TMVAssert(view().hasSubSymMatrix(i1,i2,istep)); 2638 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 2639 return const_view_type( 2640 itsm.get()+i1*(stepi()+stepj()), 2641 (i2-i1)/istep,istep*stepi(),istep*stepj(),Sym,uplo(), 2642 NonConj); 2643 } 2644 2645 inline const_uppertri_type upperTri(DiagType dt = NonUnitDiag) const 2646 { 2647 return uplo()==Upper ? 2648 const_uppertri_type( 2649 itsm.get(),size(),stepi(),stepj(),dt,NonConj) : 2650 const_uppertri_type( 2651 itsm.get(),size(),stepj(),stepi(),dt,NonConj); 2652 } 2653 2654 inline const_uppertri_type unitUpperTri() const 2655 { 2656 return uplo()==Upper ? 2657 const_uppertri_type( 2658 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj) : 2659 const_uppertri_type( 2660 itsm.get(),size(),stepj(),stepi(),UnitDiag,NonConj); 2661 } 2662 2663 inline const_lowertri_type lowerTri(DiagType dt = NonUnitDiag) const 2664 { 2665 return uplo()==Lower ? 2666 const_lowertri_type( 2667 itsm.get(),size(),stepi(),stepj(),dt,NonConj) : 2668 const_lowertri_type( 2669 itsm.get(),size(),stepj(),stepi(),dt,NonConj); 2670 } 2671 2672 inline const_lowertri_type unitLowerTri() const 2673 { 2674 return uplo()==Lower ? 2675 const_lowertri_type( 2676 itsm.get(),size(), 2677 stepi(),stepj(),UnitDiag,NonConj) : 2678 const_lowertri_type( 2679 itsm.get(),size(),stepj(),stepi(),UnitDiag,NonConj); 2680 } 2681 2682 inline const_realpart_type realPart() const 2683 { 2684 return const_realpart_type( 2685 reinterpret_cast<const RT*>(itsm.get()),size(), 2686 isReal(T()) ? stepi() : 2*stepi(), 2687 isReal(T()) ? stepj() : 2*stepj(), Sym,uplo(),NonConj); 2688 } 2689 2690 inline const_realpart_type imagPart() const 2691 { 2692 TMVAssert(isComplex(T())); 2693 return const_realpart_type( 2694 reinterpret_cast<const RT*>(itsm.get())+1,size(), 2695 2*stepi(),2*stepj(),Sym,uplo(),NonConj); 2696 } 2697 2698 inline const_view_type view() const 2699 { 2700 return const_view_type( 2701 itsm.get(),size(),stepi(),stepj(),Sym,uplo(),NonConj); 2702 } 2703 2704 inline const_view_type transpose() const 2705 { 2706 return const_view_type( 2707 itsm.get(),size(),stepj(),stepi(),Sym,TMV_UTransOf(U),NonConj); 2708 } 2709 2710 inline const_view_type conjugate() const 2711 { 2712 return const_view_type( 2713 itsm.get(),size(),stepi(),stepj(),Sym,uplo(), 2714 TMV_ConjOf(T,NonConj)); 2715 } 2716 2717 inline const_view_type adjoint() const 2718 { 2719 return const_view_type( 2720 itsm.get(),size(),stepj(),stepi(), 2721 Sym,TMV_UTransOf(U),TMV_ConjOf(T,NonConj)); 2722 } 2723 2724 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 2725 { 2726 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 2727 if (I==int(FortranStyle)) { --i1; --j1; } 2728 if ( (uplo()==Upper && i2-j1<=1) || (uplo()==Lower && j2-i1<=1) ) 2729 return rec_type( 2730 itsm.get()+i1*stepi()+j1*stepj(), 2731 i2-i1,j2-j1,stepi(),stepj(),NonConj TMV_FIRSTLAST); 2732 else 2733 return rec_type( 2734 itsm.get()+i1*stepj()+j1*stepi(), 2735 i2-i1,j2-j1,stepj(),stepi(),NonConj TMV_FIRSTLAST); 2736 } 2737 2738 inline rec_type subMatrix( 2739 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 2740 { 2741 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 2742 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 2743 if ( (uplo()==Upper && i2-j1<=istep) || 2744 (uplo()==Lower && j2-i1<=jstep) ) 2745 return rec_type( 2746 itsm.get()+i1*stepi()+j1*stepj(), 2747 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(), 2748 NonConj TMV_FIRSTLAST); 2749 else 2750 return rec_type( 2751 itsm.get()+i1*stepj()+j1*stepi(), 2752 (i2-i1)/istep,(j2-j1)/jstep,istep*stepj(),jstep*stepi(), 2753 NonConj TMV_FIRSTLAST); 2754 } 2755 2756 inline vec_type subVector( 2757 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) 2758 { 2759 TMVAssert(view().hasSubVector(i,j,istep,jstep,n)); 2760 if (I==int(FortranStyle)) { --i; --j; } 2761 if ((uplo()==Upper && i-j<=0) || (uplo()==Lower && j-i<=0)) 2762 return vec_type( 2763 itsm.get()+i*stepi()+j*stepj(),n, 2764 istep*stepi()+jstep*stepj(),NonConj TMV_FIRSTLAST); 2765 else 2766 return vec_type( 2767 itsm.get()+i*stepj()+j*stepi(),n, 2768 istep*stepj()+jstep*stepi(),NonConj TMV_FIRSTLAST); 2769 } 2770 2771 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) 2772 { 2773 TMVAssert(view().hasSubSymMatrix(i1,i2,1)); 2774 if (I==int(FortranStyle)) { --i1; } 2775 return view_type( 2776 itsm.get()+i1*(stepi()+stepj()),i2-i1, 2777 stepi(),stepj(),Sym,uplo(),NonConj TMV_FIRSTLAST); 2778 } 2779 2780 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 2781 { 2782 TMVAssert(view().hasSubSymMatrix(i1,i2,istep)); 2783 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 2784 return view_type( 2785 itsm.get()+i1*(stepi()+stepj()),(i2-i1)/istep, 2786 istep*stepi(),istep*stepj(),Sym,uplo(), NonConj TMV_FIRSTLAST); 2787 } 2788 2789 inline uppertri_type upperTri(DiagType dt = NonUnitDiag) 2790 { 2791 return uplo()==Upper ? 2792 uppertri_type( 2793 itsm.get(),size(),stepi(),stepj(),dt,NonConj 2794 TMV_FIRSTLAST) : 2795 uppertri_type( 2796 itsm.get(),size(),stepj(),stepi(),dt,NonConj 2797 TMV_FIRSTLAST); 2798 } 2799 2800 inline uppertri_type unitUpperTri() 2801 { 2802 return uplo()==Upper ? 2803 uppertri_type( 2804 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj 2805 TMV_FIRSTLAST) : 2806 uppertri_type( 2807 itsm.get(),size(),stepj(),stepi(),UnitDiag,NonConj 2808 TMV_FIRSTLAST); 2809 } 2810 2811 inline lowertri_type lowerTri(DiagType dt = NonUnitDiag) 2812 { 2813 return uplo()==Lower ? 2814 lowertri_type( 2815 itsm.get(),size(),stepi(),stepj(),dt,NonConj 2816 TMV_FIRSTLAST) : 2817 lowertri_type( 2818 itsm.get(),size(),stepj(),stepi(),dt,NonConj 2819 TMV_FIRSTLAST); 2820 } 2821 2822 inline lowertri_type unitLowerTri() 2823 { 2824 return uplo()==Lower ? 2825 lowertri_type( 2826 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj 2827 TMV_FIRSTLAST) : 2828 lowertri_type( 2829 itsm.get(),size(),stepj(),stepi(),UnitDiag,NonConj 2830 TMV_FIRSTLAST); 2831 } 2832 2833 inline realpart_type realPart() 2834 { 2835 return realpart_type( 2836 reinterpret_cast<RT*>(itsm.get()),size(), 2837 isReal(T()) ? stepi() : 2*stepi(), 2838 isReal(T()) ? stepj() : 2*stepj(), Sym,uplo(),NonConj 2839 #ifdef TMVFLDEBUG 2840 ,reinterpret_cast<const RT*>(_first) 2841 ,reinterpret_cast<const RT*>(_last) 2842 #endif 2843 ); 2844 } 2845 2846 inline realpart_type imagPart() 2847 { 2848 TMVAssert(isComplex(T())); 2849 return realpart_type( 2850 reinterpret_cast<RT*>(itsm.get())+1,size(), 2851 2*stepi(),2*stepj(),Sym,uplo(),NonConj 2852 #ifdef TMVFLDEBUG 2853 ,reinterpret_cast<const RT*>(_first)+1 2854 ,reinterpret_cast<const RT*>(_last)+1 2855 #endif 2856 ); 2857 } 2858 2859 inline view_type view() 2860 { 2861 return view_type( 2862 itsm.get(),size(),stepi(),stepj(),Sym,uplo(),NonConj 2863 TMV_FIRSTLAST); 2864 } 2865 2866 inline view_type transpose() 2867 { 2868 return view_type( 2869 itsm.get(),size(),stepj(),stepi(), 2870 Sym,TMV_UTransOf(U),NonConj TMV_FIRSTLAST); 2871 } 2872 2873 inline view_type conjugate() 2874 { 2875 return view_type( 2876 itsm.get(),size(),stepi(),stepj(), 2877 Sym,uplo(),TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 2878 } 2879 2880 inline view_type adjoint() 2881 { 2882 return view_type( 2883 itsm.get(),size(),stepj(),stepi(), 2884 Sym,TMV_UTransOf(U),TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 2885 } 2886 2887 // 2888 // I/O 2889 // 2890 2891 void read(const TMV_Reader& reader); 2892 2893 inline ptrdiff_t size() const { return itss; } 2894 inline const T* cptr() const { return itsm.get(); } 2895 inline T* ptr() { return itsm.get(); } 2896 inline ptrdiff_t stepi() const { return S==int(RowMajor) ? itss : 1; } 2897 inline ptrdiff_t stepj() const { return S==int(RowMajor) ? 1 : itss; } 2898 inline SymType sym() const { return Sym; } 2899 inline UpLoType uplo() const { return static_cast<UpLoType>(U); } 2900 inline ConjType ct() const { return NonConj; } 2901 inline bool isrm() const { return S==int(RowMajor); } 2902 inline bool iscm() const { return S==int(ColMajor); } 2903 inline bool isconj() const { return false; } 2904 inline bool isherm() const { return isReal(T()); } 2905 inline bool issym() const { return true; } 2906 inline bool isupper() const { return U==int(Upper); } 2907 2908 inline T& ref(ptrdiff_t i, ptrdiff_t j) 2909 { 2910 if ((uplo()==Upper && i <= j) || (uplo()==Lower && i>=j)) 2911 return itsm.get()[i*stepi() + j*stepj()]; 2912 else 2913 return itsm.get()[j*stepi() + i*stepj()]; 2914 } 2915 2916 inline T cref(ptrdiff_t i, ptrdiff_t j) const 2917 { 2918 if ((uplo()==Upper && i <= j) || (uplo()==Lower && i>=j)) 2919 return itsm.get()[i*stepi() + j*stepj()]; 2920 else 2921 return itsm.get()[j*stepi() + i*stepj()]; 2922 } 2923 2924 inline void resize(ptrdiff_t s) 2925 { 2926 TMVAssert(s >= 0); 2927 itslen = s*s; 2928 itsm.resize(itslen); 2929 itss = s; 2930 DivHelper<T>::resetDivType(); 2931 #ifdef TMVFLDEBUG 2932 _first = itsm.get(); 2933 _last = _first+itslen; 2934 #endif 2935 #ifdef TMV_EXTRA_DEBUG 2936 setAllTo(T(888)); 2937 #endif 2938 } 2939 2940 protected : 2941 2942 ptrdiff_t itslen; 2943 AlignedArray<T> itsm; 2944 ptrdiff_t itss; 2945 2946 #ifdef TMVFLDEBUG 2947 public: 2948 const T* _first; 2949 const T* _last; 2950 #endif 2951 2952 friend void Swap(SymMatrix<T,A>& m1, SymMatrix<T,A>& m2) 2953 { 2954 TMVAssert(m1.size() == m2.size()); 2955 m1.itsm.swapWith(m2.itsm); 2956 #ifdef TMVFLDEBUG 2957 TMV_SWAP(m1._first,m2._first); 2958 TMV_SWAP(m1._last,m2._last); 2959 #endif 2960 } 2961 2962 }; // SymMatrix 2963 2964 template <typename T, int A> 2965 class HermMatrix : public GenSymMatrix<T> 2966 { 2967 public: 2968 2969 enum { S = A & AllStorageType }; 2970 enum { U = A & Upper }; 2971 enum { I = A & FortranStyle }; 2972 2973 typedef TMV_RealType(T) RT; 2974 typedef TMV_ComplexType(T) CT; 2975 typedef HermMatrix<T,A> type; 2976 typedef type copy_type; 2977 typedef GenSymMatrix<T> base; 2978 typedef ConstSymMatrixView<T,I> const_view_type; 2979 typedef const_view_type const_transpose_type; 2980 typedef const_view_type const_conjugate_type; 2981 typedef const_view_type const_adjoint_type; 2982 typedef ConstVectorView<T,I> const_vec_type; 2983 typedef ConstMatrixView<T,I> const_rec_type; 2984 typedef ConstUpperTriMatrixView<T,I> const_uppertri_type; 2985 typedef ConstLowerTriMatrixView<T,I> const_lowertri_type; 2986 typedef ConstSymMatrixView<RT,I> const_realpart_type; 2987 typedef SymMatrixView<T,I> view_type; 2988 typedef view_type transpose_type; 2989 typedef view_type conjugate_type; 2990 typedef view_type adjoint_type; 2991 typedef VectorView<T,I> vec_type; 2992 typedef MatrixView<T,I> rec_type; 2993 typedef UpperTriMatrixView<T,I> uppertri_type; 2994 typedef LowerTriMatrixView<T,I> lowertri_type; 2995 typedef SymMatrixView<RT,I> realpart_type; 2996 typedef typename RefHelper<T>::reference reference; 2997 2998 // 2999 // Constructors 3000 // 3001 3002 #define NEW_SIZE(s) itslen((s)*(s)), itsm(itslen), itss(s) \ 3003 TMV_DEFFIRSTLAST(itsm.get(),itsm.get()+itslen) 3004 3005 explicit inline HermMatrix(ptrdiff_t _size=0) : NEW_SIZE(_size) 3006 { 3007 TMVAssert(Attrib<A>::symmatrixok); 3008 TMVAssert(_size >= 0); 3009 #ifdef TMV_EXTRA_DEBUG 3010 setAllTo(T(888)); 3011 #else 3012 if (isComplex(T())) diag().imagPart().setZero(); 3013 #endif 3014 } 3015 3016 HermMatrix(ptrdiff_t _size, const RT& x) : NEW_SIZE(_size) 3017 { 3018 TMVAssert(Attrib<A>::symmatrixok); 3019 TMVAssert(_size >= 0); 3020 setAllTo(x); 3021 } 3022 3023 HermMatrix(const type& rhs) : NEW_SIZE(rhs.size()) 3024 { 3025 TMVAssert(Attrib<A>::symmatrixok); 3026 std::copy(rhs.cptr(),rhs.cptr()+itslen,itsm.get()); 3027 } 3028 3029 HermMatrix(const GenSymMatrix<RT>& rhs) : NEW_SIZE(rhs.size()) 3030 { 3031 TMVAssert(Attrib<A>::symmatrixok); 3032 if (rhs.isherm()) rhs.assignToS(view()); 3033 else if (uplo() == Upper) upperTri() = rhs.upperTri(); 3034 else lowerTri() = rhs.lowerTri(); 3035 } 3036 3037 HermMatrix(const GenSymMatrix<CT>& rhs) : NEW_SIZE(rhs.size()) 3038 { 3039 TMVAssert(Attrib<A>::symmatrixok); 3040 TMVAssert(isComplex(T())); 3041 if (rhs.isherm()) 3042 rhs.assignToS(view()); 3043 else { 3044 if (uplo() == Upper) upperTri() = rhs.upperTri(); 3045 else lowerTri() = rhs.lowerTri(); 3046 diag().imagPart().setZero(); 3047 } 3048 } 3049 3050 template <typename T2> 3051 inline HermMatrix(const GenSymMatrix<T2>& rhs) : NEW_SIZE(rhs.size()) 3052 { 3053 TMVAssert(Attrib<A>::symmatrixok); 3054 if (uplo() == Upper) upperTri() = rhs.upperTri(); 3055 else lowerTri() = rhs.lowerTri(); 3056 if (isComplex(T()) && isComplex(T2()) && rhs.issym()) 3057 diag().imagPart().setZero(); 3058 } 3059 3060 template <typename T2> 3061 inline HermMatrix(const GenMatrix<T2>& rhs) : NEW_SIZE(rhs.rowsize()) 3062 { 3063 TMVAssert(Attrib<A>::symmatrixok); 3064 TMVAssert(rhs.isSquare()); 3065 if (uplo() == Upper) 3066 upperTri() = rhs.upperTri(); 3067 else 3068 lowerTri() = rhs.lowerTri(); 3069 if (isComplex(T()) && isComplex(T2())) diag().imagPart().setZero(); 3070 } 3071 3072 inline HermMatrix(const AssignableToSymMatrix<RT>& m2) : 3073 NEW_SIZE(m2.size()) 3074 { 3075 TMVAssert(Attrib<A>::symmatrixok); 3076 TMVAssert(m2.isherm()); 3077 m2.assignToS(view()); 3078 } 3079 3080 inline HermMatrix(const AssignableToSymMatrix<CT>& m2) : 3081 NEW_SIZE(m2.size()) 3082 { 3083 TMVAssert(Attrib<A>::symmatrixok); 3084 TMVAssert(isComplex(T())); 3085 TMVAssert(m2.isherm()); 3086 m2.assignToS(view()); 3087 } 3088 3089 inline explicit HermMatrix(const GenDiagMatrix<RT>& m2) : 3090 NEW_SIZE(m2.size()) 3091 { 3092 TMVAssert(Attrib<A>::symmatrixok); 3093 m2.assignToD(DiagMatrixViewOf(diag())); 3094 upperTri().offDiag().setZero(); 3095 } 3096 3097 inline explicit HermMatrix(const GenDiagMatrix<CT>& m2) : 3098 NEW_SIZE(m2.size()) 3099 { 3100 TMVAssert(Attrib<A>::symmatrixok); 3101 m2.assignToD(DiagMatrixViewOf(diag().realPart())); 3102 upperTri().offDiag().setZero(); 3103 if (isComplex(T())) diag().imagPart().setZero(); 3104 } 3105 3106 template <typename T2> 3107 inline HermMatrix(const OProdVV<T,T2,T2>& opvv) : 3108 NEW_SIZE(opvv.colsize()) 3109 { 3110 TMVAssert(Attrib<A>::symmatrixok); 3111 TMVAssert(opvv.colsize() == opvv.rowsize()); 3112 TMVAssert(opvv.getV1().isSameAs(opvv.getV2().conjugate())); 3113 TMVAssert(TMV_IMAG(opvv.getX()) == RT(0)); 3114 Rank1Update<false>(opvv.getX(),opvv.getV1(),view()); 3115 } 3116 3117 template <typename T2> 3118 inline HermMatrix(const ProdMM<T,T2,T2>& pmm) : 3119 NEW_SIZE(pmm.colsize()) 3120 { 3121 TMVAssert(Attrib<A>::symmatrixok); 3122 TMVAssert(pmm.colsize() == pmm.rowsize()); 3123 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().adjoint())); 3124 TMVAssert(TMV_IMAG(pmm.getX()) == RT(0)); 3125 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 3126 } 3127 3128 template <typename T2> 3129 inline HermMatrix(const ProdUL<T,T2,T2>& pmm) : 3130 NEW_SIZE(pmm.colsize()) 3131 { 3132 TMVAssert(Attrib<A>::symmatrixok); 3133 TMVAssert(pmm.colsize() == pmm.rowsize()); 3134 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().adjoint())); 3135 TMVAssert(TMV_IMAG(pmm.getX()) == RT(0)); 3136 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 3137 } 3138 3139 template <typename T2> 3140 inline HermMatrix(const ProdLU<T,T2,T2>& pmm) : 3141 NEW_SIZE(pmm.colsize()) 3142 { 3143 TMVAssert(Attrib<A>::symmatrixok); 3144 TMVAssert(pmm.colsize() == pmm.rowsize()); 3145 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().adjoint())); 3146 TMVAssert(TMV_IMAG(pmm.getX()) == RT(0)); 3147 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 3148 } 3149 3150 #undef NEW_SIZE 3151 3152 virtual inline ~HermMatrix() 3153 { 3154 TMV_SETFIRSTLAST(0,0); 3155 #ifdef TMV_EXTRA_DEBUG 3156 setAllTo(T(999)); 3157 #endif 3158 } 3159 3160 // 3161 // Op= 3162 // 3163 3164 inline type& operator=(const type& m2) 3165 { 3166 TMVAssert(m2.size() == size()); 3167 if (&m2 != this) 3168 std::copy(m2.cptr(),m2.cptr()+itslen,itsm.get()); 3169 return *this; 3170 } 3171 3172 inline type& operator=(const GenSymMatrix<RT>& m2) 3173 { 3174 TMVAssert(m2.size() == size()); 3175 TMVAssert(m2.isherm()); 3176 m2.assignToS(view()); 3177 return *this; 3178 } 3179 3180 inline type& operator=(const GenSymMatrix<CT>& m2) 3181 { 3182 TMVAssert(isComplex(T())); 3183 TMVAssert(m2.size() == size()); 3184 TMVAssert(m2.isherm()); 3185 m2.assignToS(view()); 3186 return *this; 3187 } 3188 3189 template <typename T2> 3190 inline type& operator=(const GenSymMatrix<T2>& m2) 3191 { 3192 TMVAssert(m2.size() == size()); 3193 TMVAssert(m2.isherm()); 3194 upperTri() = m2.upperTri(); 3195 return *this; 3196 } 3197 3198 inline type& operator=(const T& x) 3199 { 3200 TMVAssert(TMV_IMAG(x) == RT(0)); 3201 return setToIdentity(x); 3202 } 3203 3204 inline type& operator=(const AssignableToSymMatrix<RT>& m2) 3205 { 3206 TMVAssert(size() == m2.colsize()); 3207 TMVAssert(m2.isherm()); 3208 m2.assignToS(view()); 3209 return *this; 3210 } 3211 3212 inline type& operator=(const AssignableToSymMatrix<CT>& m2) 3213 { 3214 TMVAssert(isComplex(T())); 3215 TMVAssert(size() == m2.colsize()); 3216 TMVAssert(m2.isherm()); 3217 m2.assignToS(view()); 3218 return *this; 3219 } 3220 3221 inline type& operator=(const GenDiagMatrix<RT>& m2) 3222 { 3223 TMVAssert(size() == m2.colsize()); 3224 view() = m2; 3225 return *this; 3226 } 3227 3228 inline type& operator=(const GenDiagMatrix<CT>& m2) 3229 { 3230 TMVAssert(isComplex(T())); 3231 TMVAssert(size() == m2.colsize()); 3232 view() = m2; 3233 TMVAssert(this->isHermOK()); 3234 return *this; 3235 } 3236 3237 template <typename T2> 3238 inline type& operator=(const OProdVV<T,T2,T2>& opvv) 3239 { 3240 TMVAssert(size() == opvv.colsize()); 3241 TMVAssert(size() == opvv.rowsize()); 3242 TMVAssert(opvv.getV1().isSameAs(opvv.getV2().conjugate())); 3243 Rank1Update<false>(opvv.getX(),opvv.getV1(),view()); 3244 return *this; 3245 } 3246 3247 template <typename T2> 3248 inline type& operator=(const ProdMM<T,T2,T2>& pmm) 3249 { 3250 TMVAssert(size() == pmm.colsize()); 3251 TMVAssert(size() == pmm.rowsize()); 3252 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().adjoint())); 3253 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 3254 return *this; 3255 } 3256 3257 template <typename T2> 3258 inline type& operator=(const ProdUL<T,T2,T2>& pmm) 3259 { 3260 TMVAssert(size() == pmm.colsize()); 3261 TMVAssert(size() == pmm.rowsize()); 3262 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().adjoint())); 3263 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 3264 return *this; 3265 } 3266 3267 template <typename T2> 3268 inline type& operator=(const ProdLU<T,T2,T2>& pmm) 3269 { 3270 TMVAssert(size() == pmm.colsize()); 3271 TMVAssert(size() == pmm.rowsize()); 3272 TMVAssert(pmm.getM1().isSameAs(pmm.getM2().adjoint())); 3273 RankKUpdate<false>(pmm.getX(),pmm.getM1(),view()); 3274 return *this; 3275 } 3276 3277 3278 // 3279 // Access 3280 // 3281 3282 inline T operator()(ptrdiff_t i, ptrdiff_t j) const 3283 { 3284 if (I==int(CStyle)) { 3285 TMVAssert(i>=0 && i<size()); 3286 TMVAssert(j>=0 && j<size()); 3287 return cref(i,j); 3288 } else { 3289 TMVAssert(i>0 && i<=size()); 3290 TMVAssert(j>0 && j<=size()); 3291 return cref(i-1,j-1); 3292 } 3293 } 3294 3295 inline reference operator()(ptrdiff_t i, ptrdiff_t j) 3296 { 3297 if (I==int(CStyle)) { 3298 TMVAssert(i>=0 && i<size()); 3299 TMVAssert(j>=0 && j<size()); 3300 return ref(i,j); 3301 } else { 3302 TMVAssert(i>0 && i<=size()); 3303 TMVAssert(j>0 && j<=size()); 3304 return ref(i-1,j-1); 3305 } 3306 } 3307 3308 inline const_vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3309 { 3310 if (I==int(FortranStyle)) { 3311 TMVAssert(i>0 && i<=size()); 3312 --i; 3313 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()); 3314 --j1; 3315 } else { 3316 TMVAssert(i>=0 && i<size()); 3317 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()); 3318 } 3319 if ((uplo()==Upper && i-j1<=0) || (uplo()==Lower && j2-i<=1)) 3320 return const_vec_type( 3321 itsm.get()+i*stepi()+j1*stepj(), 3322 j2-j1,stepj(),NonConj); 3323 else 3324 return const_vec_type( 3325 itsm.get()+i*stepj()+j1*stepi(), 3326 j2-j1,stepi(),TMV_ConjOf(T,NonConj)); 3327 } 3328 3329 inline const_vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) const 3330 { 3331 if (I==int(FortranStyle)) { 3332 TMVAssert(j>0 && j<=size()); 3333 --j; 3334 TMVAssert(i1>0 && i1-i2<=0 && i2<=size()); 3335 --i1; 3336 } else { 3337 TMVAssert(j>=0 && j<size()); 3338 TMVAssert(i1>=0 && i1-i2<=0 && i2<=size()); 3339 } 3340 if ((uplo()==Upper && i2-j<=1) || (uplo()==Lower && j-i1<=0)) 3341 return const_vec_type( 3342 itsm.get()+i1*stepi()+j*stepj(), 3343 i2-i1,stepi(),NonConj); 3344 else 3345 return const_vec_type( 3346 itsm.get()+i1*stepj()+j*stepi(), 3347 i2-i1,stepj(),TMV_ConjOf(T,NonConj)); 3348 } 3349 3350 inline const_vec_type diag() const 3351 { return const_vec_type( 3352 itsm.get(),size(),stepi()+stepj(),NonConj); } 3353 3354 inline const_vec_type diag(ptrdiff_t i) const 3355 { 3356 TMVAssert(i>=-size() && i<=size()); 3357 ConjType newct = 3358 ((i>0) == (uplo()==Upper)) ? NonConj : TMV_ConjOf(T,NonConj); 3359 i = std::abs(i); 3360 TMVAssert(i<=size()); 3361 return const_vec_type( 3362 itsm.get()+i*(uplo()==Upper?stepj():stepi()), 3363 size()-i,stepi()+stepj(),newct); 3364 } 3365 3366 inline const_vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) const 3367 { 3368 TMVAssert(i>=-size() && i<=size()); 3369 if (I==int(FortranStyle)) { 3370 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()-std::abs(i)); 3371 --j1; 3372 } else { 3373 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()-std::abs(i)); 3374 } 3375 ConjType newct = 3376 ((i>0) == (uplo()==Upper)) ? NonConj : TMV_ConjOf(T,NonConj); 3377 i = std::abs(i); 3378 const ptrdiff_t ds = stepi()+stepj(); 3379 return const_vec_type( 3380 itsm.get()+i*(uplo()==Upper?stepj():stepi())+j1*ds, 3381 j2-j1,ds,newct); 3382 } 3383 3384 inline vec_type row(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3385 { 3386 if (I==int(FortranStyle)) { 3387 TMVAssert(i>0 && i<=size()); 3388 --i; 3389 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()); 3390 --j1; 3391 } else { 3392 TMVAssert(i>=0 && i<size()); 3393 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()); 3394 } 3395 if ((uplo()==Upper && i-j1<=0) || (uplo()==Lower && j2-i<=1)) 3396 return vec_type( 3397 itsm.get()+i*stepi()+j1*stepj(), 3398 j2-j1,stepj(),NonConj TMV_FIRSTLAST); 3399 else 3400 return vec_type( 3401 itsm.get()+i*stepj()+j1*stepi(), 3402 j2-j1,stepi(),TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3403 } 3404 3405 inline vec_type col(ptrdiff_t j, ptrdiff_t i1, ptrdiff_t i2) 3406 { 3407 if (I==int(FortranStyle)) { 3408 TMVAssert(j>0 && j<=size()); 3409 --j; 3410 TMVAssert(i1>0 && i1-i2<=0 && i2<=size()); 3411 --i1; 3412 } else { 3413 TMVAssert(j>=0 && j<size()); 3414 TMVAssert(i1>=0 && i1-i2<=0 && i2<=size()); 3415 } 3416 if ((uplo()==Upper && i2-j<=1) || (uplo()==Lower && j-i1<=0)) 3417 return vec_type( 3418 itsm.get()+i1*stepi()+j*stepj(), 3419 i2-i1,stepi(),NonConj TMV_FIRSTLAST); 3420 else 3421 return vec_type( 3422 itsm.get()+i1*stepj()+j*stepi(), 3423 i2-i1,stepj(),TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3424 } 3425 3426 inline vec_type diag() 3427 { 3428 return vec_type( 3429 itsm.get(),size(),stepi()+stepj(),NonConj TMV_FIRSTLAST); 3430 } 3431 3432 inline vec_type diag(ptrdiff_t i) 3433 { 3434 TMVAssert(i>=-size() && i<=size()); 3435 ConjType newct = 3436 ((i>0) == (uplo()==Upper)) ? NonConj : TMV_ConjOf(T,NonConj); 3437 i = std::abs(i); 3438 TMVAssert(i<=size()); 3439 return vec_type( 3440 itsm.get()+i*(uplo()==Upper?stepj():stepi()),size()-i, 3441 stepi()+stepj(),newct TMV_FIRSTLAST); 3442 } 3443 3444 inline vec_type diag(ptrdiff_t i, ptrdiff_t j1, ptrdiff_t j2) 3445 { 3446 TMVAssert(i>=-size() && i<=size()); 3447 if (I==int(FortranStyle)) { 3448 TMVAssert(j1>0 && j1-j2<=0 && j2<=size()-std::abs(i)); 3449 --j1; 3450 } else { 3451 TMVAssert(j1>=0 && j1-j2<=0 && j2<=size()-std::abs(i)); 3452 } 3453 ConjType newct = 3454 ((i>0) == (uplo()==Upper)) ? NonConj : TMV_ConjOf(T,NonConj); 3455 i = std::abs(i); 3456 const ptrdiff_t ds = stepi()+stepj(); 3457 return vec_type( 3458 itsm.get()+i*(uplo()==Upper?stepj():stepi())+j1*ds, 3459 j2-j1,ds,newct TMV_FIRSTLAST); 3460 } 3461 3462 // 3463 // Modifying Functions 3464 // 3465 3466 inline type& setZero() 3467 { 3468 std::fill_n(itsm.get(),itslen,T(0)); 3469 return *this; 3470 } 3471 3472 inline type& setAllTo(const T& x) 3473 { 3474 TMVAssert(TMV_IMAG(x) == RT(0)); 3475 upperTri().setAllTo(x); 3476 return *this; 3477 } 3478 3479 inline type& addToAll(const T& x) 3480 { 3481 TMVAssert(TMV_IMAG(x) == RT(0)); 3482 upperTri().addToAll(x); 3483 return *this; 3484 } 3485 3486 inline type& clip(RT thresh) 3487 { upperTri().clip(thresh); return *this; } 3488 3489 inline type& conjugateSelf() 3490 { if (isComplex(T())) upperTri().conjugateSelf(); return *this; } 3491 3492 inline type& transposeSelf() 3493 { conjugateSelf(); } 3494 3495 inline type& setToIdentity(const T& x=T(1)) 3496 { 3497 TMVAssert(TMV_IMAG(x) == RT(0)); 3498 setZero(); 3499 diag().setAllTo(x); 3500 return *this; 3501 } 3502 3503 inline type& swapRowsCols(ptrdiff_t i1, ptrdiff_t i2) 3504 { view().swapRowsCols(i1,i2); return *this; } 3505 3506 inline type& permuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 3507 { view().permuteRowsCols(p,i1,i2); return *this; } 3508 3509 inline type& reversePermuteRowsCols(const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 3510 { view().reversePermuteRowsCols(p,i1,i2); return *this; } 3511 3512 inline type& permuteRowsCols(const ptrdiff_t* p) 3513 { view().permuteRowsCols(p); return *this; } 3514 3515 inline type& reversePermuteRowsCols(const ptrdiff_t* p) 3516 { view().reversePermuteRowsCols(p); return *this; } 3517 3518 // 3519 // subMatrix 3520 // 3521 3522 inline const_rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) const 3523 { 3524 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 3525 if (I==int(FortranStyle)) { --i1; --j1; } 3526 if ((uplo()==Upper && i2-j1<=1) || (uplo()==Lower && j2-i1<=1)) 3527 return const_rec_type( 3528 itsm.get()+i1*stepi()+j1*stepj(), 3529 i2-i1,j2-j1,stepi(),stepj(),NonConj); 3530 else 3531 return const_rec_type( 3532 itsm.get()+i1*stepj()+j1*stepi(),i2-i1,j2-j1, 3533 stepj(),stepi(),TMV_ConjOf(T,NonConj)); 3534 } 3535 3536 inline const_rec_type subMatrix( 3537 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 3538 { 3539 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 3540 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 3541 if ((uplo()==Upper && i2-j1<=istep) || 3542 (uplo()==Lower && j2-i1<=jstep)) 3543 return const_rec_type( 3544 itsm.get()+i1*stepi()+j1*stepj(), 3545 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(), 3546 NonConj); 3547 else 3548 return const_rec_type( 3549 itsm.get()+i1*stepj()+j1*stepi(), 3550 (i2-i1)/istep,(j2-j1)/jstep,istep*stepj(),jstep*stepi(), 3551 TMV_ConjOf(T,NonConj)); 3552 } 3553 3554 inline const_vec_type subVector( 3555 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 3556 { 3557 TMVAssert(base::hasSubVector(i,j,istep,jstep,n)); 3558 if ((uplo()==Upper && i-j<=0) || (uplo()==Lower && j-i<=0)) 3559 return const_vec_type( 3560 itsm.get()+i*stepi()+j*stepj(),n, 3561 istep*stepi()+jstep*stepj(),NonConj); 3562 else 3563 return const_vec_type( 3564 itsm.get()+i*stepj()+j*stepi(),n, 3565 istep*stepj()+jstep*stepi(),TMV_ConjOf(T,NonConj)); 3566 } 3567 3568 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) const 3569 { 3570 TMVAssert(view().hasSubSymMatrix(i1,i2,1)); 3571 if (I==int(FortranStyle)) { --i1; } 3572 return const_view_type( 3573 itsm.get()+i1*(stepi()+stepj()), 3574 i2-i1,stepi(),stepj(),Herm,uplo(),NonConj); 3575 } 3576 3577 inline const_view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 3578 { 3579 TMVAssert(view().hasSubSymMatrix(i1,i2,istep)); 3580 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 3581 return const_view_type( 3582 itsm.get()+i1*(stepi()+stepj()), 3583 (i2-i1)/istep,istep*stepi(),istep*stepj(),Herm,uplo(), NonConj); 3584 } 3585 3586 inline const_uppertri_type upperTri(DiagType dt = NonUnitDiag) const 3587 { 3588 return uplo()==Upper ? 3589 const_uppertri_type( 3590 itsm.get(),size(),stepi(),stepj(),dt,NonConj) : 3591 const_uppertri_type( 3592 itsm.get(),size(),stepj(),stepi(), 3593 dt,TMV_ConjOf(T,NonConj)); 3594 } 3595 3596 inline const_uppertri_type unitUpperTri() const 3597 { 3598 return uplo()==Upper ? 3599 const_uppertri_type( 3600 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj) : 3601 const_uppertri_type( 3602 itsm.get(),size(),stepj(),stepi(), 3603 UnitDiag,TMV_ConjOf(T,NonConj)); 3604 } 3605 3606 inline const_lowertri_type lowerTri(DiagType dt = NonUnitDiag) const 3607 { 3608 return uplo()==Lower ? 3609 const_lowertri_type( 3610 itsm.get(),size(),stepi(),stepj(),dt,NonConj) : 3611 const_lowertri_type( 3612 itsm.get(),size(),stepj(),stepi(), 3613 dt,TMV_ConjOf(T,NonConj)); 3614 } 3615 3616 inline const_lowertri_type unitLowerTri() const 3617 { 3618 return uplo()==Lower ? 3619 const_lowertri_type( 3620 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj) : 3621 const_lowertri_type( 3622 itsm.get(),size(),stepj(),stepi(), 3623 UnitDiag,TMV_ConjOf(T,NonConj)); 3624 } 3625 3626 inline const_realpart_type realPart() const 3627 { 3628 return const_realpart_type( 3629 reinterpret_cast<const RT*>(itsm.get()),size(), 3630 isReal(T()) ? stepi() : 2*stepi(), 3631 isReal(T()) ? stepj() : 2*stepj(), Herm,uplo(),NonConj); 3632 } 3633 3634 inline const_realpart_type imagPart() const 3635 { 3636 // The imaginary part of a Hermitian matrix is anti-symmetric 3637 // so this is illegal. 3638 TMVAssert(TMV_FALSE); 3639 return const_realpart_type(0,0,0,0,Herm,uplo(),NonConj); 3640 } 3641 3642 inline const_view_type view() const 3643 { 3644 return const_view_type( 3645 itsm.get(),size(),stepi(),stepj(),Herm,uplo(),NonConj); 3646 } 3647 3648 inline const_view_type transpose() const 3649 { 3650 return const_view_type( 3651 itsm.get(),size(),stepj(),stepi(), 3652 Herm,TMV_UTransOf(U),NonConj); 3653 } 3654 3655 inline const_view_type conjugate() const 3656 { 3657 return const_view_type( 3658 itsm.get(),size(),stepi(),stepj(), 3659 Herm,uplo(),TMV_ConjOf(T,NonConj)); 3660 } 3661 3662 inline const_view_type adjoint() const 3663 { 3664 return const_view_type( 3665 itsm.get(),size(),stepj(),stepi(), 3666 Herm,TMV_UTransOf(U),TMV_ConjOf(T,NonConj)); 3667 } 3668 3669 inline rec_type subMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2) 3670 { 3671 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,1,1)); 3672 if (I==int(FortranStyle)) { --i1; --j1; } 3673 if ((uplo()==Upper && i2-j1<=1) || (uplo()==Lower && j2-i1<=1)) 3674 return rec_type( 3675 itsm.get()+i1*stepi()+j1*stepj(), 3676 i2-i1,j2-j1,stepi(),stepj(),NonConj TMV_FIRSTLAST); 3677 else 3678 return rec_type( 3679 itsm.get()+i1*stepj()+j1*stepi(), 3680 i2-i1,j2-j1,stepj(),stepi(), 3681 TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3682 } 3683 3684 inline rec_type subMatrix( 3685 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) 3686 { 3687 TMVAssert(view().hasSubMatrix(i1,i2,j1,j2,istep,jstep)); 3688 if (I==int(FortranStyle)) { --i1; --j1; i2+=istep-1; j2+=jstep-1; } 3689 if ((uplo()==Upper && i2-j1<=istep) || 3690 (uplo()==Lower && j2-i1<=jstep)) 3691 return rec_type( 3692 itsm.get()+i1*stepi()+j1*stepj(), 3693 (i2-i1)/istep,(j2-j1)/jstep,istep*stepi(),jstep*stepj(), 3694 NonConj TMV_FIRSTLAST); 3695 else 3696 return rec_type( 3697 itsm.get()+i1*stepj()+j1*stepi(), 3698 (i2-i1)/istep,(j2-j1)/jstep,istep*stepj(),jstep*stepi(), 3699 TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3700 } 3701 3702 inline vec_type subVector( 3703 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) 3704 { 3705 TMVAssert(base::hasSubVector(i,j,istep,jstep,n)); 3706 if ((uplo()==Upper && i-j<=0) || (uplo()==Lower && j-i<=0)) 3707 return vec_type( 3708 itsm.get()+i*stepi()+j*stepj(),n, 3709 istep*stepi()+jstep*stepj(),NonConj TMV_FIRSTLAST); 3710 else 3711 return vec_type( 3712 itsm.get()+i*stepj()+j*stepi(),n, 3713 istep*stepj()+jstep*stepi(),TMV_ConjOf(T,NonConj) 3714 TMV_FIRSTLAST); 3715 } 3716 3717 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2) 3718 { 3719 TMVAssert(view().hasSubSymMatrix(i1,i2,1)); 3720 if (I==int(FortranStyle)) { --i1; } 3721 return view_type( 3722 itsm.get()+i1*(stepi()+stepj()), 3723 i2-i1,stepi(),stepj(),Herm,uplo(),NonConj TMV_FIRSTLAST); 3724 } 3725 3726 inline view_type subSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) 3727 { 3728 TMVAssert(view().hasSubSymMatrix(i1,i2,istep)); 3729 if (I==int(FortranStyle)) { --i1; i2+=istep-1; } 3730 return view_type( 3731 itsm.get()+i1*(stepi()+stepj()), 3732 (i2-i1)/istep,istep*stepi(),istep*stepj(),Herm,uplo(), NonConj 3733 TMV_FIRSTLAST); 3734 } 3735 3736 inline uppertri_type upperTri(DiagType dt = NonUnitDiag) 3737 { 3738 return uplo()==Upper ? 3739 uppertri_type( 3740 itsm.get(),size(),stepi(),stepj(),dt,NonConj 3741 TMV_FIRSTLAST) : 3742 uppertri_type( 3743 itsm.get(),size(),stepj(),stepi(),dt,TMV_ConjOf(T,NonConj) 3744 TMV_FIRSTLAST); 3745 } 3746 3747 inline uppertri_type unitUpperTri() 3748 { 3749 return uplo()==Upper ? 3750 uppertri_type( 3751 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj 3752 TMV_FIRSTLAST) : 3753 uppertri_type( 3754 itsm.get(),size(),stepj(),stepi(), 3755 UnitDiag,TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3756 } 3757 3758 inline lowertri_type lowerTri(DiagType dt = NonUnitDiag) 3759 { 3760 return uplo()==Lower ? 3761 lowertri_type( 3762 itsm.get(),size(),stepi(),stepj(),dt,NonConj 3763 TMV_FIRSTLAST) : 3764 lowertri_type( 3765 itsm.get(),size(),stepj(),stepi(), 3766 dt,TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3767 } 3768 3769 inline lowertri_type unitLowerTri() 3770 { 3771 return uplo()==Lower ? 3772 lowertri_type( 3773 itsm.get(),size(),stepi(),stepj(),UnitDiag,NonConj 3774 TMV_FIRSTLAST) : 3775 lowertri_type( 3776 itsm.get(),size(),stepj(),stepi(), 3777 UnitDiag,TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3778 } 3779 3780 inline realpart_type realPart() 3781 { 3782 return realpart_type( 3783 reinterpret_cast<RT*>( 3784 itsm.get()),size(), 3785 isReal(T()) ? stepi() : 2*stepi(), 3786 isReal(T()) ? stepj() : 2*stepj(), Herm,uplo(),NonConj 3787 #ifdef TMVFLDEBUG 3788 ,reinterpret_cast<const RT*>(_first) 3789 ,reinterpret_cast<const RT*>(_last) 3790 #endif 3791 ); 3792 } 3793 3794 inline realpart_type imagPart() 3795 { 3796 // The imaginary part of a Hermitian matrix is anti-symmetric 3797 // so this is illegal. 3798 TMVAssert(TMV_FALSE); 3799 return realpart_type(0,0,0,0,Herm,uplo(),NonConj 3800 TMV_FIRSTLAST1(0,0) ); 3801 } 3802 3803 inline view_type view() 3804 { 3805 return view_type( 3806 itsm.get(),size(),stepi(),stepj(),Herm,uplo(),NonConj 3807 TMV_FIRSTLAST); 3808 } 3809 3810 inline view_type transpose() 3811 { 3812 return view_type( 3813 itsm.get(),size(),stepj(),stepi(),Herm,TMV_UTransOf(U),NonConj 3814 TMV_FIRSTLAST); 3815 } 3816 3817 inline view_type conjugate() 3818 { 3819 return view_type( 3820 itsm.get(),size(),stepi(),stepj(), 3821 Herm,uplo(),TMV_ConjOf(T,NonConj) TMV_FIRSTLAST); 3822 } 3823 3824 inline view_type adjoint() 3825 { 3826 return view_type( 3827 itsm.get(),size(),stepj(),stepi(), 3828 Herm,TMV_UTransOf(U),TMV_ConjOf(T,NonConj) 3829 TMV_FIRSTLAST); 3830 } 3831 3832 // 3833 // I/O 3834 // 3835 3836 void read(const TMV_Reader& reader); 3837 3838 inline ptrdiff_t size() const { return itss; } 3839 inline const T* cptr() const { return itsm.get(); } 3840 inline T* ptr() { return itsm.get(); } 3841 inline ptrdiff_t stepi() const { return S==int(RowMajor) ? itss : 1; } 3842 inline ptrdiff_t stepj() const { return S==int(RowMajor) ? 1 : itss; } 3843 inline SymType sym() const { return Herm; } 3844 inline UpLoType uplo() const { return static_cast<UpLoType>(U); } 3845 inline ConjType ct() const { return NonConj; } 3846 inline bool isrm() const { return S==int(RowMajor); } 3847 inline bool iscm() const { return S==int(ColMajor); } 3848 inline bool isconj() const { return false; } 3849 inline bool isherm() const { return true; } 3850 inline bool issym() const { return isReal(T()); } 3851 inline bool isupper() const { return U==int(Upper); } 3852 3853 inline reference ref(ptrdiff_t i, ptrdiff_t j) 3854 { 3855 if ((uplo()==Upper && i <= j) || (uplo()==Lower && i>=j)) 3856 return RefHelper<T>::makeRef( 3857 itsm.get() + i*stepi() + j*stepj(),NonConj); 3858 else 3859 return RefHelper<T>::makeRef( 3860 itsm.get() + j*stepi() + i*stepj(),Conj); 3861 } 3862 3863 inline T cref(ptrdiff_t i, ptrdiff_t j) const 3864 { 3865 if ((uplo()==Upper && i <= j) || (uplo()==Lower && i>=j)) 3866 return itsm.get()[i*stepi() + j*stepj()]; 3867 else 3868 return TMV_CONJ(itsm.get()[j*stepi() + i*stepj()]); 3869 } 3870 3871 inline void resize(ptrdiff_t s) 3872 { 3873 TMVAssert(s >= 0); 3874 itslen = s*s; 3875 itsm.resize(itslen); 3876 itss = s; 3877 DivHelper<T>::resetDivType(); 3878 #ifdef TMVFLDEBUG 3879 _first = itsm.get(); 3880 _last = _first+itslen; 3881 #endif 3882 #ifdef TMV_EXTRA_DEBUG 3883 setAllTo(T(888)); 3884 #else 3885 if (isComplex(T())) diag().imagPart().setZero(); 3886 #endif 3887 } 3888 3889 protected : 3890 3891 ptrdiff_t itslen; 3892 AlignedArray<T> itsm; 3893 ptrdiff_t itss; 3894 3895 #ifdef TMVFLDEBUG 3896 public: 3897 const T* _first; 3898 const T* _last; 3899 #endif 3900 3901 friend void Swap(HermMatrix<T,A>& m1, HermMatrix<T,A>& m2) 3902 { 3903 TMVAssert(m1.size() == m2.size()); 3904 m1.itsm.swapWith(m2.itsm); 3905 #ifdef TMVFLDEBUG 3906 TMV_SWAP(m1._first,m2._first); 3907 TMV_SWAP(m1._last,m2._last); 3908 #endif 3909 } 3910 3911 }; // HermMatrix 3912 3913 //------------------------------------------------------------------------- 3914 3915 // 3916 // Special Creators: 3917 // SymMatrixViewOf(m,uplo) 3918 // HermMatrixViewOf(m,uplo) 3919 // SymMatrixViewOf(mptr,size,uplo,stor) 3920 // HermMatrixViewOf(mptr,size,uplo,stor) 3921 // SymMatrixViewOf(mptr,size,uplo,stepi,stepj) 3922 // HermMatrixViewOf(mptr,size,uplo,stepi,stepj) 3923 // 3924 3925 template <typename T> 3926 inline ConstSymMatrixView<T> SymMatrixViewOf( 3927 const GenMatrix<T>& m, UpLoType uplo) 3928 { 3929 TMVAssert(m.colsize()==m.rowsize()); 3930 return ConstSymMatrixView<T>( 3931 m.cptr(),m.rowsize(),m.stepi(),m.stepj(),Sym,uplo,m.ct()); 3932 } 3933 3934 template <typename T, int A> 3935 inline ConstSymMatrixView<T,A> SymMatrixViewOf( 3936 const ConstMatrixView<T,A>& m, UpLoType uplo) 3937 { 3938 TMVAssert(m.colsize()==m.rowsize()); 3939 return ConstSymMatrixView<T,A>( 3940 m.cptr(),m.rowsize(),m.stepi(),m.stepj(),Sym,uplo,m.ct()); 3941 } 3942 3943 template <typename T, int A> 3944 inline ConstSymMatrixView<T,A&FortranStyle> SymMatrixViewOf( 3945 const Matrix<T,A>& m, UpLoType uplo) 3946 { 3947 TMVAssert(m.colsize()==m.rowsize()); 3948 return ConstSymMatrixView<T,A&FortranStyle>( 3949 m.cptr(),m.rowsize(),m.stepi(),m.stepj(),Sym,uplo,m.ct()); 3950 } 3951 3952 template <typename T, int A> 3953 inline SymMatrixView<T,A> SymMatrixViewOf( 3954 MatrixView<T,A> m, UpLoType uplo) 3955 { 3956 TMVAssert(m.colsize()==m.rowsize()); 3957 return SymMatrixView<T,A>( 3958 m.ptr(),m.rowsize(),m.stepi(),m.stepj(),Sym,uplo,m.ct() 3959 TMV_FIRSTLAST1(m._first,m._last)); 3960 } 3961 3962 template <typename T, int A> 3963 inline SymMatrixView<T,A&FortranStyle> SymMatrixViewOf( 3964 Matrix<T,A>& m, UpLoType uplo) 3965 { 3966 TMVAssert(m.colsize()==m.rowsize()); 3967 return SymMatrixView<T,A&FortranStyle>( 3968 m.ptr(),m.rowsize(),m.stepi(),m.stepj(),Sym,uplo,m.ct() 3969 TMV_FIRSTLAST1(m._first,m._last)); 3970 } 3971 3972 template <typename T> 3973 inline ConstSymMatrixView<T> HermMatrixViewOf( 3974 const GenMatrix<T>& m, UpLoType uplo) 3975 { 3976 TMVAssert(m.colsize()==m.rowsize()); 3977 TMVAssert(isReal(T()) || 3978 m.diag().imagPart().normInf() == TMV_RealType(T)(0)); 3979 return ConstSymMatrixView<T>( 3980 m.cptr(),m.rowsize(),m.stepi(),m.stepj(),Herm,uplo,m.ct()); 3981 } 3982 3983 template <typename T, int A> 3984 inline ConstSymMatrixView<T,A> HermMatrixViewOf( 3985 const ConstMatrixView<T,A>& m, UpLoType uplo) 3986 { 3987 TMVAssert(m.colsize()==m.rowsize()); 3988 TMVAssert(isReal(T()) || 3989 m.diag().imagPart().normInf() == TMV_RealType(T)(0)); 3990 return ConstSymMatrixView<T,A>( 3991 m.cptr(),m.rowsize(),m.stepi(),m.stepj(),Herm,uplo,m.ct()); 3992 } 3993 3994 template <typename T, int A> 3995 inline ConstSymMatrixView<T,A&FortranStyle> HermMatrixViewOf( 3996 const Matrix<T,A>& m, UpLoType uplo) 3997 { 3998 TMVAssert(m.colsize()==m.rowsize()); 3999 TMVAssert(isReal(T()) || 4000 m.diag().imagPart().normInf() == TMV_RealType(T)(0)); 4001 return ConstSymMatrixView<T,A&FortranStyle>( 4002 m.cptr(),m.rowsize(),m.stepi(),m.stepj(),Herm,uplo,m.ct()); 4003 } 4004 4005 template <typename T, int A> 4006 inline SymMatrixView<T,A> HermMatrixViewOf( 4007 MatrixView<T,A> m, UpLoType uplo) 4008 { 4009 TMVAssert(m.colsize()==m.rowsize()); 4010 TMVAssert(isReal(T()) || 4011 m.diag().imagPart().normInf() == TMV_RealType(T)(0)); 4012 return SymMatrixView<T,A>( 4013 m.ptr(),m.rowsize(),m.stepi(),m.stepj(),Herm,uplo,m.ct() 4014 TMV_FIRSTLAST1(m._first,m._last)); 4015 } 4016 4017 template <typename T, int A> 4018 inline SymMatrixView<T,A&FortranStyle> HermMatrixViewOf( 4019 Matrix<T,A>& m, UpLoType uplo) 4020 { 4021 TMVAssert(m.colsize()==m.rowsize()); 4022 TMVAssert(isReal(T()) || 4023 m.diag().imagPart().normInf() == TMV_RealType(T)(0)); 4024 return SymMatrixView<T,A&FortranStyle>( 4025 m.ptr(),m.rowsize(),m.stepi(),m.stepj(),Herm,uplo,m.ct() 4026 TMV_FIRSTLAST1(m._first,m._last)); 4027 } 4028 4029 template <typename T> 4030 inline ConstSymMatrixView<T> SymMatrixViewOf( 4031 const T* m, ptrdiff_t size, UpLoType uplo, StorageType stor) 4032 { 4033 TMVAssert(stor == RowMajor || stor == ColMajor); 4034 TMVAssert(size>=0); 4035 const ptrdiff_t stepi = stor == RowMajor ? size : 1; 4036 const ptrdiff_t stepj = stor == RowMajor ? 1 : size; 4037 return ConstSymMatrixView<T>( 4038 m,size,stepi,stepj,Sym,uplo,NonConj); 4039 } 4040 4041 template <typename T> 4042 inline ConstSymMatrixView<T> HermMatrixViewOf( 4043 const T* m, ptrdiff_t size, UpLoType uplo, StorageType stor) 4044 { 4045 TMVAssert(stor == RowMajor || stor == ColMajor); 4046 TMVAssert(size>=0); 4047 const ptrdiff_t stepi = stor == RowMajor ? size : 1; 4048 const ptrdiff_t stepj = stor == RowMajor ? 1 : size; 4049 return ConstSymMatrixView<T>(m,size,stepi,stepj,Herm,uplo,NonConj); 4050 } 4051 4052 template <typename T> 4053 inline SymMatrixView<T> SymMatrixViewOf( 4054 T* m, ptrdiff_t size, UpLoType uplo, StorageType stor) 4055 { 4056 TMVAssert(stor == RowMajor || stor == ColMajor); 4057 TMVAssert(size>=0); 4058 const ptrdiff_t stepi = stor == RowMajor ? size : 1; 4059 const ptrdiff_t stepj = stor == RowMajor ? 1 : size; 4060 return SymMatrixView<T>( 4061 m,size,stepi,stepj,Sym,uplo,NonConj TMV_FIRSTLAST1(m,m+size*size)); 4062 } 4063 4064 template <typename T> 4065 inline SymMatrixView<T> HermMatrixViewOf( 4066 T* m, ptrdiff_t size, UpLoType uplo, StorageType stor) 4067 { 4068 TMVAssert(stor == RowMajor || stor == ColMajor); 4069 TMVAssert(size>=0); 4070 const ptrdiff_t stepi = stor == RowMajor ? size : 1; 4071 const ptrdiff_t stepj = stor == RowMajor ? 1 : size; 4072 return SymMatrixView<T>( 4073 m,size,stepi,stepj,Herm,uplo,NonConj TMV_FIRSTLAST1(m,m+size*size)); 4074 } 4075 4076 template <typename T> 4077 inline ConstSymMatrixView<T> SymMatrixViewOf( 4078 const T* m, ptrdiff_t size, UpLoType uplo, ptrdiff_t stepi, ptrdiff_t stepj) 4079 { 4080 TMVAssert(size>=0); 4081 return ConstSymMatrixView<T>( 4082 m,size,stepi,stepj,Sym,uplo,NonConj); 4083 } 4084 4085 template <typename T> 4086 inline ConstSymMatrixView<T> HermMatrixViewOf( 4087 const T* m, ptrdiff_t size, UpLoType uplo, ptrdiff_t stepi, ptrdiff_t stepj) 4088 { 4089 TMVAssert(size>=0); 4090 return ConstSymMatrixView<T>( 4091 m,size,stepi,stepj,Herm,uplo,NonConj); 4092 } 4093 4094 template <typename T> 4095 inline SymMatrixView<T> SymMatrixViewOf( 4096 T* m, ptrdiff_t size, UpLoType uplo, ptrdiff_t stepi, ptrdiff_t stepj) 4097 { 4098 TMVAssert(size>=0); 4099 return SymMatrixView<T>( 4100 m,size,stepi,stepj,Sym,uplo,NonConj TMV_FIRSTLAST1(m,m+size*size)); 4101 } 4102 4103 template <typename T> 4104 inline SymMatrixView<T> HermMatrixViewOf( 4105 T* m, ptrdiff_t size, UpLoType uplo, ptrdiff_t stepi, ptrdiff_t stepj) 4106 { 4107 TMVAssert(size>=0); 4108 return SymMatrixView<T>( 4109 m,size,stepi,stepj,Herm,uplo,NonConj TMV_FIRSTLAST1(m,m+size*size)); 4110 } 4111 4112 // 4113 // Swap Matrices 4114 // 4115 4116 template <typename T> 4117 inline void Swap(SymMatrixView<T> m1, SymMatrixView<T> m2) 4118 { 4119 TMVAssert(m1.size() == m2.size()); 4120 TMVAssert(m1.issym() == m2.issym()); 4121 TMVAssert(m1.isherm() == m2.isherm()); 4122 Swap(m1.upperTri(),m2.upperTri()); 4123 } 4124 4125 template <typename T, int A> 4126 inline void Swap(SymMatrixView<T> m1, SymMatrix<T,A>& m2) 4127 { Swap(m1,m2.view()); } 4128 4129 template <typename T, int A> 4130 inline void Swap(SymMatrix<T,A>& m1, SymMatrixView<T> m2) 4131 { Swap(m1.view(),m2); } 4132 4133 template <typename T, int A1, int A2> 4134 inline void Swap(SymMatrix<T,A1>& m1, SymMatrix<T,A2>& m2) 4135 { Swap(m1.view(),m2.view()); } 4136 4137 template <typename T, int A> 4138 inline void Swap(SymMatrixView<T> m1, HermMatrix<T,A>& m2) 4139 { Swap(m1,m2.view()); } 4140 4141 template <typename T, int A> 4142 inline void Swap(HermMatrix<T,A>& m1, SymMatrixView<T> m2) 4143 { Swap(m1.view(),m2); } 4144 4145 template <typename T, int A1, int A2> 4146 inline void Swap(HermMatrix<T,A1>& m1, HermMatrix<T,A2>& m2) 4147 { Swap(m1.view(),m2.view()); } 4148 4149 4150 4151 // 4152 // Views: 4153 // 4154 4155 template <typename T> 4156 inline ConstSymMatrixView<T> Transpose(const GenSymMatrix<T>& m) 4157 { return m.transpose(); } 4158 4159 template <typename T, int A> 4160 inline ConstSymMatrixView<T,A> Transpose(const ConstSymMatrixView<T,A>& m) 4161 { return m.transpose(); } 4162 4163 template <typename T, int A> 4164 inline ConstSymMatrixView<T,A&FortranStyle> Transpose( 4165 const SymMatrix<T,A>& m) 4166 { return m.transpose(); } 4167 4168 template <typename T, int A> 4169 inline ConstSymMatrixView<T,A&FortranStyle> Transpose( 4170 const HermMatrix<T,A>& m) 4171 { return m.transpose(); } 4172 4173 template <typename T, int A> 4174 inline SymMatrixView<T,A> Transpose(SymMatrixView<T,A> m) 4175 { return m.transpose(); } 4176 4177 template <typename T, int A> 4178 inline SymMatrixView<T,A&FortranStyle> Transpose(SymMatrix<T,A>& m) 4179 { return m.transpose(); } 4180 4181 template <typename T, int A> 4182 inline SymMatrixView<T,A&FortranStyle> Transpose(HermMatrix<T,A>& m) 4183 { return m.transpose(); } 4184 4185 template <typename T> 4186 inline ConstSymMatrixView<T> Conjugate(const GenSymMatrix<T>& m) 4187 { return m.conjugate(); } 4188 4189 template <typename T, int A> 4190 inline ConstSymMatrixView<T,A> Conjugate(const ConstSymMatrixView<T,A>& m) 4191 { return m.conjugate(); } 4192 4193 template <typename T, int A> 4194 inline ConstSymMatrixView<T,A&FortranStyle> Conjugate( 4195 const SymMatrix<T,A>& m) 4196 { return m.conjugate(); } 4197 4198 template <typename T, int A> 4199 inline ConstSymMatrixView<T,A&FortranStyle> Conjugate( 4200 const HermMatrix<T,A>& m) 4201 { return m.conjugate(); } 4202 4203 template <typename T, int A> 4204 inline SymMatrixView<T,A> Conjugate(SymMatrixView<T,A> m) 4205 { return m.conjugate(); } 4206 4207 template <typename T, int A> 4208 inline SymMatrixView<T,A&FortranStyle> Conjugate(SymMatrix<T,A>& m) 4209 { return m.conjugate(); } 4210 4211 template <typename T, int A> 4212 inline SymMatrixView<T,A&FortranStyle> Conjugate(HermMatrix<T,A>& m) 4213 { return m.conjugate(); } 4214 4215 template <typename T> 4216 inline ConstSymMatrixView<T> Adjoint(const GenSymMatrix<T>& m) 4217 { return m.adjoint(); } 4218 4219 template <typename T, int A> 4220 inline ConstSymMatrixView<T,A> Adjoint(const ConstSymMatrixView<T,A>& m) 4221 { return m.adjoint(); } 4222 4223 template <typename T, int A> 4224 inline ConstSymMatrixView<T,A&FortranStyle> Adjoint( 4225 const SymMatrix<T,A>& m) 4226 { return m.adjoint(); } 4227 4228 template <typename T, int A> 4229 inline ConstSymMatrixView<T,A&FortranStyle> Adjoint( 4230 const HermMatrix<T,A>& m) 4231 { return m.adjoint(); } 4232 4233 template <typename T, int A> 4234 inline SymMatrixView<T,A> Adjoint(SymMatrixView<T,A> m) 4235 { return m.adjoint(); } 4236 4237 template <typename T, int A> 4238 inline SymMatrixView<T,A&FortranStyle> Adjoint(SymMatrix<T,A>& m) 4239 { return m.adjoint(); } 4240 4241 template <typename T, int A> 4242 inline SymMatrixView<T,A&FortranStyle> Adjoint(HermMatrix<T,A>& m) 4243 { return m.adjoint(); } 4244 4245 template <typename T> 4246 inline QuotXS<T,T> Inverse(const GenSymMatrix<T>& m) 4247 { return m.inverse(); } 4248 4249 // 4250 // SymMatrix ==, != SymMatrix 4251 // 4252 4253 template <typename T1, typename T2> 4254 inline bool operator==( 4255 const GenSymMatrix<T1>& m1, const GenSymMatrix<T2>& m2) 4256 { return m1.upperTri() == m2.upperTri(); } 4257 4258 template <typename T1, typename T2> 4259 inline bool operator!=( 4260 const GenSymMatrix<T1>& m1, const GenSymMatrix<T2>& m2) 4261 { return !(m1 == m2); } 4262 4263 template <typename T1, typename T2> 4264 inline bool operator==( 4265 const GenSymMatrix<T1>& m1, const GenMatrix<T2>& m2) 4266 { 4267 return 4268 m1.upperTri() == m2.upperTri() && 4269 m1.lowerTri() == m2.lowerTri(); 4270 } 4271 4272 template <typename T1, typename T2> 4273 inline bool operator==( 4274 const GenMatrix<T1>& m1, const GenSymMatrix<T2>& m2) 4275 { return m2 == m1; } 4276 4277 template <typename T1, typename T2> 4278 inline bool operator!=( 4279 const GenSymMatrix<T1>& m1, const GenMatrix<T2>& m2) 4280 { return !(m1 == m2); } 4281 4282 template <typename T1, typename T2> 4283 inline bool operator!=( 4284 const GenMatrix<T1>& m1, const GenSymMatrix<T2>& m2) 4285 { return !(m1 == m2); } 4286 4287 4288 // 4289 // I/O 4290 // 4291 4292 template <typename T> 4293 inline std::istream& operator>>(std::istream& is, SymMatrixView<T> m) 4294 { return is >> IOStyle() >> m; } 4295 4296 template <typename T, int A> 4297 inline std::istream& operator>>(std::istream& is, SymMatrix<T,A>& m) 4298 { return is >> IOStyle() >> m; } 4299 4300 template <typename T, int A> 4301 inline std::istream& operator>>(std::istream& is, HermMatrix<T,A>& m) 4302 { return is >> IOStyle() >> m; } 4303 4304 template <typename T> 4305 inline std::istream& operator>>( 4306 const TMV_Reader& reader, SymMatrixView<T> m) 4307 { m.read(reader); return reader.getis(); } 4308 4309 template <typename T, int A> 4310 inline std::istream& operator>>( 4311 const TMV_Reader& reader, SymMatrix<T,A>& m) 4312 { m.read(reader); return reader.getis(); } 4313 4314 template <typename T, int A> 4315 inline std::istream& operator>>( 4316 const TMV_Reader& reader, HermMatrix<T,A>& m) 4317 { m.read(reader); return reader.getis(); } 4318 4319 } // namespace tmv 4320 4321 #endif 4322