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 //#define XDEBUG 24 25 26 #include "TMV_Blas.h" 27 #include "tmv/TMV_SymMatrix.h" 28 #include "tmv/TMV_Matrix.h" 29 #include "tmv/TMV_Vector.h" 30 #include "tmv/TMV_SymLDLD.h" 31 #include "tmv/TMV_SymCHD.h" 32 #include "tmv/TMV_SymSVD.h" 33 #include "tmv/TMV_VIt.h" 34 #include "tmv/TMV_SymMatrixArith.h" 35 #include "tmv/TMV_DiagMatrix.h" 36 #include "TMV_IntegerDet.h" 37 #include <iostream> 38 39 #ifdef XDEBUG 40 using std::cerr; 41 using std::endl; 42 #endif 43 44 namespace tmv { 45 46 #define RT TMV_RealType(T) 47 48 // 49 // Access 50 // 51 52 template <class T> cref(ptrdiff_t i,ptrdiff_t j) const53 T GenSymMatrix<T>::cref(ptrdiff_t i, ptrdiff_t j) const 54 { 55 if ((uplo() == Upper && i<=j) || (uplo() == Lower && i>=j)) { 56 const T* mi = cptr() + i*stepi() + j*stepj(); 57 return isconj() ? TMV_CONJ(*mi) : *mi; 58 } else { 59 const T* mi = cptr() + j*stepi() + i*stepj(); 60 return issym() != isconj() ? *mi : TMV_CONJ(*mi); 61 } 62 } 63 64 template <class T, int A> ref(ptrdiff_t i,ptrdiff_t j)65 typename SymMatrixView<T,A>::reference SymMatrixView<T,A>::ref( 66 ptrdiff_t i, ptrdiff_t j) 67 { 68 if ((uplo() == Upper && i<=j) || (uplo() == Lower && i>=j)) { 69 T* mi = ptr() + i*stepi() + j*stepj(); 70 return RefHelper<T>::makeRef(mi,ct()); 71 } else { 72 T* mi = ptr() + j*stepi() + i*stepj(); 73 return RefHelper<T>::makeRef( 74 mi, this->issym() != this->isconj() ? NonConj : Conj); 75 } 76 } 77 78 template <class T> setDiv() const79 void GenSymMatrix<T>::setDiv() const 80 { 81 if (!this->divIsSet()) { 82 DivType dt = this->getDivType(); 83 TMVAssert(dt == tmv::LU || dt == tmv::CH || dt == tmv::SV); 84 TMVAssert(isherm() || dt != tmv::CH); 85 switch (dt) { 86 case LU : 87 this->divider.reset( 88 new SymLDLDiv<T>(*this,this->divIsInPlace())); 89 break; 90 case CH : 91 this->divider.reset( 92 new HermCHDiv<T>(*this,this->divIsInPlace())); 93 break; 94 case SV : 95 if (isherm()) 96 this->divider.reset( 97 new HermSVDiv<T>(*this,this->divIsInPlace())); 98 else 99 this->divider.reset( 100 new SymSVDiv<T>(*this,this->divIsInPlace())); 101 break; 102 default : 103 // The above assert should have already failed 104 // so go ahead and fall through. 105 break; 106 } 107 } 108 } 109 110 #ifdef INST_INT 111 template <> setDiv() const112 void GenSymMatrix<int>::setDiv() const 113 { TMVAssert(TMV_FALSE); } 114 template <> setDiv() const115 void GenSymMatrix<std::complex<int> >::setDiv() const 116 { TMVAssert(TMV_FALSE); } 117 #endif 118 119 template <class T> QInverse() const120 QuotXS<T,T> GenSymMatrix<T>::QInverse() const 121 { return QuotXS<T,T>(T(1),*this); } 122 123 template <class T> template <class T1> doMakeInverse(SymMatrixView<T1> sinv) const124 void GenSymMatrix<T>::doMakeInverse(SymMatrixView<T1> sinv) const 125 { 126 TMVAssert(issym() == sinv.issym()); 127 TMVAssert(isherm() == sinv.isherm()); 128 129 this->setDiv(); 130 const SymDivider<T>* sdiv = dynamic_cast<const SymDivider<T>*>( 131 this->getDiv()); 132 TMVAssert(sdiv); 133 sdiv->makeInverse(sinv); 134 } 135 136 template <class T> divIsLUDiv() const137 bool GenSymMatrix<T>::divIsLUDiv() const 138 { return dynamic_cast<const SymLDLDiv<T>*>(this->getDiv()); } 139 140 template <class T> divIsCHDiv() const141 bool GenSymMatrix<T>::divIsCHDiv() const 142 { return dynamic_cast<const HermCHDiv<T>*>(this->getDiv()); } 143 144 template <class T> divIsHermSVDiv() const145 bool GenSymMatrix<T>::divIsHermSVDiv() const 146 { return dynamic_cast<const HermSVDiv<T>*>(this->getDiv()); } 147 148 template <class T> divIsSymSVDiv() const149 bool GenSymMatrix<T>::divIsSymSVDiv() const 150 { return dynamic_cast<const SymSVDiv<T>*>(this->getDiv()); } 151 152 #ifdef INST_INT 153 template <> divIsLUDiv() const154 bool GenSymMatrix<int>::divIsLUDiv() const 155 { return false; } 156 template <> divIsCHDiv() const157 bool GenSymMatrix<int>::divIsCHDiv() const 158 { return false; } 159 template <> divIsHermSVDiv() const160 bool GenSymMatrix<int>::divIsHermSVDiv() const 161 { return false; } 162 template <> divIsSymSVDiv() const163 bool GenSymMatrix<int>::divIsSymSVDiv() const 164 { return false; } 165 166 template <> divIsLUDiv() const167 bool GenSymMatrix<std::complex<int> >::divIsLUDiv() const 168 { return false; } 169 template <> divIsCHDiv() const170 bool GenSymMatrix<std::complex<int> >::divIsCHDiv() const 171 { return false; } 172 template <> divIsHermSVDiv() const173 bool GenSymMatrix<std::complex<int> >::divIsHermSVDiv() const 174 { return false; } 175 template <> divIsSymSVDiv() const176 bool GenSymMatrix<std::complex<int> >::divIsSymSVDiv() const 177 { return false; } 178 #endif 179 180 // 181 // OK? (SubMatrix, etc.) 182 // 183 184 template <class T> hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const185 bool GenSymMatrix<T>::hasSubMatrix( 186 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 187 { 188 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 189 bool ok = true; 190 ptrdiff_t i2x = i2-istep; 191 ptrdiff_t j2x = j2-jstep; 192 if (istep == 0) { 193 ok = false; 194 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 195 } 196 if (i1 < 0 || i1 >= size()) { 197 ok = false; 198 std::cerr<<"first col element ("<<i1<<") must be in 0 -- "; 199 std::cerr<<size()-1<<std::endl; 200 } 201 if (i2x < 0 || i2x >= size()) { 202 ok = false; 203 std::cerr<<"last col element ("<<i2x<<") must be in 0 -- "; 204 std::cerr<<size()-1<<std::endl; 205 } 206 if ((i2-i1)%istep != 0) { 207 ok = false; 208 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 209 std::cerr<<istep<<")\n"; 210 } 211 if ((i2-i1)/istep < 0) { 212 ok = false; 213 std::cerr<<"n col elements ("<<(i2-i1)/istep<<") must be nonnegative\n"; 214 } 215 if (jstep == 0) { 216 ok = false; 217 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 218 } 219 if (j1 < 0 || j1 >= size()) { 220 ok = false; 221 std::cerr<<"first row element ("<<j1<<") must be in 0 -- "; 222 std::cerr<<size()-1<<std::endl; 223 } 224 if (j2-jstep < 0 || j2-jstep >= size()) { 225 ok = false; 226 std::cerr<<"last row element ("<<j2-jstep<<") must be in 0 -- "; 227 std::cerr<<size()-1<<std::endl; 228 } 229 if ((j2-j1)%jstep != 0) { 230 ok = false; 231 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 232 std::cerr<<jstep<<")\n"; 233 } 234 if ((j2-j1)/jstep < 0) { 235 ok = false; 236 std::cerr<<"n row elements ("<<(j2-j1)/jstep<<") must be nonnegative\n"; 237 } 238 if ((i1<j1 && i2x>j2x) || (i1>j1 && i2x<j2x)) { 239 ok = false; 240 std::cerr<<"Upper left ("<<i1<<','<<j1<<") and lower right ("; 241 std::cerr<<i2x<<','<<j2x<<") corners must be in same triangle\n"; 242 } 243 if ((i2x<j1 && i1>j2x) || (i2x>j1 && i1<j2x)) { 244 ok = false; 245 std::cerr<<"Upper right ("<<i1<<','<<j2x<<") and lower left ("; 246 std::cerr<<i2x<<','<<j1<<") corners must be in same triangle\n"; 247 } 248 return ok; 249 } 250 251 template <class T> hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t n) const252 bool GenSymMatrix<T>::hasSubVector( 253 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 254 { 255 if (n==0) return true; 256 bool ok = true; 257 if (istep == 0 && jstep == 0) { 258 ok = false; 259 std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep; 260 std::cerr<<") can not both be 0\n"; 261 } 262 if (i<0 || i >= size()) { 263 ok = false; 264 std::cerr<<"i ("<<i<<") must be in 0 -- "<<size()-1<<std::endl; 265 } 266 if (j<0 || j >= size()) { 267 ok = false; 268 std::cerr<<"j ("<<j<<") must be in 0 -- "<<size()-1<<std::endl; 269 } 270 ptrdiff_t i2 = i+istep*(n-1); 271 ptrdiff_t j2 = j+jstep*(n-1); 272 if (i2 < 0 || i2 >= size()) { 273 ok = false; 274 std::cerr<<"last element's i ("<<i2<<") must be in 0 -- "; 275 std::cerr<<size()-1<<std::endl; 276 } 277 if (j2 < 0 || j2 >= size()) { 278 ok = false; 279 std::cerr<<"last element's j ("<<j2<<") must be in 0 -- "; 280 std::cerr<<size()-1<<std::endl; 281 } 282 if ((i<j && i2>j2) || (i>j && i2<j2)) { 283 ok = false; 284 std::cerr<<"First ("<<i<<','<<j<<") and last ("<<i2<<','<<j2; 285 std::cerr<<") elements must be in same triangle\n"; 286 } 287 return ok; 288 } 289 290 template <class T> hasSubSymMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep) const291 bool GenSymMatrix<T>::hasSubSymMatrix(ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 292 { 293 if (i1==i2) return true; 294 bool ok=true; 295 if (istep == 0) { 296 ok = false; 297 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 298 } 299 if (i1<0 || i1 >= size()) { 300 ok = false; 301 std::cerr<<"first diag element ("<<i1<<") must be in 0 -- "; 302 std::cerr<<size()-1<<std::endl; 303 } 304 if (i2-istep<0 || i2-istep >= size()) { 305 ok = false; 306 std::cerr<<"last diag element ("<<i2-istep<<") must be in 0 -- "; 307 std::cerr<<size()-1<<std::endl; 308 } 309 if ((i2-i1)%istep != 0) { 310 ok = false; 311 std::cerr<<"range ("<<i2-i1<<") must be multiple of istep ("; 312 std::cerr<<istep<<")\n"; 313 } 314 if ((i2-i1)/istep < 0) { 315 ok = false; 316 std::cerr<<"n diag elements ("<<(i2-i1)/istep<<") must be nonnegative\n"; 317 } 318 return ok; 319 } 320 321 template <class T> hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const322 bool ConstSymMatrixView<T,FortranStyle>::hasSubMatrix( 323 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 324 { 325 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 326 bool ok = true; 327 if (istep == 0) { 328 ok = false; 329 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 330 } 331 if (i1 < 1 || i1 > this->size()) { 332 ok = false; 333 std::cerr<<"first col element ("<<i1<<") must be in 1 -- "; 334 std::cerr<<this->size()<<std::endl; 335 } 336 if (i2 < 1 || i2 > this->size()) { 337 ok = false; 338 std::cerr<<"last col element ("<<i2-istep<<") must be in 1 -- "; 339 std::cerr<<this->size()<<std::endl; 340 } 341 if ((i2-i1)%istep != 0) { 342 ok = false; 343 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 344 std::cerr<<istep<<")\n"; 345 } 346 if ((i2-i1)/istep < 0) { 347 ok = false; 348 std::cerr<<"n col elements ("<<(i2-i1)/istep+1<<") must be positive\n"; 349 } 350 if (jstep == 0) { 351 ok = false; 352 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 353 } 354 if (j1 < 1 || j1 > this->size()) { 355 ok = false; 356 std::cerr<<"first row element ("<<j1<<") must be in 1 -- "; 357 std::cerr<<this->size()<<std::endl; 358 } 359 if (j2 < 1 || j2 > this->size()) { 360 ok = false; 361 std::cerr<<"last row element ("<<j2-jstep<<") must be in 1 -- "; 362 std::cerr<<this->size()<<std::endl; 363 } 364 if ((j2-j1)%jstep != 0) { 365 ok = false; 366 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 367 std::cerr<<jstep<<")\n"; 368 } 369 if ((j2-j1)/jstep < 0) { 370 ok = false; 371 std::cerr<<"n row elements ("<<(j2-j1)/jstep+1<<") must be positive\n"; 372 } 373 if ((i1<j1 && i2>j2) || (i1>j1 && i2<j2)) { 374 ok = false; 375 std::cerr<<"Upper left ("<<i1<<','<<j1<<") and lower right ("; 376 std::cerr<<i2<<','<<j2<<") corners must be in same triangle\n"; 377 } 378 if ((i2<j1 && i1>j2) || (i2>j1 && i1<j2)) { 379 ok = false; 380 std::cerr<<"Upper right ("<<i1<<','<<j2<<") and lower left ("; 381 std::cerr<<i2<<','<<j1<<") corners must be in same triangle\n"; 382 } 383 return ok; 384 } 385 386 template <class T> hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t n) const387 bool ConstSymMatrixView<T,FortranStyle>::hasSubVector( 388 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t n) const 389 { 390 if (n==0) return true; 391 bool ok = true; 392 if (istep == 0 && jstep == 0) { 393 ok = false; 394 std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep; 395 std::cerr<<") can not both be 0\n"; 396 } 397 if (i < 1 || i > this->size()) { 398 ok = false; 399 std::cerr<<"i ("<<i<<") must be in 1 -- "<<this->size()<<std::endl; 400 } 401 if (i < 1 || j > this->size()) { 402 ok = false; 403 std::cerr<<"j ("<<j<<") must be in 1 -- "<<this->size()<<std::endl; 404 } 405 ptrdiff_t i2 = i+istep*(n-1); 406 ptrdiff_t j2 = j+jstep*(n-1); 407 if (i2 < 1 || i2 > this->size()) { 408 ok = false; 409 std::cerr<<"last element's i ("<<i2<<") must be in 1 -- "; 410 std::cerr<<this->size()<<std::endl; 411 } 412 if (j2 < 1 || j2 > this->size()) { 413 ok = false; 414 std::cerr<<"last element's j ("<<j2<<") must be in 1 -- "; 415 std::cerr<<this->size()<<std::endl; 416 } 417 if ((i<j && i2>j2) || (i>j && i2<j2)) { 418 ok = false; 419 std::cerr<<"First ("<<i<<','<<j<<") and last ("<<i2<<','<<j2; 420 std::cerr<<") elements must be in same triangle\n"; 421 } 422 return ok; 423 } 424 425 template <class T> hasSubSymMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t istep) const426 bool ConstSymMatrixView<T,FortranStyle>::hasSubSymMatrix( 427 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t istep) const 428 { 429 if (i1==i2) return true; 430 bool ok=true; 431 if (istep == 0) { 432 ok = false; 433 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 434 } 435 if (i1<1 || i1 > this->size()) { 436 ok = false; 437 std::cerr<<"first diag element ("<<i1<<") must be in 1 -- "; 438 std::cerr<<this->size()<<std::endl; 439 } 440 if (i2-istep<1 || i2-istep > this->size()) { 441 ok = false; 442 std::cerr<<"last diag element ("<<i2-istep<<") must be in 1 -- "; 443 std::cerr<<this->size()<<std::endl; 444 } 445 if ((i2-i1)%istep != 0) { 446 ok = false; 447 std::cerr<<"range ("<<i2-i1<<") must be multiple of istep ("<<istep; 448 std::cerr<<")\n"; 449 } 450 if ((i2-i1)/istep < 0) { 451 ok = false; 452 std::cerr<<"n diag elements ("<<(i2-i1)/istep+1<<") must be positive\n"; 453 } 454 return ok; 455 } 456 457 // 458 // SwapRowsCols 459 // 460 461 template <class T, int A> swapRowsCols(ptrdiff_t i1,ptrdiff_t i2)462 SymMatrixView<T,A>& SymMatrixView<T,A>::swapRowsCols(ptrdiff_t i1, ptrdiff_t i2) 463 { 464 TMVAssert(i1<size()); 465 TMVAssert(i2<size()); 466 if (i1 == i2) return *this; 467 else { 468 #ifdef XDEBUG 469 Matrix<T> M(*this); 470 Matrix<T> M0(*this); 471 M.swapRows(i1,i2); 472 M.swapCols(i1,i2); 473 #endif 474 if (i1 > i2) TMV_SWAP(i1,i2); 475 // Now i1 < i2 476 if (uplo() == Upper) transpose().swapRowsCols(i1,i2); 477 else { 478 Swap(row(i1,0,i1),row(i2,0,i1)); 479 Swap(row(i2,i1+1,i2),col(i1,i1+1,i2)); 480 if (!this->issym()) { 481 row(i2,i1,i2).conjugateSelf(); // Conj m(i2,i1) as well 482 col(i1,i1+1,i2).conjugateSelf(); 483 } 484 Swap(col(i1,i2+1,size()),col(i2,i2+1,size())); 485 diag().swap(i1,i2); 486 } 487 #ifdef XDEBUG 488 if (Norm(M-*this) > 1.e-5*TMV_MAX(RT(1),Norm(M))) { 489 cerr<<"swapRowsCols: i1,i2 = "<<i1<<','<<i2<<endl; 490 cerr<<"M0 = "<<TMV_Text(*this)<<" "<<M0<<endl; 491 cerr<<"M = "<<M<<endl; 492 cerr<<"S = "<<*this<<endl; 493 abort(); 494 } 495 #endif 496 return *this; 497 } 498 } 499 500 template <class T, int A> permuteRowsCols(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)501 SymMatrixView<T,A>& SymMatrixView<T,A>::permuteRowsCols( 502 const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 503 { 504 const ptrdiff_t* pi = p+i1; 505 for(ptrdiff_t i=i1;i<i2;++i,++pi) { 506 TMVAssert(*pi < size()); 507 swapRowsCols(i,*pi); 508 } 509 return *this; 510 } 511 512 template <class T, int A> reversePermuteRowsCols(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)513 SymMatrixView<T,A>& SymMatrixView<T,A>::reversePermuteRowsCols( 514 const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 515 { 516 const ptrdiff_t* pi = p+i2; 517 for(ptrdiff_t i=i2;i>i1;) { 518 --i; --pi; 519 TMVAssert(*pi < size()); 520 swapRowsCols(i,*pi); 521 } 522 return *this; 523 } 524 525 526 // 527 // Norms 528 // 529 530 template <class T> det() const531 T GenSymMatrix<T>::det() const 532 { return DivHelper<T>::det(); } 533 534 template <class T> logDet(T * sign) const535 RT GenSymMatrix<T>::logDet(T* sign) const 536 { return DivHelper<T>::logDet(sign); } 537 538 template <class T> isSingular() const539 bool GenSymMatrix<T>::isSingular() const 540 { return DivHelper<T>::isSingular(); } 541 542 #ifdef INST_INT 543 template <> det() const544 int GenSymMatrix<int>::det() const 545 { return IntegerDet(*this); } 546 547 template <> det() const548 std::complex<int> GenSymMatrix<std::complex<int> >::det() const 549 { return IntegerDet(*this); } 550 551 template <> logDet(int *) const552 int GenSymMatrix<int>::logDet(int* ) const 553 { TMVAssert(TMV_FALSE); return 0; } 554 555 template <> logDet(std::complex<int> *) const556 int GenSymMatrix<std::complex<int> >::logDet(std::complex<int>* ) const 557 { TMVAssert(TMV_FALSE); return 0; } 558 559 template <> isSingular() const560 bool GenSymMatrix<int>::isSingular() const 561 { return det() == 0; } 562 563 template <> isSingular() const564 bool GenSymMatrix<std::complex<int> >::isSingular() const 565 { return det() == 0; } 566 #endif 567 568 template <class T> sumElements() const569 T GenSymMatrix<T>::sumElements() const 570 { 571 T sum = diag().sumElements(); 572 if (size() > 1) { 573 T temp = upperTri().offDiag().sumElements(); 574 if (issym()) { 575 sum += RT(2) * temp; 576 } else { 577 // temp + conj(temp) = 2*real(temp) 578 sum += RT(2) * TMV_REAL(temp); 579 } 580 } 581 return sum; 582 } 583 584 template <class T> DoSumAbsElements(const GenSymMatrix<T> & m)585 static RT DoSumAbsElements(const GenSymMatrix<T>& m) 586 { 587 RT sum = m.diag().sumAbsElements(); 588 if (m.size() > 1) 589 sum += RT(2) * m.upperTri().offDiag().sumAbsElements(); 590 return sum; 591 } 592 593 template <class T> DoSumAbs2Elements(const GenSymMatrix<T> & m)594 static RT DoSumAbs2Elements(const GenSymMatrix<T>& m) 595 { 596 RT sum = m.diag().sumAbs2Elements(); 597 if (m.size() > 1) 598 sum += RT(2) * m.upperTri().offDiag().sumAbs2Elements(); 599 return sum; 600 } 601 602 #ifdef INST_INT DoSumAbsElements(const GenSymMatrix<std::complex<int>> &)603 static int DoSumAbsElements(const GenSymMatrix<std::complex<int> >& ) 604 { TMVAssert(TMV_FALSE); return 0; } 605 #endif 606 607 template <class T> sumAbsElements() const608 RT GenSymMatrix<T>::sumAbsElements() const 609 { return DoSumAbsElements(*this); } 610 611 template <class T> sumAbs2Elements() const612 RT GenSymMatrix<T>::sumAbs2Elements() const 613 { return DoSumAbs2Elements(*this); } 614 615 template <class T> normSq(const RT scale) const616 RT GenSymMatrix<T>::normSq(const RT scale) const 617 { 618 RT sum = diag().normSq(scale); 619 if (size() > 1) sum += RT(2) * upperTri().offDiag().normSq(scale); 620 return sum; 621 } 622 623 template <class T> NonLapNorm1(const GenSymMatrix<T> & m)624 static RT NonLapNorm1(const GenSymMatrix<T>& m) 625 { 626 RT max(0); 627 for(ptrdiff_t j=0;j<m.size();++j) { 628 RT temp = m.col(j,0,j).norm1(); 629 temp += m.col(j,j,m.size()).norm1(); 630 if (temp > max) max = temp; 631 } 632 return max; 633 } 634 635 template <class T> NonLapNormF(const GenSymMatrix<T> & m)636 static RT NonLapNormF(const GenSymMatrix<T>& m) 637 { 638 const RT eps = TMV_Epsilon<T>(); 639 640 RT mmax = m.maxAbs2Element(); 641 if (mmax == RT(0)) return RT(0); 642 else if (TMV_Underflow(mmax * mmax)) { 643 // Then we need to rescale, since underflow has caused 644 // rounding errors. 645 // Epsilon is a pure power of 2, so no rounding errors from 646 // rescaling. 647 const RT inveps = RT(1)/eps; 648 RT scale = inveps; 649 mmax *= scale; 650 const RT eps2 = eps*eps; 651 while (mmax < eps2) { scale *= inveps; mmax *= inveps; } 652 return TMV_SQRT(m.normSq(scale))/scale; 653 } else if (RT(1) / mmax == RT(0)) { 654 // Then mmax is already inf, so no hope of making it more accurate. 655 return mmax; 656 } else if (RT(1) / (mmax*mmax) == RT(0)) { 657 // Then we have overflow, so we need to rescale: 658 const RT inveps = RT(1)/eps; 659 RT scale = eps; 660 mmax *= scale; 661 while (mmax > inveps) { scale *= eps; mmax *= eps; } 662 return TMV_SQRT(m.normSq(scale))/scale; 663 } else { 664 return TMV_SQRT(m.normSq()); 665 } 666 } 667 668 template <class T> NonLapMaxAbsElement(const GenSymMatrix<T> & m)669 static inline RT NonLapMaxAbsElement(const GenSymMatrix<T>& m) 670 { return m.upperTri().maxAbsElement(); } 671 672 template <class T> NonLapMaxAbs2Element(const GenSymMatrix<T> & m)673 static inline RT NonLapMaxAbs2Element(const GenSymMatrix<T>& m) 674 { return m.upperTri().maxAbs2Element(); } 675 676 #ifdef XLAP 677 template <class T> LapNorm(const char c,const GenSymMatrix<T> & m)678 static RT LapNorm(const char c, const GenSymMatrix<T>& m) 679 { 680 switch(c) { 681 case 'M' : return NonLapMaxAbsElement(m); 682 case '1' : return NonLapNorm1(m); 683 case 'F' : return NonLapNormF(m); 684 default : TMVAssert(TMV_FALSE); 685 } 686 return RT(0); 687 } 688 #ifdef INST_DOUBLE 689 template <> LapNorm(const char c,const GenSymMatrix<double> & m)690 double LapNorm(const char c, const GenSymMatrix<double>& m) 691 { 692 TMVAssert(m.iscm() || m.isrm()); 693 char cc = c; 694 int N = m.size(); 695 int lda = m.iscm() ? m.stepj() : m.stepi(); 696 #ifndef LAPNOWORK 697 int lwork = c=='1' ? N : 0; 698 AlignedArray<double> work(lwork); 699 VectorViewOf(work.get(),lwork).setZero(); 700 #endif 701 double norm = LAPNAME(dlansy) ( 702 LAPCM LAPV(cc), 703 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO), 704 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1); 705 return norm; 706 } 707 template <> LapNorm(const char c,const GenSymMatrix<std::complex<double>> & m)708 double LapNorm(const char c, const GenSymMatrix<std::complex<double> >& m) 709 { 710 TMVAssert(m.iscm() || m.isrm()); 711 char cc = c; 712 int N = m.size(); 713 int lda = m.iscm() ? m.stepj() : m.stepi(); 714 #ifndef LAPNOWORK 715 int lwork = c=='1' ? N : 0; 716 AlignedArray<double> work(lwork); 717 VectorViewOf(work.get(),lwork).setZero(); 718 #endif 719 double norm = m.isherm() ? 720 LAPNAME(zlanhe) ( 721 LAPCM LAPV(cc), 722 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO), 723 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1) : 724 LAPNAME(zlansy) ( 725 LAPCM LAPV(cc), 726 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO), 727 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1); 728 return norm; 729 } 730 #endif 731 #ifdef INST_FLOAT 732 template <> LapNorm(const char c,const GenSymMatrix<float> & m)733 float LapNorm(const char c, const GenSymMatrix<float>& m) 734 { 735 TMVAssert(m.iscm() || m.isrm()); 736 char cc = c; 737 int N = m.size(); 738 int lda = m.iscm() ? m.stepj() : m.stepi(); 739 #ifndef LAPNOWORK 740 int lwork = c=='1' ? N : 0; 741 AlignedArray<float> work(lwork); 742 VectorViewOf(work.get(),lwork).setZero(); 743 #endif 744 float norm = LAPNAME(slansy) ( 745 LAPCM LAPV(cc), 746 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO), 747 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1); 748 return norm; 749 } 750 template <> LapNorm(const char c,const GenSymMatrix<std::complex<float>> & m)751 float LapNorm(const char c, const GenSymMatrix<std::complex<float> >& m) 752 { 753 TMVAssert(m.iscm() || m.isrm()); 754 char cc = c; 755 int N = m.size(); 756 int lda = m.iscm() ? m.stepj() : m.stepi(); 757 #ifndef LAPNOWORK 758 int lwork = c=='1' ? N : 0; 759 AlignedArray<float> work(lwork); 760 VectorViewOf(work.get(),lwork).setZero(); 761 #endif 762 float norm = m.isherm() ? 763 LAPNAME(clanhe) ( 764 LAPCM LAPV(cc), 765 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO), 766 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1) : 767 LAPNAME(clansy) ( 768 LAPCM LAPV(cc), 769 (m.iscm() == (m.uplo()==Upper) ? LAPCH_UP : LAPCH_LO), 770 LAPV(N),LAPP(m.cptr()),LAPV(lda) LAPWK(work.get()) LAP1 LAP1); 771 return norm; 772 } 773 #endif 774 #endif // XLAP 775 776 #ifdef INST_INT NonLapNormF(const GenSymMatrix<int> &)777 static int NonLapNormF(const GenSymMatrix<int>& ) 778 { TMVAssert(TMV_FALSE); return 0; } NonLapNormF(const GenSymMatrix<std::complex<int>> &)779 static int NonLapNormF(const GenSymMatrix<std::complex<int> >& ) 780 { TMVAssert(TMV_FALSE); return 0; } NonLapNorm1(const GenSymMatrix<std::complex<int>> &)781 static int NonLapNorm1(const GenSymMatrix<std::complex<int> >& ) 782 { TMVAssert(TMV_FALSE); return 0; } NonLapMaxAbsElement(const GenSymMatrix<std::complex<int>> &)783 static int NonLapMaxAbsElement(const GenSymMatrix<std::complex<int> >& ) 784 { TMVAssert(TMV_FALSE); return 0; } 785 #endif 786 787 template <class T> maxAbsElement() const788 RT GenSymMatrix<T>::maxAbsElement() const 789 { 790 #ifdef XLAP 791 if ((isrm() && stepi()>0) || (iscm() && stepj()>0)) 792 return LapNorm('M',*this); 793 else 794 #endif 795 return NonLapMaxAbsElement(*this); 796 } 797 template <class T> maxAbs2Element() const798 RT GenSymMatrix<T>::maxAbs2Element() const 799 { 800 #ifdef XLAP 801 if (Traits<T>::iscomplex) return NonLapMaxAbs2Element(*this); 802 else if ((isrm() && stepi()>0) || (iscm() && stepj()>0)) 803 return LapNorm('M',*this); 804 else 805 #endif 806 return NonLapMaxAbs2Element(*this); 807 } 808 template <class T> norm1() const809 RT GenSymMatrix<T>::norm1() const 810 { 811 #ifdef XLAP 812 if ((isrm() && stepi()>0) || (iscm() && stepj()>0)) 813 return LapNorm('1',*this); 814 else 815 #endif 816 return NonLapNorm1(*this); 817 } 818 template <class T> normF() const819 RT GenSymMatrix<T>::normF() const 820 { 821 #ifdef XLAP 822 if ((isrm() && stepi()>0) || (iscm() && stepj()>0)) 823 return LapNorm('F',*this); 824 else 825 #endif 826 return NonLapNormF(*this); 827 } 828 829 template <class T> DoNorm2(const GenSymMatrix<T> & m)830 static RT DoNorm2(const GenSymMatrix<T>& m) 831 { 832 if (m.size() == 0) return RT(0); 833 DiagMatrix<RT> S(m.size()); 834 if (m.isherm()) { 835 HermMatrix<T> m2(m); 836 SV_Decompose(m2.view(),S.view()); 837 } else { 838 SymMatrix<T> m2(m); 839 SV_Decompose(m2.view(),S.view()); 840 } 841 return S(0); 842 } 843 844 template <class T> DoCondition(const GenSymMatrix<T> & m)845 static RT DoCondition(const GenSymMatrix<T>& m) 846 { 847 if (m.size() == 0) return RT(1); 848 DiagMatrix<RT> S(m.size()); 849 if (m.isherm()) { 850 HermMatrix<T> m2(m); 851 SV_Decompose(m2.view(),S.view()); 852 } else { 853 SymMatrix<T> m2(m); 854 SV_Decompose(m2.view(),S.view()); 855 } 856 return TMV_ABS(S(0)/S(S.size()-1)); 857 } 858 859 #ifdef INST_INT DoNorm2(const GenSymMatrix<int> &)860 static int DoNorm2(const GenSymMatrix<int>& ) 861 { TMVAssert(TMV_FALSE); return 0; } DoCondition(const GenSymMatrix<int> &)862 static int DoCondition(const GenSymMatrix<int>& ) 863 { TMVAssert(TMV_FALSE); return 0; } DoNorm2(const GenSymMatrix<std::complex<int>> &)864 static int DoNorm2(const GenSymMatrix<std::complex<int> >& ) 865 { TMVAssert(TMV_FALSE); return 0; } DoCondition(const GenSymMatrix<std::complex<int>> &)866 static int DoCondition(const GenSymMatrix<std::complex<int> >& ) 867 { TMVAssert(TMV_FALSE); return 0; } 868 #endif 869 870 template <class T> doNorm2() const871 RT GenSymMatrix<T>::doNorm2() const 872 { return tmv::DoNorm2(*this); } 873 template <class T> doCondition() const874 RT GenSymMatrix<T>::doCondition() const 875 { return tmv::DoCondition(*this); } 876 877 878 // 879 // I/O 880 // 881 882 template <class T> write(const TMV_Writer & writer) const883 void GenSymMatrix<T>::write(const TMV_Writer& writer) const 884 { 885 const ptrdiff_t N = size(); 886 writer.begin(); 887 writer.writeCode(issym() ? "S" : "H"); 888 writer.writeSize(N); 889 writer.writeSimpleSize(N); 890 writer.writeStart(); 891 for(ptrdiff_t i=0;i<N;++i) { 892 writer.writeLParen(); 893 for(ptrdiff_t j=0;j<i+1;++j) { 894 if (j > 0) writer.writeSpace(); 895 writer.writeValue(cref(i,j)); 896 } 897 if (!writer.isCompact()) { 898 for(ptrdiff_t j=i+1;j<N;++j) { 899 writer.writeSpace(); 900 writer.writeValue(cref(i,j)); 901 } 902 } 903 writer.writeRParen(); 904 if (i < N-1) writer.writeRowEnd(); 905 } 906 writer.writeFinal(); 907 writer.end(); 908 } 909 910 #ifndef NOTHROW 911 template <class T> 912 class SymMatrixReadError : public ReadError 913 { 914 public : 915 SymMatrix<T> m; 916 ptrdiff_t i,j; 917 std::string exp,got; 918 ptrdiff_t s; 919 T v1, v2; 920 bool is, iseof, isbad; 921 SymMatrixReadError(std::istream & _is)922 SymMatrixReadError(std::istream& _is) throw() : 923 ReadError("SymMatrix."), 924 i(0), j(0), exp(0), got(0), s(0), v1(0), v2(0), 925 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} SymMatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)926 SymMatrixReadError( 927 std::istream& _is, 928 const std::string& _e, const std::string& _g) throw() : 929 ReadError("SymMatrix."), 930 i(0), j(0), exp(_e), got(_g), s(0), v1(0), v2(0), 931 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 932 SymMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)933 SymMatrixReadError( 934 ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m, 935 std::istream& _is, 936 const std::string& _e, const std::string& _g) throw() : 937 ReadError("SymMatrix."), 938 m(_m), i(_i), j(_j), exp(_e), got(_g), s(m.size()), v1(0), v2(0), 939 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} SymMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,T _v1=0,T _v2=0)940 SymMatrixReadError( 941 ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m, 942 std::istream& _is, T _v1=0, T _v2=0) throw() : 943 ReadError("SymMatrix."), 944 m(_m), i(_i), j(_j), s(m.size()), v1(_v1), v2(_v2), 945 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} SymMatrixReadError(const GenSymMatrix<T> & _m,std::istream & _is,ptrdiff_t _s)946 SymMatrixReadError( 947 const GenSymMatrix<T>& _m, std::istream& _is, ptrdiff_t _s) throw() : 948 ReadError("SymMatrix."), 949 m(_m), i(0), j(0), exp(0), got(0), s(_s), v1(0), v2(0), 950 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 951 SymMatrixReadError(const SymMatrixReadError<T> & rhs)952 SymMatrixReadError(const SymMatrixReadError<T>& rhs) : 953 m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got), 954 s(rhs.s), v1(rhs.v1), v2(rhs.v2), 955 is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {} ~SymMatrixReadError()956 virtual ~SymMatrixReadError() throw() {} 957 write(std::ostream & os) const958 virtual void write(std::ostream& os) const throw() 959 { 960 os<<"TMV Read Error: Reading istream input for SymMatrix\n"; 961 if (exp != got) { 962 os<<"Wrong format: expected '"<<exp<<"'"; 963 if (isReal(T()) && exp == "S") os<<" (or 'H')"; 964 os<<", got '"<<got<<"'.\n"; 965 } 966 if (s != m.size()) { 967 os<<"Wrong size: expected "<<m.size()<<", got "<<s<<".\n"; 968 } 969 if (!is) { 970 if (iseof) { 971 os<<"Input stream reached end-of-file prematurely.\n"; 972 } else if (isbad) { 973 os<<"Input stream is corrupted.\n"; 974 } else { 975 os<<"Input stream cannot read next character.\n"; 976 } 977 } 978 if (v1 != v2) { 979 os<<"Input matrix is not symmetric.\n"; 980 os<<"Lower triangle has the value "<<v1<<" at ("<<i<<","<<j<<")\n"; 981 os<<"Upper triangle has the value "<<v2<<" at ("<<j<<","<<i<<")\n"; 982 } 983 if (m.size() > 0) { 984 os<<"The portion of the SymMatrix which was successfully " 985 "read is: \n"; 986 const ptrdiff_t N = m.size(); 987 for(ptrdiff_t ii=0;ii<i;++ii) { 988 os<<"( "; 989 for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' '; 990 os<<" )\n"; 991 } 992 os<<"( "; 993 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' '; 994 os<<" )\n"; 995 } 996 } 997 }; 998 999 template <class T> 1000 class HermMatrixReadError : public ReadError 1001 { 1002 public : 1003 HermMatrix<T> m; 1004 ptrdiff_t i,j; 1005 std::string exp,got; 1006 ptrdiff_t s; 1007 T v1, v2; 1008 bool is, iseof, isbad; 1009 HermMatrixReadError(std::istream & _is)1010 HermMatrixReadError(std::istream& _is) throw() : 1011 ReadError("HermMatrix."), 1012 i(0), j(0), exp(0), got(0), s(0), v1(0), v2(0), 1013 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} HermMatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)1014 HermMatrixReadError( 1015 std::istream& _is, 1016 const std::string& _e, const std::string& _g) throw() : 1017 ReadError("HermMatrix."), 1018 i(0), j(0), exp(_e), got(_g), s(0), v1(0), v2(0), 1019 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 1020 HermMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)1021 HermMatrixReadError( 1022 ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m, 1023 std::istream& _is, 1024 const std::string& _e, const std::string& _g) throw() : 1025 ReadError("HermMatrix."), 1026 m(_m), i(_i), j(_j), exp(_e), got(_g), 1027 s(m.size()), v1(0), v2(0), 1028 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} HermMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenSymMatrix<T> & _m,std::istream & _is,T _v1=0,T _v2=0)1029 HermMatrixReadError( 1030 ptrdiff_t _i, ptrdiff_t _j, const GenSymMatrix<T>& _m, 1031 std::istream& _is, T _v1=0, T _v2=0) throw() : 1032 ReadError("HermMatrix."), 1033 m(_m), i(_i), j(_j), exp(0), got(0), 1034 s(m.size()), v1(_v1), v2(_v2), 1035 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} HermMatrixReadError(const GenSymMatrix<T> & _m,std::istream & _is,ptrdiff_t _s)1036 HermMatrixReadError( 1037 const GenSymMatrix<T>& _m, std::istream& _is, ptrdiff_t _s) throw() : 1038 ReadError("HermMatrix."), 1039 m(_m), i(0), j(0), exp(0), got(0), s(_s), v1(0), v2(0), 1040 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 1041 HermMatrixReadError(const HermMatrixReadError<T> & rhs)1042 HermMatrixReadError(const HermMatrixReadError<T>& rhs) : 1043 m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got), 1044 s(rhs.s), v1(rhs.v1), v2(rhs.v2), 1045 is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {} ~HermMatrixReadError()1046 virtual ~HermMatrixReadError() throw() {} 1047 write(std::ostream & os) const1048 virtual void write(std::ostream& os) const throw() 1049 { 1050 os<<"TMV Read Error: Reading istream input for HermMatrix\n"; 1051 if (exp != got) { 1052 os<<"Wrong format: expected '"<<exp<<"'"; 1053 if (isReal(T()) && exp == "H") os<<" (or 'S')"; 1054 os<<", got '"<<got<<"'.\n"; 1055 } 1056 if (s != m.size()) { 1057 os<<"Wrong size: expected "<<m.size()<<", got "<<s<<".\n"; 1058 } 1059 if (!is) { 1060 if (iseof) { 1061 os<<"Input stream reached end-of-file prematurely.\n"; 1062 } else if (isbad) { 1063 os<<"Input stream is corrupted.\n"; 1064 } else { 1065 os<<"Input stream cannot read next character.\n"; 1066 } 1067 } 1068 if (i==j && v1 != T(0)) { 1069 os<<"Non-real value found on diagonal: "<<v1<<std::endl; 1070 } 1071 if (i!=j && v1 != v2) { 1072 os<<"Input matrix is not Hermitian.\n"; 1073 os<<"Lower triangle has the value "<<v1<<" at ("<<i<<","<<j<<")\n"; 1074 os<<"Upper triangle has the value "<<v2<<" at ("<<j<<","<<i<<")\n"; 1075 } 1076 if (m.size() > 0) { 1077 os<<"The portion of the HermMatrix which was successfully " 1078 "read is: \n"; 1079 const ptrdiff_t N = m.size(); 1080 for(ptrdiff_t ii=0;ii<i;++ii) { 1081 os<<"( "; 1082 for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' '; 1083 os<<" )\n"; 1084 } 1085 os<<"( "; 1086 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' '; 1087 os<<" )\n"; 1088 } 1089 } 1090 }; 1091 #endif 1092 1093 template <class T> FinishRead(const TMV_Reader & reader,SymMatrixView<T> m)1094 static void FinishRead(const TMV_Reader& reader, SymMatrixView<T> m) 1095 { 1096 const ptrdiff_t N = m.size(); 1097 std::string exp, got; 1098 T temp; 1099 if (!reader.readStart(exp,got)) { 1100 #ifdef NOTHROW 1101 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1102 exit(1); 1103 #else 1104 if (m.issym()) 1105 throw SymMatrixReadError<T>(0,0,m,reader.getis(),exp,got); 1106 else 1107 throw HermMatrixReadError<T>(0,0,m,reader.getis(),exp,got); 1108 #endif 1109 } 1110 for(ptrdiff_t i=0;i<N;++i) { 1111 if (!reader.readLParen(exp,got)) { 1112 #ifdef NOTHROW 1113 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1114 exit(1); 1115 #else 1116 if (m.issym()) 1117 throw SymMatrixReadError<T>(i,0,m,reader.getis(),exp,got); 1118 else 1119 throw HermMatrixReadError<T>(i,0,m,reader.getis(),exp,got); 1120 #endif 1121 } 1122 for(ptrdiff_t j=0;j<i+1;++j) { 1123 if (j>0 && !reader.readSpace(exp,got)) { 1124 #ifdef NOTHROW 1125 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1126 exit(1); 1127 #else 1128 if (m.issym()) 1129 throw SymMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1130 else 1131 throw HermMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1132 #endif 1133 } 1134 if (!reader.readValue(temp)) { 1135 #ifdef NOTHROW 1136 std::cerr<<"SymMatrix Read Error: reading value\n"; 1137 exit(1); 1138 #else 1139 if (m.issym()) 1140 throw SymMatrixReadError<T>(i,j,m,reader.getis()); 1141 else 1142 throw HermMatrixReadError<T>(i,j,m,reader.getis()); 1143 #endif 1144 } 1145 if (j == i && !m.issym() && TMV_IMAG(temp) != RT(0)) { 1146 #ifdef NOTHROW 1147 std::cerr<<"SymMatrix Read Error: "<<temp<<" not real\n"; 1148 exit(1); 1149 #else 1150 throw HermMatrixReadError<T>(i,j,m,reader.getis(),temp); 1151 #endif 1152 } 1153 if (j < i && !reader.isCompact()) { 1154 T temp2 = m.issym() ? m.cref(j,i) : TMV_CONJ(m.cref(j,i)); 1155 if (temp != temp2) { 1156 #ifdef NOTHROW 1157 std::cerr<<"SymMatrix Read Error: "<<temp<<" != "<<temp2<<"\n"; 1158 exit(1); 1159 #else 1160 if (m.issym()) 1161 throw SymMatrixReadError<T>(i,j,m,reader.getis(),temp,m.cref(j,i)); 1162 else 1163 throw HermMatrixReadError<T>(i,j,m,reader.getis(),temp,m.cref(j,i)); 1164 #endif 1165 } 1166 } else { 1167 m.ref(i,j) = temp; 1168 } 1169 } 1170 if (!reader.isCompact()) { 1171 for(ptrdiff_t j=i+1;j<N;++j) { 1172 if (!reader.readSpace(exp,got)) { 1173 #ifdef NOTHROW 1174 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1175 exit(1); 1176 #else 1177 if (m.issym()) 1178 throw SymMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1179 else 1180 throw HermMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1181 #endif 1182 } 1183 if (!reader.readValue(temp)) { 1184 #ifdef NOTHROW 1185 std::cerr<<"SymMatrix Read Error: reading value\n"; 1186 exit(1); 1187 #else 1188 if (m.issym()) 1189 throw SymMatrixReadError<T>(i,j,m,reader.getis()); 1190 else 1191 throw HermMatrixReadError<T>(i,j,m,reader.getis()); 1192 #endif 1193 } 1194 m.ref(i,j) = temp; 1195 } 1196 } 1197 if (!reader.readRParen(exp,got)) { 1198 #ifdef NOTHROW 1199 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1200 exit(1); 1201 #else 1202 if (m.issym()) 1203 throw SymMatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1204 else 1205 throw HermMatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1206 #endif 1207 } 1208 if (i < N-1 && !reader.readRowEnd(exp,got)) { 1209 #ifdef NOTHROW 1210 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1211 exit(1); 1212 #else 1213 if (m.issym()) 1214 throw SymMatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1215 else 1216 throw HermMatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1217 #endif 1218 } 1219 } 1220 if (!reader.readFinal(exp,got)) { 1221 #ifdef NOTHROW 1222 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1223 exit(1); 1224 #else 1225 if (m.issym()) 1226 throw SymMatrixReadError<T>(N,0,m,reader.getis(),exp,got); 1227 else 1228 throw HermMatrixReadError<T>(N,0,m,reader.getis(),exp,got); 1229 #endif 1230 } 1231 } 1232 1233 template <class T, int A> read(const TMV_Reader & reader)1234 void SymMatrix<T,A>::read(const TMV_Reader& reader) 1235 { 1236 std::string exp,got; 1237 bool ok = isReal(T()) ? 1238 reader.readCode("S","H",exp,got) : reader.readCode("S",exp,got); 1239 if (!ok) { 1240 #ifdef NOTHROW 1241 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1242 exit(1); 1243 #else 1244 throw SymMatrixReadError<T>(reader.getis(),exp,got); 1245 #endif 1246 } 1247 ptrdiff_t s=size(); 1248 if (!reader.readSize(s,exp,got)) { 1249 #ifdef NOTHROW 1250 std::cerr<<"SymMatrix Read Error: reading size\n"; 1251 exit(1); 1252 #else 1253 throw SymMatrixReadError<T>(reader.getis(),exp,got); 1254 #endif 1255 } 1256 if (s != size()) resize(s); 1257 s=size(); 1258 if (!reader.readSimpleSize(s,exp,got)) { 1259 #ifdef NOTHROW 1260 std::cerr<<"SymMatrix Read Error: reading size\n"; 1261 exit(1); 1262 #else 1263 throw SymMatrixReadError<T>(reader.getis(),exp,got); 1264 #endif 1265 } 1266 if (s != size()) { 1267 #ifdef NOTHROW 1268 std::cerr<<"SymMatrix Read Error: Wrong size\n"; 1269 exit(1); 1270 #else 1271 throw SymMatrixReadError<T>(*this,reader.getis(),s); 1272 #endif 1273 } 1274 SymMatrixView<T> v = view(); 1275 FinishRead(reader,v); 1276 } 1277 1278 template <class T, int A> read(const TMV_Reader & reader)1279 void HermMatrix<T,A>::read(const TMV_Reader& reader) 1280 { 1281 std::string exp,got; 1282 bool ok = isReal(T()) ? 1283 reader.readCode("S","H",exp,got) : reader.readCode("H",exp,got); 1284 if (!ok) { 1285 #ifdef NOTHROW 1286 std::cerr<<"HermMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1287 exit(1); 1288 #else 1289 throw HermMatrixReadError<T>(reader.getis(),exp,got); 1290 #endif 1291 } 1292 ptrdiff_t s=size(); 1293 if (!reader.readSize(s,exp,got)) { 1294 #ifdef NOTHROW 1295 std::cerr<<"HermMatrix Read Error: reading size\n"; 1296 exit(1); 1297 #else 1298 throw HermMatrixReadError<T>(reader.getis(),exp,got); 1299 #endif 1300 } 1301 if (s != size()) resize(s); 1302 s=size(); 1303 if (!reader.readSimpleSize(s,exp,got)) { 1304 #ifdef NOTHROW 1305 std::cerr<<"HermMatrix Read Error: reading size\n"; 1306 exit(1); 1307 #else 1308 throw HermMatrixReadError<T>(reader.getis(),exp,got); 1309 #endif 1310 } 1311 if (s != size()) { 1312 #ifdef NOTHROW 1313 std::cerr<<"HermMatrix Read Error: Wrong size\n"; 1314 exit(1); 1315 #else 1316 throw HermMatrixReadError<T>(*this,reader.getis(),s); 1317 #endif 1318 } 1319 SymMatrixView<T> v = view(); 1320 FinishRead(reader,v); 1321 } 1322 1323 template <class T, int A> read(const TMV_Reader & reader)1324 void SymMatrixView<T,A>::read(const TMV_Reader& reader) 1325 { 1326 std::string exp,got; 1327 bool ok = isReal(T()) ? 1328 reader.readCode("S","H",exp,got) : 1329 reader.readCode(issym()?"S":"H",exp,got); 1330 if (!ok) { 1331 #ifdef NOTHROW 1332 std::cerr<<"SymMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1333 exit(1); 1334 #else 1335 if (issym()) 1336 throw SymMatrixReadError<T>(reader.getis(),exp,got); 1337 else 1338 throw HermMatrixReadError<T>(reader.getis(),exp,got); 1339 #endif 1340 } 1341 ptrdiff_t s=size(); 1342 if (!reader.readSize(s,exp,got)) { 1343 #ifdef NOTHROW 1344 std::cerr<<"SymMatrix Read Error: reading size\n"; 1345 exit(1); 1346 #else 1347 if (issym()) 1348 throw SymMatrixReadError<T>(reader.getis(),exp,got); 1349 else 1350 throw HermMatrixReadError<T>(reader.getis(),exp,got); 1351 #endif 1352 } 1353 if (s != size()) { 1354 #ifdef NOTHROW 1355 std::cerr<<"SymMatrix Read Error: Wrong size\n"; 1356 exit(1); 1357 #else 1358 if (issym()) 1359 throw SymMatrixReadError<T>(*this,reader.getis(),s); 1360 else 1361 throw HermMatrixReadError<T>(*this,reader.getis(),s); 1362 #endif 1363 } 1364 s=size(); 1365 if (!reader.readSimpleSize(s,exp,got)) { 1366 #ifdef NOTHROW 1367 std::cerr<<"SymMatrix Read Error: reading size\n"; 1368 exit(1); 1369 #else 1370 if (issym()) 1371 throw SymMatrixReadError<T>(reader.getis(),exp,got); 1372 else 1373 throw HermMatrixReadError<T>(reader.getis(),exp,got); 1374 #endif 1375 } 1376 if (s != size()) { 1377 #ifdef NOTHROW 1378 std::cerr<<"SymMatrix Read Error: Wrong size\n"; 1379 exit(1); 1380 #else 1381 if (issym()) 1382 throw SymMatrixReadError<T>(*this,reader.getis(),s); 1383 else 1384 throw HermMatrixReadError<T>(*this,reader.getis(),s); 1385 #endif 1386 } 1387 FinishRead(reader,*this); 1388 } 1389 1390 1391 1392 #define InstFile "TMV_SymMatrix.inst" 1393 #include "TMV_Inst.h" 1394 #undef InstFile 1395 1396 } // namespace tmv 1397 1398 1399