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 #include "TMV_Blas.h" 25 #include "tmv/TMV_Matrix.h" 26 #include "tmv/TMV_Vector.h" 27 #include "tmv/TMV_DiagMatrix.h" 28 #include "tmv/TMV_Divider.h" 29 #include "tmv/TMV_LUD.h" 30 #include "tmv/TMV_QRD.h" 31 #include "tmv/TMV_QRPD.h" 32 #include "tmv/TMV_SVD.h" 33 #include "tmv/TMV_MatrixArith.h" 34 #include "tmv/TMV_VIt.h" 35 #include "TMV_IntegerDet.h" 36 #include <iostream> 37 38 namespace tmv { 39 40 #define RT TMV_RealType(T) 41 42 #ifdef TMV_BLOCKSIZE 43 #define PERM_BLOCKSIZE TMV_BLOCKSIZE/2 44 #else 45 #define PERM_BLOCKSIZE 32 46 #endif 47 48 // 49 // Access 50 // 51 52 template <class T> cref(ptrdiff_t i,ptrdiff_t j) const53 T GenMatrix<T>::cref(ptrdiff_t i, ptrdiff_t j) const 54 { 55 const T* mi = cptr() + i*stepi() + j*stepj(); 56 return isconj() ? TMV_CONJ(*mi) : *mi; 57 } 58 59 template <class T, int A> ref(ptrdiff_t i,ptrdiff_t j)60 typename MatrixView<T,A>::reference MatrixView<T,A>::ref(ptrdiff_t i, ptrdiff_t j) 61 { 62 T* mi = ptr() + i*itssi + j*stepj(); 63 return RefHelper<T>::makeRef(mi,ct()); 64 } 65 66 template <class T> setDiv() const67 void GenMatrix<T>::setDiv() const 68 { 69 if (!this->divIsSet()) { 70 DivType dt = this->getDivType(); 71 TMVAssert(dt == tmv::LU || dt == tmv::QR || 72 dt == tmv::QRP || dt == tmv::SV); 73 switch (dt) { 74 case LU : 75 this->divider.reset( 76 new LUDiv<T>(*this,this->divIsInPlace())); 77 break; 78 case QR : 79 this->divider.reset( 80 new QRDiv<T>(*this,this->divIsInPlace())); 81 break; 82 case QRP : 83 this->divider.reset( 84 new QRPDiv<T>(*this,this->divIsInPlace())); 85 break; 86 case SV : 87 this->divider.reset( 88 new SVDiv<T>(*this,this->divIsInPlace())); 89 break; 90 default : 91 // The above assert should have already failed 92 // so go ahead and fall through. 93 break; 94 } 95 } 96 } 97 98 #ifdef INST_INT 99 template <> setDiv() const100 void GenMatrix<int>::setDiv() const 101 { TMVAssert(TMV_FALSE); } 102 template <> setDiv() const103 void GenMatrix<std::complex<int> >::setDiv() const 104 { TMVAssert(TMV_FALSE); } 105 #endif 106 107 // Note: These need to be in the .cpp file, not the .h file for 108 // dynamic libraries. Apparently the typeinfo used for dynamic_cast 109 // doesn't get shared correctly by different modules, so the 110 // dynamic_cast fails when called in one module for an object that 111 // was created in a different module. 112 // 113 // So putting these functions here puts the dynamic cast in the shared 114 // library, which is also where it is created (by setDiv() above). 115 template <class T> divIsLUDiv() const116 bool GenMatrix<T>::divIsLUDiv() const 117 { return dynamic_cast<const LUDiv<T>*>(this->getDiv()); } 118 119 template <class T> divIsQRDiv() const120 bool GenMatrix<T>::divIsQRDiv() const 121 { return dynamic_cast<const QRDiv<T>*>(this->getDiv()); } 122 123 template <class T> divIsQRPDiv() const124 bool GenMatrix<T>::divIsQRPDiv() const 125 { return dynamic_cast<const QRPDiv<T>*>(this->getDiv()); } 126 127 template <class T> divIsSVDiv() const128 bool GenMatrix<T>::divIsSVDiv() const 129 { return dynamic_cast<const SVDiv<T>*>(this->getDiv()); } 130 131 132 #ifdef INST_INT 133 template <> divIsLUDiv() const134 bool GenMatrix<int>::divIsLUDiv() const 135 { return false; } 136 template <> divIsQRDiv() const137 bool GenMatrix<int>::divIsQRDiv() const 138 { return false; } 139 template <> divIsQRPDiv() const140 bool GenMatrix<int>::divIsQRPDiv() const 141 { return false; } 142 template <> divIsSVDiv() const143 bool GenMatrix<int>::divIsSVDiv() const 144 { return false; } 145 146 template <> divIsLUDiv() const147 bool GenMatrix<std::complex<int> >::divIsLUDiv() const 148 { return false; } 149 template <> divIsQRDiv() const150 bool GenMatrix<std::complex<int> >::divIsQRDiv() const 151 { return false; } 152 template <> divIsQRPDiv() const153 bool GenMatrix<std::complex<int> >::divIsQRPDiv() const 154 { return false; } 155 template <> divIsSVDiv() const156 bool GenMatrix<std::complex<int> >::divIsSVDiv() const 157 { return false; } 158 #endif 159 160 // 161 // OK? (SubMatrix, SubVector) 162 // 163 164 template <class T> hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const165 bool GenMatrix<T>::hasSubMatrix( 166 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 167 { 168 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 169 bool ok = true; 170 if (istep == 0) { 171 ok = false; 172 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 173 } 174 if (i1 < 0 || i1 >= colsize()) { 175 ok = false; 176 std::cerr<<"first col element ("<<i1<<") must be in 0 -- "; 177 std::cerr<<colsize()-1<<std::endl; 178 } 179 if (i2-istep < 0 || i2-istep >= colsize()) { 180 ok = false; 181 std::cerr<<"last col element ("<<i2-istep<<") must be in 0 -- "; 182 std::cerr<<colsize()-1<<std::endl; 183 } 184 if ((i2-i1)%istep != 0) { 185 ok = false; 186 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 187 std::cerr<<istep<<")\n"; 188 } 189 if ((i2-i1)/istep < 0) { 190 ok = false; 191 std::cerr<<"n col elements ("<<(i2-i1)/istep<<") must be nonnegative\n"; 192 } 193 if (jstep == 0) { 194 ok = false; 195 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 196 } 197 if (j1 < 0 || j1 >= rowsize()) { 198 ok = false; 199 std::cerr<<"first row element ("<<j1<<") must be in 0 -- "; 200 std::cerr<<rowsize()-1<<std::endl; 201 } 202 if (j2-jstep < 0 || j2-jstep >= rowsize()) { 203 ok = false; 204 std::cerr<<"last row element ("<<j2-jstep<<") must be in 0 -- "; 205 std::cerr<<rowsize()-1<<std::endl; 206 } 207 if ((j2-j1)%jstep != 0) { 208 ok = false; 209 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 210 std::cerr<<jstep<<")\n"; 211 } 212 if ((j2-j1)/jstep < 0) { 213 ok = false; 214 std::cerr<<"n row elements ("<<(j2-j1)/jstep<<") must be nonnegative\n"; 215 } 216 return ok; 217 } 218 219 template <class T> hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t size) const220 bool GenMatrix<T>::hasSubVector( 221 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 222 { 223 if (size==0) return true; 224 bool ok = true; 225 if (istep == 0 && jstep == 0) { 226 ok = false; 227 std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep; 228 std::cerr<<") can not both be 0\n"; 229 } 230 if (i < 0 || i >= colsize()) { 231 ok = false; 232 std::cerr<<"i ("<<i<<") must be in 0 -- "<<colsize()-1<<std::endl; 233 } 234 if (j < 0 || j >= rowsize()) { 235 ok = false; 236 std::cerr<<"j ("<<j<<") must be in 0 -- "<<rowsize()-1<<std::endl; 237 } 238 ptrdiff_t i2 = i+istep*(size-1); 239 ptrdiff_t j2 = j+jstep*(size-1); 240 if (i2 < 0 || i2 >= colsize()) { 241 ok = false; 242 std::cerr<<"last element's i ("<<i2<<") must be in 0 -- "; 243 std::cerr<<colsize()-1<<std::endl; 244 } 245 if (j2 < 0 || j2 >= rowsize()) { 246 ok = false; 247 std::cerr<<"last element's j ("<<j2<<") must be in 0 -- "; 248 std::cerr<<rowsize()-1<<std::endl; 249 } 250 return ok; 251 } 252 253 template <class T> hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const254 bool ConstMatrixView<T,FortranStyle>::hasSubMatrix( 255 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 256 { 257 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 258 bool ok = true; 259 if (istep == 0) { 260 ok = false; 261 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 262 } 263 if (i1 < 1 || i1 > this->colsize()) { 264 ok = false; 265 std::cerr<<"first col element ("<<i1<<") must be in 1 -- "; 266 std::cerr<<this->colsize()<<std::endl; 267 } 268 if (i2 < 1 || i2 > this->colsize()) { 269 ok = false; 270 std::cerr<<"last col element ("<<i2<<") must be in 1 -- "; 271 std::cerr<<this->colsize()<<std::endl; 272 } 273 if ((i2-i1)%istep != 0) { 274 ok = false; 275 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 276 std::cerr<<istep<<")\n"; 277 } 278 if ((i2-i1)/istep < 0) { 279 ok = false; 280 std::cerr<<"n col elements ("<<(i2-i1)/istep+1<<") must be positive\n"; 281 } 282 if (jstep == 0) { 283 ok = false; 284 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 285 } 286 if (j1 < 1 || j1 > this->rowsize()) { 287 ok = false; 288 std::cerr<<"first row element ("<<j1<<") must be in 1 -- "; 289 std::cerr<<this->rowsize()<<std::endl; 290 } 291 if (j2 < 1 || j2 > this->rowsize()) { 292 ok = false; 293 std::cerr<<"last row element ("<<j2<<") must be in 1 -- "; 294 std::cerr<<this->rowsize()<<std::endl; 295 } 296 if ((j2-j1)%jstep != 0) { 297 ok = false; 298 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 299 std::cerr<<jstep<<")\n"; 300 } 301 if ((j2-j1)/jstep < 0) { 302 ok = false; 303 std::cerr<<"n row elements ("<<(j2-j1)/jstep+1<<") must be positive\n"; 304 } 305 return ok; 306 } 307 308 template <class T> hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t size) const309 bool ConstMatrixView<T,FortranStyle>::hasSubVector( 310 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 311 { 312 if (size==0) return true; 313 bool ok = true; 314 if (istep == 0 && jstep == 0) { 315 ok = false; 316 std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep; 317 std::cerr<<") can not both be 0\n"; 318 } 319 if (i < 1 || i > this->colsize()) { 320 ok = false; 321 std::cerr<<"i ("<<i<<") must be in 1 -- "<<this->colsize()<<std::endl; 322 } 323 if (j < 1 || j > this->rowsize()) { 324 ok = false; 325 std::cerr<<"j ("<<j<<") must be in 1 -- "<<this->rowsize()<<std::endl; 326 } 327 ptrdiff_t i2 = i+istep*(size-1); 328 ptrdiff_t j2 = j+jstep*(size-1); 329 if (i2 < 1 || i2 > this->colsize()) { 330 ok = false; 331 std::cerr<<"last element's i ("<<i2<<") must be in 1 -- "; 332 std::cerr<<this->colsize()<<std::endl; 333 } 334 if (j2 < 1 || j2 > this->rowsize()) { 335 ok = false; 336 std::cerr<<"last element's j ("<<j2<<") must be in 1 -- "; 337 std::cerr<<this->rowsize()<<std::endl; 338 } 339 return ok; 340 } 341 342 343 // 344 // Norms, det, etc. 345 // 346 347 template <class T> det() const348 T GenMatrix<T>::det() const 349 { return DivHelper<T>::det(); } 350 351 template <class T> logDet(T * sign) const352 RT GenMatrix<T>::logDet(T* sign) const 353 { return DivHelper<T>::logDet(sign); } 354 355 template <class T> isSingular() const356 bool GenMatrix<T>::isSingular() const 357 { return DivHelper<T>::isSingular(); } 358 359 #ifdef INST_INT 360 template <> det() const361 int GenMatrix<int>::det() const 362 { return IntegerDet(*this); } 363 364 template <> det() const365 std::complex<int> GenMatrix<std::complex<int> >::det() const 366 { return IntegerDet(*this); } 367 368 template <> logDet(int *) const369 int GenMatrix<int>::logDet(int* ) const 370 { TMVAssert(TMV_FALSE); return 0; } 371 372 template <> logDet(std::complex<int> *) const373 int GenMatrix<std::complex<int> >::logDet(std::complex<int>* ) const 374 { TMVAssert(TMV_FALSE); return 0; } 375 376 template <> isSingular() const377 bool GenMatrix<int>::isSingular() const 378 { return det() == 0; } 379 380 template <> isSingular() const381 bool GenMatrix<std::complex<int> >::isSingular() const 382 { return det() == 0; } 383 #endif 384 385 template <class T> sumElements() const386 T GenMatrix<T>::sumElements() const 387 { 388 if (canLinearize()) return constLinearView().sumElements(); 389 else { 390 T sum(0); 391 if (iscm()) { 392 const ptrdiff_t N = rowsize(); 393 for(ptrdiff_t j=0;j<N;++j) sum += col(j).sumElements(); 394 } else { 395 const ptrdiff_t M = colsize(); 396 for(ptrdiff_t i=0;i<M;++i) sum += row(i).sumElements(); 397 } 398 return sum; 399 } 400 } 401 402 template <class T> DoSumAbsElements(const GenMatrix<T> & m)403 static RT DoSumAbsElements(const GenMatrix<T>& m) 404 { 405 if (m.canLinearize()) return m.constLinearView().sumAbsElements(); 406 else { 407 RT sum(0); 408 if (m.iscm()) { 409 const ptrdiff_t N = m.rowsize(); 410 for(ptrdiff_t j=0;j<N;++j) sum += m.col(j).sumAbsElements(); 411 } else { 412 const ptrdiff_t M = m.colsize(); 413 for(ptrdiff_t i=0;i<M;++i) sum += m.row(i).sumAbsElements(); 414 } 415 return sum; 416 } 417 } 418 419 template <class T> DoSumAbs2Elements(const GenMatrix<T> & m)420 static RT DoSumAbs2Elements(const GenMatrix<T>& m) 421 { 422 if (m.canLinearize()) return m.constLinearView().sumAbs2Elements(); 423 else { 424 RT sum(0); 425 if (m.iscm()) { 426 const ptrdiff_t N = m.rowsize(); 427 for(ptrdiff_t j=0;j<N;++j) sum += m.col(j).sumAbs2Elements(); 428 } else { 429 const ptrdiff_t M = m.colsize(); 430 for(ptrdiff_t i=0;i<M;++i) sum += m.row(i).sumAbs2Elements(); 431 } 432 return sum; 433 } 434 } 435 436 #ifdef INST_INT DoSumAbsElements(const GenMatrix<std::complex<int>> &)437 static int DoSumAbsElements(const GenMatrix<std::complex<int> >& ) 438 { TMVAssert(TMV_FALSE); return 0; } 439 #endif 440 441 template <class T> sumAbsElements() const442 RT GenMatrix<T>::sumAbsElements() const 443 { return DoSumAbsElements(*this); } 444 445 template <class T> sumAbs2Elements() const446 RT GenMatrix<T>::sumAbs2Elements() const 447 { return DoSumAbs2Elements(*this); } 448 449 template <class T> normSq(const RT scale) const450 RT GenMatrix<T>::normSq(const RT scale) const 451 { 452 if (canLinearize()) return constLinearView().normSq(scale); 453 else { 454 RT sum(0); 455 if (isrm()) { 456 const ptrdiff_t M = colsize(); 457 for(ptrdiff_t i=0;i<M;++i) sum += row(i).normSq(scale); 458 } else { 459 const ptrdiff_t N = rowsize(); 460 for(ptrdiff_t j=0;j<N;++j) sum += col(j).normSq(scale); 461 } 462 return sum; 463 } 464 } 465 466 template <class T> NonLapMaxAbsElement(const GenMatrix<T> & m)467 static RT NonLapMaxAbsElement(const GenMatrix<T>& m) 468 { 469 if (m.canLinearize()) return m.constLinearView().maxAbsElement(); 470 else { 471 RT max(0); 472 if (m.iscm()) { 473 const ptrdiff_t N = m.rowsize(); 474 for(ptrdiff_t j=0;j<N;++j) { 475 RT temp = m.col(j).normInf(); 476 if (temp > max) max = temp; 477 } 478 } else { 479 const ptrdiff_t M = m.colsize(); 480 for(ptrdiff_t i=0;i<M;++i) { 481 RT temp = m.row(i).normInf(); 482 if (temp > max) max = temp; 483 } 484 } 485 return max; 486 } 487 } 488 489 template <class T> NonLapMaxAbs2Element(const GenMatrix<T> & m)490 static RT NonLapMaxAbs2Element(const GenMatrix<T>& m) 491 { 492 if (m.canLinearize()) return m.constLinearView().maxAbs2Element(); 493 else { 494 RT max(0); 495 if (m.iscm()) { 496 const ptrdiff_t N = m.rowsize(); 497 for(ptrdiff_t j=0;j<N;++j) { 498 RT temp = m.col(j).maxAbs2Element(); 499 if (temp > max) max = temp; 500 } 501 } else { 502 const ptrdiff_t M = m.colsize(); 503 for(ptrdiff_t i=0;i<M;++i) { 504 RT temp = m.row(i).maxAbs2Element(); 505 if (temp > max) max = temp; 506 } 507 } 508 return max; 509 } 510 } 511 512 template <class T> NonLapNorm1(const GenMatrix<T> & m)513 static RT NonLapNorm1(const GenMatrix<T>& m) 514 { 515 RT max(0); 516 const ptrdiff_t N = m.rowsize(); 517 for(ptrdiff_t j=0;j<N;++j) { 518 RT temp = m.col(j).norm1(); 519 if (temp > max) max = temp; 520 } 521 return max; 522 } 523 524 template <class T> NonLapNormF(const GenMatrix<T> & m)525 static RT NonLapNormF(const GenMatrix<T>& m) 526 { 527 const RT eps = TMV_Epsilon<T>(); 528 529 RT mmax = m.maxAbs2Element(); 530 RT norm; 531 if (mmax == RT(0)) { 532 // Then norm is also 0 533 norm = RT(0); 534 } else if (TMV_Underflow(mmax * mmax)) { 535 // Then we need to rescale, since underflow has caused 536 // rounding errors. 537 // Epsilon is a pure power of 2, so no rounding errors from 538 // rescaling. 539 const RT inveps = RT(1)/eps; 540 RT scale = inveps; 541 mmax *= scale; 542 const RT eps2 = eps*eps; 543 while (mmax < eps2) { scale *= inveps; mmax *= inveps; } 544 norm = TMV_SQRT(m.normSq(scale))/scale; 545 } else if (RT(1) / mmax == RT(0)) { 546 // Then mmax is already inf, so no hope of making it more accurate. 547 norm = mmax; 548 } else if (RT(1) / (mmax*mmax) == RT(0)) { 549 // Then we have overflow, so we need to rescale: 550 const RT inveps = RT(1)/eps; 551 RT scale = eps; 552 mmax *= scale; 553 while (mmax > inveps) { scale *= eps; mmax *= eps; } 554 norm = TMV_SQRT(m.normSq(scale))/scale; 555 } else { 556 norm = TMV_SQRT(m.normSq()); 557 } 558 return norm; 559 } 560 561 template <class T> NonLapNormInf(const GenMatrix<T> & m)562 static inline RT NonLapNormInf(const GenMatrix<T>& m) 563 { return NonLapNorm1(m.transpose()); } 564 565 #ifdef XLAP 566 template <class T> LapNorm(const char c,const GenMatrix<T> & m)567 static RT LapNorm(const char c, const GenMatrix<T>& m) 568 { 569 switch(c) { 570 case 'M' : return NonLapMaxAbsElement(m); 571 case '1' : return NonLapNorm1(m); 572 case 'F' : return NonLapNormF(m); 573 case 'I' : return NonLapNormInf(m); 574 default : TMVAssert(TMV_FALSE); 575 } 576 return RT(0); 577 } 578 #ifdef INST_DOUBLE 579 template <> LapNorm(const char c,const GenMatrix<double> & m)580 double LapNorm(const char c, const GenMatrix<double>& m) 581 { 582 TMVAssert(m.iscm()); 583 char cc = c; 584 int M = m.colsize(); 585 int N = m.rowsize(); 586 int lda = m.stepj(); 587 #ifndef LAPNOWORK 588 int lwork = c=='I' ? M : 0; 589 AlignedArray<double> work(lwork); 590 VectorViewOf(work.get(),lwork).setZero(); 591 #endif 592 double norm = LAPNAME(dlange) ( 593 LAPCM LAPV(cc),LAPV(M),LAPV(N), 594 LAPP(m.cptr()),LAPV(lda) LAPWK(work.get())); 595 return norm; 596 } 597 template <> LapNorm(const char c,const GenMatrix<std::complex<double>> & m)598 double LapNorm(const char c, const GenMatrix<std::complex<double> >& m) 599 { 600 TMVAssert(m.iscm()); 601 char cc = c; 602 int M = m.colsize(); 603 int N = m.rowsize(); 604 int lda = m.stepj(); 605 #ifndef LAPNOWORK 606 int lwork = c=='I' ? M : 0; 607 AlignedArray<double> work(lwork); 608 VectorViewOf(work.get(),lwork).setZero(); 609 #endif 610 double norm = LAPNAME(zlange) ( 611 LAPCM LAPV(cc),LAPV(M),LAPV(N), 612 LAPP(m.cptr()),LAPV(lda) LAPWK(work.get())); 613 return norm; 614 } 615 #endif 616 #ifdef INST_FLOAT 617 template <> LapNorm(const char c,const GenMatrix<float> & m)618 float LapNorm(const char c, const GenMatrix<float>& m) 619 { 620 TMVAssert(m.iscm()); 621 char cc = c; 622 int M = m.colsize(); 623 int N = m.rowsize(); 624 int lda = m.stepj(); 625 #ifndef LAPNOWORK 626 int lwork = c=='I' ? M : 0; 627 AlignedArray<float> work(lwork); 628 VectorViewOf(work.get(),lwork).setZero(); 629 #endif 630 double norm = LAPNAME(slange) ( 631 LAPCM LAPV(cc),LAPV(M),LAPV(N), 632 LAPP(m.cptr()),LAPV(lda) LAPWK(work.get())); 633 return norm; 634 } 635 template <> LapNorm(const char c,const GenMatrix<std::complex<float>> & m)636 float LapNorm(const char c, const GenMatrix<std::complex<float> >& m) 637 { 638 TMVAssert(m.iscm()); 639 char cc = c; 640 int M = m.colsize(); 641 int N = m.rowsize(); 642 int lda = m.stepj(); 643 #ifndef LAPNOWORK 644 int lwork = c=='I' ? M : 0; 645 AlignedArray<float> work(lwork); 646 VectorViewOf(work.get(),lwork).setZero(); 647 #endif 648 double norm = LAPNAME(clange) ( 649 LAPCM LAPV(cc),LAPV(M),LAPV(N), 650 LAPP(m.cptr()),LAPV(lda) LAPWK(work.get())); 651 return norm; 652 } 653 #endif 654 #endif // XLAP 655 656 #ifdef INST_INT NonLapNormF(const GenMatrix<int> &)657 static int NonLapNormF(const GenMatrix<int>& ) 658 { TMVAssert(TMV_FALSE); return 0; } NonLapNormF(const GenMatrix<std::complex<int>> &)659 static int NonLapNormF(const GenMatrix<std::complex<int> >& ) 660 { TMVAssert(TMV_FALSE); return 0; } NonLapNorm1(const GenMatrix<std::complex<int>> &)661 static int NonLapNorm1(const GenMatrix<std::complex<int> >& ) 662 { TMVAssert(TMV_FALSE); return 0; } NonLapMaxAbsElement(const GenMatrix<std::complex<int>> &)663 static int NonLapMaxAbsElement(const GenMatrix<std::complex<int> >& ) 664 { TMVAssert(TMV_FALSE); return 0; } 665 #endif 666 667 template <class T> maxAbsElement() const668 RT GenMatrix<T>::maxAbsElement() const 669 { 670 #ifdef XLAP 671 if (isrm() && stepi() > 0) return LapNorm('M',transpose()); 672 else if (iscm() && stepj() > 0) return LapNorm('M',*this); 673 else 674 #endif 675 return NonLapMaxAbsElement(*this); 676 } 677 template <class T> maxAbs2Element() const678 RT GenMatrix<T>::maxAbs2Element() const 679 { 680 #ifdef XLAP 681 if (Traits<T>::iscomplex) return NonLapMaxAbs2Element(*this); 682 else if (isrm() && stepi() > 0) return LapNorm('M',transpose()); 683 else if (iscm() && stepj() > 0) return LapNorm('M',*this); 684 else 685 #endif 686 return NonLapMaxAbs2Element(*this); 687 } 688 template <class T> norm1() const689 RT GenMatrix<T>::norm1() const 690 { 691 #ifdef XLAP 692 if (isrm() && stepi() > 0) return LapNorm('I',transpose()); 693 else if (iscm() && stepj() > 0) return LapNorm('1',*this); 694 else 695 #endif 696 return NonLapNorm1(*this); 697 } 698 template <class T> normF() const699 RT GenMatrix<T>::normF() const 700 { 701 #ifdef XLAP 702 if (isrm() && stepi() > 0) return LapNorm('F',transpose()); 703 else if (iscm() && stepj() > 0) return LapNorm('F',*this); 704 else 705 #endif 706 return NonLapNormF(*this); 707 } 708 709 template <class T> DoNorm2(const GenMatrix<T> & m)710 static RT DoNorm2(const GenMatrix<T>& m) 711 { 712 if (m.colsize() < m.rowsize()) return DoNorm2(m.transpose()); 713 if (m.rowsize() == 0) return RT(0); 714 Matrix<T> m2(m); 715 DiagMatrix<RT> S(m.rowsize()); 716 SV_Decompose(m2.view(),S.view(),false); 717 return S(0); 718 } 719 720 template <class T> DoCondition(const GenMatrix<T> & m)721 static RT DoCondition(const GenMatrix<T>& m) 722 { 723 if (m.colsize() < m.rowsize()) return DoCondition(m.transpose()); 724 if (m.rowsize() == 0) return RT(1); 725 Matrix<T> m2(m); 726 DiagMatrix<RT> S(m.rowsize()); 727 SV_Decompose(m2.view(),S.view(),false); 728 return S(0)/S(S.size()-1); 729 } 730 731 #ifdef INST_INT DoNorm2(const GenMatrix<int> &)732 static int DoNorm2(const GenMatrix<int>& ) 733 { TMVAssert(TMV_FALSE); return 0; } DoCondition(const GenMatrix<int> &)734 static int DoCondition(const GenMatrix<int>& ) 735 { TMVAssert(TMV_FALSE); return 0; } DoNorm2(const GenMatrix<std::complex<int>> &)736 static int DoNorm2(const GenMatrix<std::complex<int> >& ) 737 { TMVAssert(TMV_FALSE); return 0; } DoCondition(const GenMatrix<std::complex<int>> &)738 static int DoCondition(const GenMatrix<std::complex<int> >& ) 739 { TMVAssert(TMV_FALSE); return 0; } 740 #endif 741 742 template <class T> doNorm2() const743 RT GenMatrix<T>::doNorm2() const 744 { return tmv::DoNorm2(*this); } 745 746 template <class T> doCondition() const747 RT GenMatrix<T>::doCondition() const 748 { return tmv::DoCondition(*this); } 749 750 template <class T> QInverse() const751 QuotXM<T,T> GenMatrix<T>::QInverse() const 752 { return QuotXM<T,T>(T(1),*this); } 753 754 // 755 // Modifying Functions 756 // 757 758 template <class T, int A> clip(RT thresh)759 MatrixView<T,A>& MatrixView<T,A>::clip(RT thresh) 760 { 761 TMVAssert(A==CStyle); 762 if (this->canLinearize()) linearView().clip(thresh); 763 else { 764 if (this->isrm()) { 765 const ptrdiff_t M = colsize(); 766 for(ptrdiff_t i=0;i<M;++i) row(i).clip(thresh); 767 } else { 768 const ptrdiff_t N = rowsize(); 769 for(ptrdiff_t j=0;j<N;++j) col(j).clip(thresh); 770 } 771 } 772 return *this; 773 } 774 775 template <class T, int A> setZero()776 MatrixView<T,A>& MatrixView<T,A>::setZero() 777 { 778 TMVAssert(A==CStyle); 779 if (this->canLinearize()) linearView().setZero(); 780 else { 781 if (this->isrm()) { 782 const ptrdiff_t M = colsize(); 783 for(ptrdiff_t i=0;i<M;++i) row(i).setZero(); 784 } else { 785 const ptrdiff_t N = rowsize(); 786 for(ptrdiff_t j=0;j<N;++j) col(j).setZero(); 787 } 788 } 789 return *this; 790 } 791 792 template <class T, int A> setAllTo(const T & x)793 MatrixView<T,A>& MatrixView<T,A>::setAllTo(const T& x) 794 { 795 TMVAssert(A==CStyle); 796 if (this->canLinearize()) linearView().setAllTo(x); 797 else { 798 if (this->isrm()) { 799 const ptrdiff_t M = colsize(); 800 for(ptrdiff_t i=0;i<M;++i) row(i).setAllTo(x); 801 } else { 802 const ptrdiff_t N = rowsize(); 803 for(ptrdiff_t j=0;j<N;++j) col(j).setAllTo(x); 804 } 805 } 806 return *this; 807 } 808 809 template <class T, int A> addToAll(const T & x)810 MatrixView<T,A>& MatrixView<T,A>::addToAll(const T& x) 811 { 812 TMVAssert(A==CStyle); 813 if (this->canLinearize()) linearView().addToAll(x); 814 else { 815 if (this->isrm()) { 816 const ptrdiff_t M = colsize(); 817 for(ptrdiff_t i=0;i<M;++i) row(i).addToAll(x); 818 } else { 819 const ptrdiff_t N = rowsize(); 820 for(ptrdiff_t j=0;j<N;++j) col(j).addToAll(x); 821 } 822 } 823 return *this; 824 } 825 826 template <class T, int A> transposeSelf()827 MatrixView<T,A>& MatrixView<T,A>::transposeSelf() 828 { 829 TMVAssert(A==CStyle); 830 TMVAssert(colsize() == rowsize()); 831 const ptrdiff_t M = colsize(); 832 for(ptrdiff_t i=1;i<M;++i) Swap(row(i,0,i),col(i,0,i)); 833 return *this; 834 } 835 836 template <class T, int A> conjugateSelf()837 MatrixView<T,A>& MatrixView<T,A>::conjugateSelf() 838 { 839 TMVAssert(A==CStyle); 840 if (isComplex(T())) { 841 if (this->canLinearize()) linearView().conjugateSelf(); 842 else { 843 if (this->isrm()) { 844 const ptrdiff_t M = colsize(); 845 for(ptrdiff_t i=0;i<M;++i) row(i).conjugateSelf(); 846 } else { 847 const ptrdiff_t N = rowsize(); 848 for(ptrdiff_t j=0;j<N;++j) col(j).conjugateSelf(); 849 } 850 } 851 } 852 return *this; 853 } 854 855 template <class T, int A> setToIdentity(const T & x)856 MatrixView<T,A>& MatrixView<T,A>::setToIdentity(const T& x) 857 { 858 TMVAssert(A==CStyle); 859 TMVAssert(colsize() == rowsize()); 860 setZero(); 861 diag().setAllTo(x); 862 return *this; 863 } 864 865 template <class T, int A> permuteRows(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)866 MatrixView<T,A>& MatrixView<T,A>::permuteRows( 867 const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 868 { 869 TMVAssert(A==CStyle); 870 TMVAssert(i2<=colsize()); 871 TMVAssert(i1<=i2); 872 // This idea of doing the permutations a block at a time 873 // is cribbed from the LAPack code. It does speed things up 874 // quite a bit for large matrices. On my machine where BLOCKSIZE=64 875 // is optimal for most routines, blocks of 32 were optimal here, 876 // so I use BLOCKSIZE/2 in general. 877 const ptrdiff_t N = rowsize(); 878 const ptrdiff_t Nx = N/PERM_BLOCKSIZE*PERM_BLOCKSIZE; 879 if (Nx != 0) { 880 for(ptrdiff_t j=0;j<Nx;) { 881 ptrdiff_t j2 = j+PERM_BLOCKSIZE; 882 const ptrdiff_t* pi = p+i1; 883 for(ptrdiff_t i=i1;i<i2;++i,++pi) { 884 TMVAssert(*pi < colsize()); 885 colRange(j,j2).swapRows(i,*pi); 886 } 887 j = j2; 888 } 889 } 890 if (Nx != N) { 891 const ptrdiff_t* pi = p+i1; 892 for(ptrdiff_t i=i1;i<i2;++i,++pi) { 893 TMVAssert(*pi < colsize()); 894 colRange(Nx,N).swapRows(i,*pi); 895 } 896 } 897 return *this; 898 } 899 900 template <class T, int A> reversePermuteRows(const ptrdiff_t * p,ptrdiff_t i1,ptrdiff_t i2)901 MatrixView<T,A>& MatrixView<T,A>::reversePermuteRows( 902 const ptrdiff_t* p, ptrdiff_t i1, ptrdiff_t i2) 903 { 904 TMVAssert(A==CStyle); 905 TMVAssert(i2<=colsize()); 906 TMVAssert(i1<=i2); 907 const ptrdiff_t N = rowsize(); 908 const ptrdiff_t Nx = N/PERM_BLOCKSIZE*PERM_BLOCKSIZE; 909 if (Nx != 0) { 910 for(ptrdiff_t j=0;j<Nx;) { 911 ptrdiff_t j2 = j+PERM_BLOCKSIZE; 912 const ptrdiff_t* pi = p+i2; 913 for(ptrdiff_t i=i2;i>i1;) { 914 --i; --pi; 915 TMVAssert(*pi < colsize()); 916 colRange(j,j2).swapRows(i,*pi); 917 } 918 j = j2; 919 } 920 } 921 if (Nx != N) { 922 const ptrdiff_t* pi = p+i2; 923 for(ptrdiff_t i=i2;i>i1;) { 924 --i; --pi; 925 TMVAssert(*pi < colsize()); 926 colRange(Nx,N).swapRows(i,*pi); 927 } 928 } 929 return *this; 930 } 931 932 // 933 // Copy Matrices 934 // 935 936 template <class T> NonLapCopy(const GenMatrix<T> & m1,MatrixView<T> m2)937 static void NonLapCopy(const GenMatrix<T>& m1, MatrixView<T> m2) 938 { 939 TMVAssert(m2.rowsize() == m1.rowsize()); 940 TMVAssert(m2.colsize() == m1.colsize()); 941 TMVAssert(m1.ct()==NonConj); 942 TMVAssert(m2.ct()==NonConj); 943 const ptrdiff_t M = m2.colsize(); 944 const ptrdiff_t N = m2.rowsize(); 945 946 if (m1.iscm() && m2.iscm()) { 947 const T* p1 = m1.cptr(); 948 T* p2 = m2.ptr(); 949 const ptrdiff_t s1 = m1.stepj(); 950 const ptrdiff_t s2 = m2.stepj(); 951 for(ptrdiff_t j=0;j<N;++j,p1+=s1,p2+=s2) { 952 std::copy(p1,p1+M,p2); 953 } 954 } else if (M > N) { 955 if (shouldReverse(m1.stepi(),m2.stepi())) 956 for(ptrdiff_t j=0;j<N;++j) 957 DoCopySameType(m1.col(j).reverse(),m2.col(j).reverse()); 958 else 959 for(ptrdiff_t j=0;j<N;++j) 960 DoCopySameType(m1.col(j),m2.col(j)); 961 } else { 962 if (shouldReverse(m1.stepj(),m2.stepj())) 963 for(ptrdiff_t i=0;i<M;++i) 964 DoCopySameType(m1.row(i).reverse(),m2.row(i).reverse()); 965 else 966 for(ptrdiff_t i=0;i<M;++i) 967 DoCopySameType(m1.row(i),m2.row(i)); 968 } 969 } 970 #ifdef ELAP 971 template <class T> LapCopy(const GenMatrix<T> & m1,MatrixView<T> m2)972 static inline void LapCopy(const GenMatrix<T>& m1, MatrixView<T> m2) 973 { NonLapCopy(m1,m2); } 974 #ifdef INST_DOUBLE 975 template <> LapCopy(const GenMatrix<double> & m1,MatrixView<double> m2)976 void LapCopy(const GenMatrix<double>& m1, MatrixView<double> m2) 977 { 978 TMVAssert(m2.rowsize() == m1.rowsize()); 979 TMVAssert(m2.colsize() == m1.colsize()); 980 TMVAssert(m1.iscm()); 981 TMVAssert(m2.iscm()); 982 char c = 'A'; 983 int m = m1.colsize(); 984 int n = m1.rowsize(); 985 int ld1 = m1.stepj(); 986 int ld2 = m2.stepj(); 987 LAPNAME(dlacpy) ( 988 LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1), 989 LAPP(m2.ptr()),LAPV(ld2)); 990 } 991 template <> LapCopy(const GenMatrix<std::complex<double>> & m1,MatrixView<std::complex<double>> m2)992 void LapCopy( 993 const GenMatrix<std::complex<double> >& m1, 994 MatrixView<std::complex<double> > m2) 995 { 996 TMVAssert(m2.rowsize() == m1.rowsize()); 997 TMVAssert(m2.colsize() == m1.colsize()); 998 TMVAssert(m1.iscm()); 999 TMVAssert(m2.iscm()); 1000 TMVAssert(m1.ct() == NonConj); 1001 TMVAssert(m2.ct() == NonConj); 1002 char c = 'A'; 1003 int m = m1.colsize(); 1004 int n = m1.rowsize(); 1005 int ld1 = m1.stepj(); 1006 int ld2 = m2.stepj(); 1007 LAPNAME(zlacpy) ( 1008 LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1), 1009 LAPP(m2.ptr()),LAPV(ld2)); 1010 } 1011 #endif 1012 #ifdef INST_FLOAT 1013 template <> LapCopy(const GenMatrix<float> & m1,MatrixView<float> m2)1014 void LapCopy(const GenMatrix<float>& m1, MatrixView<float> m2) 1015 { 1016 TMVAssert(m2.rowsize() == m1.rowsize()); 1017 TMVAssert(m2.colsize() == m1.colsize()); 1018 TMVAssert(m1.iscm()); 1019 TMVAssert(m2.iscm()); 1020 char c = 'A'; 1021 int m = m1.colsize(); 1022 int n = m1.rowsize(); 1023 int ld1 = m1.stepj(); 1024 int ld2 = m2.stepj(); 1025 LAPNAME(slacpy) ( 1026 LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1), 1027 LAPP(m2.ptr()),LAPV(ld2)); 1028 } 1029 template <> LapCopy(const GenMatrix<std::complex<float>> & m1,MatrixView<std::complex<float>> m2)1030 void LapCopy( 1031 const GenMatrix<std::complex<float> >& m1, 1032 MatrixView<std::complex<float> > m2) 1033 { 1034 TMVAssert(m2.rowsize() == m1.rowsize()); 1035 TMVAssert(m2.colsize() == m1.colsize()); 1036 TMVAssert(m1.iscm()); 1037 TMVAssert(m2.iscm()); 1038 TMVAssert(m1.ct() == NonConj); 1039 TMVAssert(m2.ct() == NonConj); 1040 char c = 'A'; 1041 int m = m1.colsize(); 1042 int n = m1.rowsize(); 1043 int ld1 = m1.stepj(); 1044 int ld2 = m2.stepj(); 1045 LAPNAME(clacpy) ( 1046 LAPCM LAPV(c),LAPV(m),LAPV(n),LAPP(m1.cptr()),LAPV(ld1), 1047 LAPP(m2.ptr()),LAPV(ld2)); 1048 } 1049 #endif 1050 #endif 1051 template <class T> DoCopySameType(const GenMatrix<T> & m1,MatrixView<T> m2)1052 void DoCopySameType(const GenMatrix<T>& m1, MatrixView<T> m2) 1053 { 1054 TMVAssert(m2.rowsize() == m1.rowsize()); 1055 TMVAssert(m2.colsize() == m1.colsize()); 1056 TMVAssert(m2.rowsize() > 0); 1057 TMVAssert(m2.colsize() > 0); 1058 TMVAssert(m1.ct() == NonConj); 1059 TMVAssert(m2.ct() == NonConj); 1060 TMVAssert(!m2.isSameAs(m1)); 1061 1062 #ifdef ELAP 1063 if (m1.iscm() && m2.iscm() && m1.stepj()>0 && m2.stepj()>0) 1064 LapCopy(m1,m2); 1065 else 1066 #endif 1067 NonLapCopy(m1,m2); 1068 } 1069 1070 // 1071 // Swap 1072 // 1073 1074 template <class T> Swap(MatrixView<T> m1,MatrixView<T> m2)1075 void Swap(MatrixView<T> m1, MatrixView<T> m2) 1076 { 1077 TMVAssert(m1.colsize() == m2.colsize()); 1078 TMVAssert(m1.rowsize() == m2.rowsize()); 1079 if (m1.canLinearize() && m2.canLinearize() && 1080 m1.stepi() == m2.stepi() && m1.stepj() == m2.stepj()) { 1081 TMVAssert(m1.linearView().size() == m2.linearView().size()); 1082 TMVAssert(m1.stepi()==m2.stepi() && m1.stepj()==m2.stepj()); 1083 Swap(m1.linearView(),m2.linearView()); 1084 } else { 1085 if (m1.isrm() && m2.isrm()) { 1086 const ptrdiff_t M = m1.colsize(); 1087 for(ptrdiff_t i=0;i<M;++i) Swap(m1.row(i),m2.row(i)); 1088 } else { 1089 const ptrdiff_t N = m1.rowsize(); 1090 for(ptrdiff_t j=0;j<N;++j) Swap(m1.col(j),m2.col(j)); 1091 } 1092 } 1093 } 1094 1095 // 1096 // m1 == m2 1097 // 1098 1099 template <class T1, class T2> operator ==(const GenMatrix<T1> & m1,const GenMatrix<T2> & m2)1100 bool operator==(const GenMatrix<T1>& m1, const GenMatrix<T2>& m2) 1101 { 1102 if (m1.colsize() != m2.colsize()) return false; 1103 else if (m1.rowsize() != m2.rowsize()) return false; 1104 else if (m1.isSameAs(m2)) return true; 1105 else if (m1.stepi()==m2.stepi() && m1.stepj()==m2.stepj() && 1106 m1.canLinearize() && m2.canLinearize()) 1107 return m1.constLinearView() == m2.constLinearView(); 1108 else { 1109 const ptrdiff_t M = m1.colsize(); 1110 for(ptrdiff_t i=0;i<M;++i) 1111 if (m1.row(i) != m2.row(i)) return false; 1112 return true; 1113 } 1114 } 1115 1116 // 1117 // I/O 1118 // 1119 1120 template <class T> write(const TMV_Writer & writer) const1121 void GenMatrix<T>::write(const TMV_Writer& writer) const 1122 { 1123 const ptrdiff_t M = colsize(); 1124 const ptrdiff_t N = rowsize(); 1125 writer.begin(); 1126 writer.writeCode("M"); 1127 writer.writeSize(M); 1128 writer.writeSize(N); 1129 writer.writeStart(); 1130 for(ptrdiff_t i=0;i<M;++i) { 1131 writer.writeLParen(); 1132 for(ptrdiff_t j=0;j<N;++j) { 1133 if (j > 0) writer.writeSpace(); 1134 writer.writeValue(cref(i,j)); 1135 } 1136 writer.writeRParen(); 1137 if (i < M-1) writer.writeRowEnd(); 1138 } 1139 writer.writeFinal(); 1140 writer.end(); 1141 } 1142 1143 #ifndef NOTHROW 1144 template <class T> 1145 class MatrixReadError : public ReadError 1146 { 1147 public : 1148 Matrix<T> m; 1149 ptrdiff_t i,j; 1150 std::string exp,got; 1151 ptrdiff_t cs,rs; 1152 bool is,iseof,isbad; 1153 MatrixReadError(std::istream & _is)1154 MatrixReadError(std::istream& _is) throw() : 1155 ReadError("Matrix."), 1156 i(0), j(0), cs(0), rs(0), 1157 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} MatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)1158 MatrixReadError( 1159 std::istream& _is, 1160 const std::string& _e, const std::string& _g) throw() : 1161 ReadError("Matrix."), 1162 i(0), j(0), exp(_e), got(_g), cs(0), rs(0), 1163 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 1164 MatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenMatrix<T> & _m,std::istream & _is)1165 MatrixReadError( 1166 ptrdiff_t _i, ptrdiff_t _j, const GenMatrix<T>& _m, std::istream& _is) throw() : 1167 ReadError("Matrix."), 1168 m(_m), i(_i), j(_j), 1169 cs(m.colsize()), rs(m.rowsize()), 1170 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} MatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)1171 MatrixReadError( 1172 ptrdiff_t _i, ptrdiff_t _j, const GenMatrix<T>& _m, std::istream& _is, 1173 const std::string& _e, const std::string& _g) throw() : 1174 ReadError("Matrix."), 1175 m(_m), i(_i), j(_j), exp(_e), got(_g), 1176 cs(m.colsize()), rs(m.rowsize()), 1177 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} MatrixReadError(const GenMatrix<T> & _m,std::istream & _is,ptrdiff_t _cs,ptrdiff_t _rs)1178 MatrixReadError( 1179 const GenMatrix<T>& _m, 1180 std::istream& _is, ptrdiff_t _cs, ptrdiff_t _rs) throw() : 1181 ReadError("Matrix."), 1182 m(_m), i(0), j(0), cs(_cs), rs(_rs), 1183 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 1184 MatrixReadError(const MatrixReadError<T> & rhs)1185 MatrixReadError(const MatrixReadError<T>& rhs) : 1186 m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got), 1187 cs(rhs.cs), rs(rhs.rs), 1188 is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {} ~MatrixReadError()1189 ~MatrixReadError() throw() {} 1190 write(std::ostream & os) const1191 void write(std::ostream& os) const throw() 1192 { 1193 os<<"TMV Read Error: Reading istream input for Matrix\n"; 1194 if (exp != got) { 1195 os<<"Wrong format: expected '"<<exp<<"', got '"<<got<<"'.\n"; 1196 } 1197 if (cs != m.colsize()) { 1198 os<<"Wrong column size: expected "<<m.colsize()<< 1199 ", got "<<cs<<".\n"; 1200 } 1201 if (rs != m.rowsize()) { 1202 os<<"Wrong row size: expected "<<m.rowsize()<< 1203 ", got "<<rs<<".\n"; 1204 } 1205 if (!is) { 1206 if (iseof) { 1207 os<<"Input stream reached end-of-file prematurely.\n"; 1208 } else if (isbad) { 1209 os<<"Input stream is corrupted.\n"; 1210 } else { 1211 os<<"Input stream cannot read next character.\n"; 1212 } 1213 } 1214 if (m.colsize() > 0 || m.rowsize() > 0) { 1215 const ptrdiff_t N = m.rowsize(); 1216 os<<"The portion of the Matrix which was successfully " 1217 "read is: \n"; 1218 for(ptrdiff_t ii=0;ii<i;++ii) { 1219 os<<"( "; 1220 for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' '; 1221 os<<" )\n"; 1222 } 1223 os<<"( "; 1224 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' '; 1225 os<<" )\n"; 1226 } 1227 } 1228 }; 1229 #endif 1230 1231 template <class T> FinishRead(const TMV_Reader & reader,MatrixView<T> m)1232 static void FinishRead(const TMV_Reader& reader, MatrixView<T> m) 1233 { 1234 const ptrdiff_t M = m.colsize(); 1235 const ptrdiff_t N = m.rowsize(); 1236 std::string exp, got; 1237 T temp; 1238 if (!reader.readStart(exp,got)) { 1239 #ifdef NOTHROW 1240 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1241 exit(1); 1242 #else 1243 throw MatrixReadError<T>(0,0,m,reader.getis(),exp,got); 1244 #endif 1245 } 1246 for(ptrdiff_t i=0;i<M;++i) { 1247 if (!reader.readLParen(exp,got)) { 1248 #ifdef NOTHROW 1249 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1250 exit(1); 1251 #else 1252 throw MatrixReadError<T>(i,0,m,reader.getis(),exp,got); 1253 #endif 1254 } 1255 for(ptrdiff_t j=0;j<N;++j) { 1256 if (j>0) { 1257 if (!reader.readSpace(exp,got)) { 1258 #ifdef NOTHROW 1259 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1260 exit(1); 1261 #else 1262 throw MatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1263 #endif 1264 } 1265 } 1266 if (!reader.readValue(temp)) { 1267 #ifdef NOTHROW 1268 std::cerr<<"Matrix Read Error: reading value\n"; 1269 exit(1); 1270 #else 1271 throw MatrixReadError<T>(i,j,m,reader.getis()); 1272 #endif 1273 } 1274 m.ref(i,j) = temp; 1275 } 1276 if (!reader.readRParen(exp,got)) { 1277 #ifdef NOTHROW 1278 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1279 exit(1); 1280 #else 1281 throw MatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1282 #endif 1283 } 1284 if (i < M-1 && !reader.readRowEnd(exp,got)) { 1285 #ifdef NOTHROW 1286 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1287 exit(1); 1288 #else 1289 throw MatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1290 #endif 1291 } 1292 } 1293 if (!reader.readFinal(exp,got)) { 1294 #ifdef NOTHROW 1295 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1296 exit(1); 1297 #else 1298 throw MatrixReadError<T>(M,0,m,reader.getis(),exp,got); 1299 #endif 1300 } 1301 } 1302 1303 template <class T, int A> read(const TMV_Reader & reader)1304 void Matrix<T,A>::read(const TMV_Reader& reader) 1305 { 1306 std::string exp,got; 1307 if (!reader.readCode("M",exp,got)) { 1308 #ifdef NOTHROW 1309 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1310 exit(1); 1311 #else 1312 throw MatrixReadError<T>(reader.getis(),exp,got); 1313 #endif 1314 } 1315 ptrdiff_t cs=colsize(), rs=rowsize(); 1316 if (!reader.readSize(cs,exp,got) || 1317 !reader.readSize(rs,exp,got)) { 1318 #ifdef NOTHROW 1319 std::cerr<<"Matrix Read Error: reading size\n"; 1320 exit(1); 1321 #else 1322 throw MatrixReadError<T>(reader.getis(),exp,got); 1323 #endif 1324 } 1325 if (cs != colsize() || rs != rowsize()) resize(cs,rs); 1326 MatrixView<T> v = view(); 1327 FinishRead(reader,v); 1328 } 1329 1330 template <class T, int A> read(const TMV_Reader & reader)1331 void MatrixView<T,A>::read(const TMV_Reader& reader) 1332 { 1333 std::string exp,got; 1334 if (!reader.readCode("M",exp,got)) { 1335 #ifdef NOTHROW 1336 std::cerr<<"Matrix Read Error: "<<got<<" != "<<exp<<std::endl; 1337 exit(1); 1338 #else 1339 throw MatrixReadError<T>(reader.getis(),exp,got); 1340 #endif 1341 } 1342 ptrdiff_t cs=colsize(), rs=rowsize(); 1343 if (!reader.readSize(cs,exp,got) || 1344 !reader.readSize(rs,exp,got)) { 1345 #ifdef NOTHROW 1346 std::cerr<<"Matrix Read Error: reading size\n"; 1347 exit(1); 1348 #else 1349 throw MatrixReadError<T>(reader.getis(),exp,got); 1350 #endif 1351 } 1352 if (cs != colsize() || rs != rowsize()) { 1353 #ifdef NOTHROW 1354 std::cerr<<"Matrix Read Error: wrong size\n"; 1355 exit(1); 1356 #else 1357 throw MatrixReadError<T>(*this,reader.getis(),cs,rs); 1358 #endif 1359 } 1360 FinishRead(reader,*this); 1361 } 1362 1363 #undef RT 1364 1365 #define InstFile "TMV_Matrix.inst" 1366 #include "TMV_Inst.h" 1367 #undef InstFile 1368 1369 } // namespace tmv 1370 1371 1372