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_BandMatrix.h" 26 #include "tmv/TMV_BandLUD.h" 27 #include "tmv/TMV_BandQRD.h" 28 #include "tmv/TMV_BandSVD.h" 29 #include "tmv/TMV_VIt.h" 30 #include "tmv/TMV_BandMatrixArith.h" 31 #include "tmv/TMV_DiagMatrix.h" 32 #include "TMV_IntegerDet.h" 33 #include <iostream> 34 35 namespace tmv { 36 37 #define RT TMV_RealType(T) 38 39 // 40 // Access 41 // 42 43 template <class T> cref(ptrdiff_t i,ptrdiff_t j) const44 T GenBandMatrix<T>::cref(ptrdiff_t i, ptrdiff_t j) const 45 { 46 if (okij(i,j)) { 47 const T* mi = cptr()+i*stepi()+j*stepj(); 48 return isconj() ? TMV_CONJ(*mi) : *mi; 49 } else { 50 return T(0); 51 } 52 } 53 54 template <class T, int A> ref(ptrdiff_t i,ptrdiff_t j)55 typename BandMatrixView<T,A>::reference BandMatrixView<T,A>::ref( 56 ptrdiff_t i, ptrdiff_t j) 57 { 58 T* mi = ptr()+i*stepi()+j*stepj(); 59 return RefHelper<T>::makeRef(mi,ct()); 60 } 61 62 template <class T> setDiv() const63 void GenBandMatrix<T>::setDiv() const 64 { 65 if (!this->divIsSet()) { 66 DivType dt = this->getDivType(); 67 TMVAssert(dt == tmv::LU || dt == tmv::QR || dt == tmv::SV); 68 switch (dt) { 69 case LU : 70 this->divider.reset( 71 new BandLUDiv<T>(*this,this->divIsInPlace())); break; 72 case QR : 73 this->divider.reset( 74 new BandQRDiv<T>(*this,this->divIsInPlace())); break; 75 case SV : 76 this->divider.reset(new BandSVDiv<T>(*this)); break; 77 default : 78 // The above assert should have already failed 79 // so go ahead and fall through. 80 break; 81 } 82 } 83 } 84 85 #ifdef INST_INT 86 template <> setDiv() const87 void GenBandMatrix<int>::setDiv() const 88 { TMVAssert(TMV_FALSE); } 89 template <> setDiv() const90 void GenBandMatrix<std::complex<int> >::setDiv() const 91 { TMVAssert(TMV_FALSE); } 92 #endif 93 94 template <class T> divIsLUDiv() const95 bool GenBandMatrix<T>::divIsLUDiv() const 96 { return dynamic_cast<const BandLUDiv<T>*>(this->getDiv()); } 97 98 template <class T> divIsQRDiv() const99 bool GenBandMatrix<T>::divIsQRDiv() const 100 { return dynamic_cast<const BandQRDiv<T>*>(this->getDiv()); } 101 102 template <class T> divIsSVDiv() const103 bool GenBandMatrix<T>::divIsSVDiv() const 104 { return dynamic_cast<const BandSVDiv<T>*>(this->getDiv()); } 105 106 #ifdef INST_INT 107 template <> divIsLUDiv() const108 bool GenBandMatrix<int>::divIsLUDiv() const 109 { return false; } 110 template <> divIsQRDiv() const111 bool GenBandMatrix<int>::divIsQRDiv() const 112 { return false; } 113 template <> divIsSVDiv() const114 bool GenBandMatrix<int>::divIsSVDiv() const 115 { return false; } 116 117 template <> divIsLUDiv() const118 bool GenBandMatrix<std::complex<int> >::divIsLUDiv() const 119 { return false; } 120 template <> divIsQRDiv() const121 bool GenBandMatrix<std::complex<int> >::divIsQRDiv() const 122 { return false; } 123 template <> divIsSVDiv() const124 bool GenBandMatrix<std::complex<int> >::divIsSVDiv() const 125 { return false; } 126 #endif 127 BandStorageLength(StorageType s,ptrdiff_t cs,ptrdiff_t rs,ptrdiff_t lo,ptrdiff_t hi)128 ptrdiff_t BandStorageLength(StorageType s, ptrdiff_t cs, ptrdiff_t rs, ptrdiff_t lo, ptrdiff_t hi) 129 { 130 TMVAssert(s == RowMajor || s == ColMajor || s == DiagMajor); 131 if (cs == 0 || rs == 0) return 0; 132 else if (cs == rs) return (cs-1)*(lo+hi)+cs; 133 else { 134 // correct cs, rs to be actual end of data 135 if (cs > rs+lo) cs = rs+lo; 136 if (rs > cs+hi) rs = cs+hi; 137 138 if (s == RowMajor) 139 // si = lo+hi, sj = 1, size = (cs-1)*si + (rs-1)*sj + 1 140 return (cs-1)*(lo+hi) + rs; 141 else if (s == ColMajor) 142 // si = 1, sj = lo+hi, size = (cs-1)*si + (rs-1)*sj + 1 143 return (rs-1)*(lo+hi) + cs; 144 else if (cs > rs) 145 // si = -rs, sj = 1+rs, size = (rs-lo-hi-1)*si + (rs-1)*sj + 1 146 // size = (rs-lo-hi-1)(-rs) + (rs-1)(1+rs) + 1 147 // = rs(lo+hi+1) 148 return rs*(lo+hi+1); 149 else 150 // si = 1-cs, sj = cs, size = (rs-lo-hi-1)*si + (rs-1)*sj + 1 151 // size = (cs-lo-hi-1+rs-cs)(1-cs) + (rs-1)(cs) + 1 152 // = cs(lo+hi+1-rs+rs-1) + rs-lo-hi-1 + 1 153 // = (cs-1)(lo+hi) + rs 154 return (cs-1)*(lo+hi) + rs; 155 } 156 } 157 BandNumElements(ptrdiff_t cs,ptrdiff_t rs,ptrdiff_t lo,ptrdiff_t hi)158 ptrdiff_t BandNumElements(ptrdiff_t cs, ptrdiff_t rs, ptrdiff_t lo, ptrdiff_t hi) 159 { 160 if (cs == 0 || rs == 0) return 0; 161 else if (cs == rs) { 162 return cs*(lo+hi+1) - (lo*(lo+1)/2) - (hi*(hi+1)/2); 163 } else if (cs > rs) { 164 // lox = number of subdiagonals that are clipped. 165 ptrdiff_t lox = TMV_MAX(rs+lo-cs,ptrdiff_t(0)); 166 return rs*(lo+hi+1) - (lox*(lox+1)/2) - (hi*(hi+1)/2); 167 } else { 168 // hix = number of superdiagonals that are clipped. 169 ptrdiff_t hix = TMV_MAX(cs+hi-rs,ptrdiff_t(0)); 170 return cs*(lo+hi+1) - (lo*(lo+1)/2) - (hix*(hix+1)/2); 171 } 172 } 173 174 // 175 // OK? (SubMatrix, etc.) 176 // 177 178 template <class T> hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const179 bool GenBandMatrix<T>::hasSubMatrix( 180 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 181 { 182 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 183 bool ok = true; 184 if (istep == 0) { 185 ok = false; 186 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 187 } 188 if (i1 < 0 || i1 >= colsize()) { 189 ok = false; 190 std::cerr<<"first col element ("<<i1<<") must be in 0 -- "; 191 std::cerr<<colsize()-1<<std::endl; 192 } 193 if (i2-istep < 0 || i2-istep >= colsize()) { 194 ok = false; 195 std::cerr<<"last col element ("<<i2-istep<<") must be in 0 -- "; 196 std::cerr<<colsize()-1<<std::endl; 197 } 198 if ((i2-i1)%istep != 0) { 199 ok = false; 200 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 201 std::cerr<<istep<<")\n"; 202 } 203 if ((i2-i1)/istep < 0) { 204 ok = false; 205 std::cerr<<"n col elements ("<<(i2-i1)/istep<<") must be nonnegative\n"; 206 } 207 if (jstep == 0) { 208 ok = false; 209 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 210 } 211 if (j1 < 0 || j1 >= rowsize()) { 212 ok = false; 213 std::cerr<<"first row element ("<<j1<<") must be in 0 -- "; 214 std::cerr<<rowsize()-1<<std::endl; 215 } 216 if (j2-jstep < 0 || j2-jstep >= rowsize()) { 217 ok = false; 218 std::cerr<<"last row element ("<<j2-jstep<<") must be in 0 -- "; 219 std::cerr<<rowsize()-1<<std::endl; 220 } 221 if ((j2-j1)%jstep != 0) { 222 ok = false; 223 std::cerr<<"row range ("<<j2-j1<<") must be multiple of jstep ("; 224 std::cerr<<jstep<<")\n"; 225 } 226 if ((j2-j1)/jstep < 0) { 227 ok = false; 228 std::cerr<<"n row elements ("<<(j2-j1)/jstep<<") must be nonnegative\n"; 229 } 230 if (!okij(i1,j1)) { 231 ok = false; 232 std::cerr<<"Upper left corner ("<<i1<<','<<j1<<") must be in band\n"; 233 } 234 if (!okij(i1,j2-jstep)) { 235 ok = false; 236 std::cerr<<"Upper right corner ("<<i1<<','<<j2-jstep; 237 std::cerr<<") must be in band\n"; 238 } 239 if (!okij(i2-istep,j1)) { 240 ok = false; 241 std::cerr<<"Lower left corner ("<<i2-istep<<','<<j1; 242 std::cerr<<") must be in band\n"; 243 } 244 if (!okij(i2-istep,j2-jstep)) { 245 ok = false; 246 std::cerr<<"Lower right corner ("<<i2-istep<<','<<j2-jstep; 247 std::cerr<<") must be in band\n"; 248 } 249 return ok; 250 } 251 252 template <class T> hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t size) const253 bool GenBandMatrix<T>::hasSubVector( 254 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 255 { 256 if (size==0) return true; 257 bool ok = true; 258 if (istep == 0 && jstep == 0) { 259 ok = false; 260 std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep; 261 std::cerr<<") can not both be 0\n"; 262 } 263 if (i<0 || i >= colsize()) { 264 ok = false; 265 std::cerr<<"i ("<<i<<") must be in 0 -- "<<colsize()-1<<std::endl; 266 } 267 if (j<0 || j >= rowsize()) { 268 ok = false; 269 std::cerr<<"j ("<<j<<") must be in 0 -- "<<rowsize()-1<<std::endl; 270 } 271 ptrdiff_t i2 = i+istep*(size-1); 272 ptrdiff_t j2 = j+jstep*(size-1); 273 if (i2 < 0 || i2 >= colsize()) { 274 ok = false; 275 std::cerr<<"last element's i ("<<i2<<") must be in 0 -- "; 276 std::cerr<<colsize()-1<<std::endl; 277 } 278 if (j2 < 0 || j2 >= rowsize()) { 279 ok = false; 280 std::cerr<<"last element's j ("<<j2<<") must be in 0 -- "; 281 std::cerr<<rowsize()-1<<std::endl; 282 } 283 if (!okij(i,j)) { 284 ok = false; 285 std::cerr<<"First element ("<<i<<','<<j<<") must be in band\n"; 286 } 287 if (!okij(i2,j2)) { 288 ok = false; 289 std::cerr<<"Last element ("<<i2<<','<<j2<<") must be in band\n"; 290 } 291 return ok; 292 } 293 template <class T> hasSubBandMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t newnlo,ptrdiff_t newnhi,ptrdiff_t istep,ptrdiff_t jstep) const294 bool GenBandMatrix<T>::hasSubBandMatrix( 295 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t newnlo, ptrdiff_t newnhi, 296 ptrdiff_t istep, ptrdiff_t jstep) const 297 { 298 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 299 bool ok = true; 300 if (istep == 0) { 301 ok = false; 302 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 303 } 304 if (i1 < 0 || i1 >= colsize()) { 305 ok = false; 306 std::cerr<<"first col element ("<<i1<<") must be in 0 -- "; 307 std::cerr<<colsize()-1<<std::endl; 308 } 309 if (i2-istep < 0 || i2-istep >= colsize()) { 310 ok = false; 311 std::cerr<<"last col element ("<<i2-istep<<") must be in 0 -- "; 312 std::cerr<<colsize()-1<<std::endl; 313 } 314 if ((i2-i1)%istep != 0) { 315 ok = false; 316 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 317 std::cerr<<istep<<")\n"; 318 } 319 if ((i2-i1)/istep < 0) { 320 ok = false; 321 std::cerr<<"n col elements ("<<(i2-i1)/istep<<") must be nonnegative\n"; 322 } 323 if (jstep == 0) { 324 ok = false; 325 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 326 } 327 if (j1 < 0 || j1 >= rowsize()) { 328 ok = false; 329 std::cerr<<"first row element ("<<j1<<") must be in 0 -- "; 330 std::cerr<<rowsize()-1<<std::endl; 331 } 332 if (j2-jstep < 0 || j2-jstep >= rowsize()) { 333 ok = false; 334 std::cerr<<"last row element ("<<j2-jstep<<") must be in 0 -- "; 335 std::cerr<<rowsize()-1<<std::endl; 336 } 337 if ((j2-j1)%jstep != 0) { 338 ok = false; 339 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 340 std::cerr<<jstep<<")\n"; 341 } 342 if ((j2-j1)/jstep < 0) { 343 ok = false; 344 std::cerr<<"n row elements ("<<(j2-j1)/jstep<<") must be nonnegative\n"; 345 } 346 if (!okij(i1,j1)) { 347 ok = false; 348 std::cerr<<"Upper left corner ("<<i1<<','<<j1<<") must be in band\n"; 349 } 350 if (!okij(i1,j1+newnhi)) { 351 ok = false; 352 std::cerr<<"Start of top diagonal ("<<i1<<','<<j1+newnhi; 353 std::cerr<<") must be in band\n"; 354 } 355 if (!okij(i1+newnlo,j1)) { 356 ok = false; 357 std::cerr<<"Start of bottom diagonal ("<<i1+newnlo<<','<<j1; 358 std::cerr<<") must be in band\n"; 359 } 360 if (newnhi >= j2-j1) { 361 ok = false; 362 std::cerr<<"new nhi ("<<newnhi<<") must be less than the new rowsize ("; 363 std::cerr<<j2-j1<<")\n"; 364 } 365 if (newnlo >= i2-i1) { 366 ok = false; 367 std::cerr<<"new nlo ("<<newnlo<<") must be less than the new colsize ("; 368 std::cerr<<i2-i1<<")\n"; 369 } 370 return ok; 371 } 372 373 template <class T> hasSubMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t istep,ptrdiff_t jstep) const374 bool ConstBandMatrixView<T,FortranStyle>::hasSubMatrix( 375 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t istep, ptrdiff_t jstep) const 376 { 377 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 378 bool ok = true; 379 if (istep == 0) { 380 ok = false; 381 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 382 } 383 if (i1 < 1 || i1 > this->colsize()) { 384 ok = false; 385 std::cerr<<"first col element ("<<i1<<") must be in 1 -- "; 386 std::cerr<<this->colsize()<<std::endl; 387 } 388 if (i2 < 1 || i2 > this->colsize()) { 389 ok = false; 390 std::cerr<<"last col element ("<<i2<<") must be in 1 -- "; 391 std::cerr<<this->colsize()<<std::endl; 392 } 393 if ((i2-i1)%istep != 0) { 394 ok = false; 395 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 396 std::cerr<<istep<<")\n"; 397 } 398 if ((i2-i1)/istep < 0) { 399 ok = false; 400 std::cerr<<"n col elements ("<<(i2-i1)/istep+1<<") must be positive\n"; 401 } 402 if (jstep == 0) { 403 ok = false; 404 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 405 } 406 if (j1 < 0 || j1 >= this->rowsize()) { 407 ok = false; 408 std::cerr<<"first row element ("<<j1<<") must be in 1 -- "; 409 std::cerr<<this->rowsize()<<std::endl; 410 } 411 if (j2 < 0 || j2 >= this->rowsize()) { 412 ok = false; 413 std::cerr<<"last row element ("<<j2<<") must be in 1 -- "; 414 std::cerr<<this->rowsize()<<std::endl; 415 } 416 if ((j2-j1)%jstep != 0) { 417 ok = false; 418 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 419 std::cerr<<jstep<<")\n"; 420 } 421 if ((j2-j1)/jstep < 0) { 422 ok = false; 423 std::cerr<<"n row elements ("<<(j2-j1)/jstep+1<<") must be positive\n"; 424 } 425 if (!this->okij(i1-1,j1-1)) { 426 ok = false; 427 std::cerr<<"Upper left corner ("<<i1<<','<<j1<<") must be in band\n"; 428 } 429 if (!this->okij(i1-1,j2-1)) { 430 ok = false; 431 std::cerr<<"Upper right corner ("<<i1<<','<<j2<<") must be in band\n"; 432 } 433 if (!this->okij(i2-1,j1-1)) { 434 ok = false; 435 std::cerr<<"Lower left corner ("<<i2<<','<<j1<<") must be in band\n"; 436 } 437 if (!this->okij(i2-1,j2-1)) { 438 ok = false; 439 std::cerr<<"Lower right corner ("<<i2<<','<<j2<<") must be in band\n"; 440 } 441 return ok; 442 } 443 444 template <class T> hasSubVector(ptrdiff_t i,ptrdiff_t j,ptrdiff_t istep,ptrdiff_t jstep,ptrdiff_t size) const445 bool ConstBandMatrixView<T,FortranStyle>::hasSubVector( 446 ptrdiff_t i, ptrdiff_t j, ptrdiff_t istep, ptrdiff_t jstep, ptrdiff_t size) const 447 { 448 if (size==0) return true; 449 bool ok = true; 450 if (istep == 0 && jstep == 0) { 451 ok = false; 452 std::cerr<<"istep ("<<istep<<") and jstep ("<<jstep; 453 std::cerr<<") can not both be 0\n"; 454 } 455 if (i < 1 || i > this->colsize()) { 456 ok = false; 457 std::cerr<<"i ("<<i<<") must be in 1 -- "<<this->colsize()<<std::endl; 458 } 459 if (j < 1 || j > this->rowsize()) { 460 ok = false; 461 std::cerr<<"j ("<<j<<") must be in 1 -- "<<this->rowsize()<<std::endl; 462 } 463 ptrdiff_t i2 = i+istep*(size-1); 464 ptrdiff_t j2 = j+jstep*(size-1); 465 if (i2 < 1 || i2 > this->colsize()) { 466 ok = false; 467 std::cerr<<"last element's i ("<<i2<<") must be in 1 -- "; 468 std::cerr<<this->colsize()<<std::endl; 469 } 470 if (j2 < 1 || j2 > this->rowsize()) { 471 ok = false; 472 std::cerr<<"last element's j ("<<j2<<") must be in 1 -- "; 473 std::cerr<<this->rowsize()<<std::endl; 474 } 475 if (!this->okij(i-1,j-1)) { 476 ok = false; 477 std::cerr<<"First element ("<<i<<','<<j<<") must be in band\n"; 478 } 479 if (!this->okij(i2-1,j2-1)) { 480 ok = false; 481 std::cerr<<"Last element ("<<i2<<','<<j2<<") must be in band\n"; 482 } 483 return ok; 484 } 485 486 template <class T> hasSubBandMatrix(ptrdiff_t i1,ptrdiff_t i2,ptrdiff_t j1,ptrdiff_t j2,ptrdiff_t newnlo,ptrdiff_t newnhi,ptrdiff_t istep,ptrdiff_t jstep) const487 bool ConstBandMatrixView<T,FortranStyle>::hasSubBandMatrix( 488 ptrdiff_t i1, ptrdiff_t i2, ptrdiff_t j1, ptrdiff_t j2, ptrdiff_t newnlo, ptrdiff_t newnhi, 489 ptrdiff_t istep, ptrdiff_t jstep) const 490 { 491 if (i1==i2 || j1==j2) return true; // no elements, so whatever... 492 bool ok = true; 493 if (istep == 0) { 494 ok = false; 495 std::cerr<<"istep ("<<istep<<") can not be 0\n"; 496 } 497 if (i1 < 1 || i1 > this->colsize()) { 498 ok = false; 499 std::cerr<<"first col element ("<<i1<<") must be in 1 -- "; 500 std::cerr<<this->colsize()<<std::endl; 501 } 502 if (i2 < 1 || i2 > this->colsize()) { 503 ok = false; 504 std::cerr<<"last col element ("<<i2<<") must be in 1 -- "; 505 std::cerr<<this->colsize()<<std::endl; 506 } 507 if ((i2-i1)%istep != 0) { 508 ok = false; 509 std::cerr<<"col range ("<<i2-i1<<") must be multiple of istep ("; 510 std::cerr<<istep<<")\n"; 511 } 512 if ((i2-i1)/istep < 0) { 513 ok = false; 514 std::cerr<<"n col elements ("<<(i2-i1)/istep+1<<") must be positive\n"; 515 } 516 if (jstep == 0) { 517 ok = false; 518 std::cerr<<"jstep ("<<jstep<<") can not be 0\n"; 519 } 520 if (j1 < 1 || j1 > this->rowsize()) { 521 ok = false; 522 std::cerr<<"first row element ("<<j1<<") must be in 1 -- "; 523 std::cerr<<this->rowsize()<<std::endl; 524 } 525 if (j2 < 1 || j2 > this->rowsize()) { 526 ok = false; 527 std::cerr<<"last row element ("<<j2<<") must be in 1 -- "; 528 std::cerr<<this->rowsize()<<std::endl; 529 } 530 if ((j2-j1)%jstep != 0) { 531 ok = false; 532 std::cerr<<"row range ("<<j2-j1<<") must be multiple of istep ("; 533 std::cerr<<jstep<<")\n"; 534 } 535 if ((j2-j1)/jstep < 0) { 536 ok = false; 537 std::cerr<<"n row elements ("<<(j2-j1)/jstep+1<<") must be positive\n"; 538 } 539 if (!this->okij(i1-1,j1-1)) { 540 ok = false; 541 std::cerr<<"Upper left corner ("<<i1<<','<<j1<<") must be in band\n"; 542 } 543 if (!this->okij(i1-1,j1-1+newnhi)) { 544 ok = false; 545 std::cerr<<"Start of top diagonal ("<<i1<<','<<j1+newnhi; 546 std::cerr<<") must be in band\n"; 547 } 548 if (!this->okij(i1-1+newnlo,j1-1)) { 549 ok = false; 550 std::cerr<<"Start of bottom diagonal ("<<i1+newnlo<<','<<j1; 551 std::cerr<<") must be in band\n"; 552 } 553 if (newnhi >= j2-j1+1) { 554 ok = false; 555 std::cerr<<"new nhi ("<<newnhi<<") must be less than the new rowsize ("; 556 std::cerr<<j2-j1+1<<")\n"; 557 } 558 if (newnlo >= i2-i1+1) { 559 ok = false; 560 std::cerr<<"new nlo ("<<newnlo<<") must be less than the new colsize ("; 561 std::cerr<<i2-i1+1<<")\n"; 562 } 563 return ok; 564 } 565 566 template <class T, int A> canLinearize() const567 bool ConstBandMatrixView<T,A>::canLinearize() const 568 { 569 if (linsize == -1) { 570 ptrdiff_t rs = this->rowsize(); 571 ptrdiff_t cs = this->colsize(); 572 ptrdiff_t lo = this->nlo(); 573 ptrdiff_t hi = this->nhi(); 574 if (rs > cs+this->nhi()) rs = cs+this->nhi(); 575 if (cs > rs+this->nlo()) cs = rs+this->nlo(); 576 577 if (rs == 0 || cs == 0) linsize = 0; 578 else if (stepi() == 1 && stepj() == (lo+hi)) 579 linsize = cs + (rs-1)*(lo+hi); 580 else if (stepj() == 1 && stepi() == (lo+hi)) 581 linsize = rs + (cs-1)*(lo+hi); 582 } 583 return linsize > 0; 584 } 585 586 template <class T, int A> canLinearize() const587 bool BandMatrixView<T,A>::canLinearize() const 588 { 589 if (linsize == -1) { 590 ptrdiff_t rs = this->rowsize(); 591 ptrdiff_t cs = this->colsize(); 592 ptrdiff_t lo = this->nlo(); 593 ptrdiff_t hi = this->nhi(); 594 if (rs > cs+this->nhi()) rs = cs+this->nhi(); 595 if (cs > rs+this->nlo()) cs = rs+this->nlo(); 596 597 if (rs == 0 || cs == 0) linsize = 0; 598 else if (stepi() == 1 && stepj() == (lo+hi)) 599 linsize = cs + (rs-1)*(lo+hi); 600 else if (stepj() == 1 && stepi() == (lo+hi)) 601 linsize = rs + (cs-1)*(lo+hi); 602 } 603 return linsize > 0; 604 } 605 QInverse() const606 template <class T> QuotXB<T,T> GenBandMatrix<T>::QInverse() const 607 { return QuotXB<T,T>(T(1),*this); } 608 609 // 610 // Norms 611 // 612 613 template <class T> det() const614 T GenBandMatrix<T>::det() const 615 { return DivHelper<T>::det(); } 616 617 template <class T> logDet(T * sign) const618 RT GenBandMatrix<T>::logDet(T* sign) const 619 { return DivHelper<T>::logDet(sign); } 620 621 template <class T> isSingular() const622 bool GenBandMatrix<T>::isSingular() const 623 { return DivHelper<T>::isSingular(); } 624 625 #ifdef INST_INT 626 template <> det() const627 int GenBandMatrix<int>::det() const 628 { return IntegerDet(*this); } 629 630 template <> det() const631 std::complex<int> GenBandMatrix<std::complex<int> >::det() const 632 { return IntegerDet(*this); } 633 634 template <> logDet(int *) const635 int GenBandMatrix<int>::logDet(int* ) const 636 { TMVAssert(TMV_FALSE); return 0; } 637 638 template <> logDet(std::complex<int> *) const639 int GenBandMatrix<std::complex<int> >::logDet(std::complex<int>* ) const 640 { TMVAssert(TMV_FALSE); return 0; } 641 642 template <> isSingular() const643 bool GenBandMatrix<int>::isSingular() const 644 { return det() == 0; } 645 646 template <> isSingular() const647 bool GenBandMatrix<std::complex<int> >::isSingular() const 648 { return det() == 0; } 649 #endif 650 651 template <class T> sumElements() const652 T GenBandMatrix<T>::sumElements() const 653 { 654 const ptrdiff_t M = colsize(); 655 const ptrdiff_t N = rowsize(); 656 T sum = 0; 657 if (M > 0 && N > 0) { 658 if (isrm()) { 659 ptrdiff_t j1=0; 660 ptrdiff_t j2=nhi()+1; 661 ptrdiff_t k=nlo(); 662 for(ptrdiff_t i=0;i<M;++i) { 663 sum += row(i,j1,j2).sumElements(); 664 if (k>0) --k; else ++j1; 665 if (j2<N) ++j2; 666 else if (j1==N) break; 667 } 668 } else if (iscm()) { 669 ptrdiff_t i1=0; 670 ptrdiff_t i2=nlo()+1; 671 ptrdiff_t k=nhi(); 672 for(ptrdiff_t j=0;j<N;++j) { 673 sum += col(j,i1,i2).sumElements(); 674 if (k>0) --k; else ++i1; 675 if (i2<M) ++i2; 676 else if (i1==M) break; 677 } 678 } else { 679 for(ptrdiff_t i=-nlo();i<=nhi();++i) sum += diag(i).sumElements(); 680 } 681 } 682 return sum; 683 } 684 685 template <class T> DoSumAbsElements(const GenBandMatrix<T> & m)686 static RT DoSumAbsElements(const GenBandMatrix<T>& m) 687 { 688 const ptrdiff_t M = m.colsize(); 689 const ptrdiff_t N = m.rowsize(); 690 RT sum = 0; 691 if (M > 0 && N > 0) { 692 if (m.isrm()) { 693 ptrdiff_t j1=0; 694 ptrdiff_t j2=m.nhi()+1; 695 ptrdiff_t k=m.nlo(); 696 for(ptrdiff_t i=0;i<M;++i) { 697 sum += m.row(i,j1,j2).sumAbsElements(); 698 if (k>0) --k; else ++j1; 699 if (j2<N) ++j2; 700 else if (j1==N) break; 701 } 702 } else if (m.iscm()) { 703 ptrdiff_t i1=0; 704 ptrdiff_t i2=m.nlo()+1; 705 ptrdiff_t k=m.nhi(); 706 for(ptrdiff_t j=0;j<N;++j) { 707 sum += m.col(j,i1,i2).sumAbsElements(); 708 if (k>0) --k; else ++i1; 709 if (i2<M) ++i2; 710 else if (i1==M) break; 711 } 712 } else { 713 for(ptrdiff_t i=-m.nlo();i<=m.nhi();++i) 714 sum += m.diag(i).sumAbsElements(); 715 } 716 } 717 return sum; 718 } 719 720 template <class T> DoSumAbs2Elements(const GenBandMatrix<T> & m)721 static RT DoSumAbs2Elements(const GenBandMatrix<T>& m) 722 { 723 const ptrdiff_t M = m.colsize(); 724 const ptrdiff_t N = m.rowsize(); 725 RT sum = 0; 726 if (M > 0 && N > 0) { 727 if (m.isrm()) { 728 ptrdiff_t j1=0; 729 ptrdiff_t j2=m.nhi()+1; 730 ptrdiff_t k=m.nlo(); 731 for(ptrdiff_t i=0;i<M;++i) { 732 sum += m.row(i,j1,j2).sumAbs2Elements(); 733 if (k>0) --k; else ++j1; 734 if (j2<N) ++j2; 735 else if (j1==N) break; 736 } 737 } else if (m.iscm()) { 738 ptrdiff_t i1=0; 739 ptrdiff_t i2=m.nlo()+1; 740 ptrdiff_t k=m.nhi(); 741 for(ptrdiff_t j=0;j<N;++j) { 742 sum += m.col(j,i1,i2).sumAbs2Elements(); 743 if (k>0) --k; else ++i1; 744 if (i2<M) ++i2; 745 else if (i1==M) break; 746 } 747 } else { 748 for(ptrdiff_t i=-m.nlo();i<=m.nhi();++i) 749 sum += m.diag(i).sumAbs2Elements(); 750 } 751 } 752 return sum; 753 } 754 755 #ifdef INST_INT DoSumAbsElements(const GenBandMatrix<std::complex<int>> &)756 static int DoSumAbsElements(const GenBandMatrix<std::complex<int> >& ) 757 { TMVAssert(TMV_FALSE); return 0; } 758 #endif 759 760 template <class T> sumAbsElements() const761 RT GenBandMatrix<T>::sumAbsElements() const 762 { return DoSumAbsElements(*this); } 763 764 template <class T> sumAbs2Elements() const765 RT GenBandMatrix<T>::sumAbs2Elements() const 766 { return DoSumAbs2Elements(*this); } 767 768 template <class T> normSq(RT scale) const769 RT GenBandMatrix<T>::normSq(RT scale) const 770 { 771 const ptrdiff_t M = colsize(); 772 const ptrdiff_t N = rowsize(); 773 RT sum = 0; 774 if (M > 0 && N > 0) { 775 if (isrm()) { 776 ptrdiff_t j1=0; 777 ptrdiff_t j2=nhi()+1; 778 ptrdiff_t k=nlo(); 779 for(ptrdiff_t i=0;i<M;++i) { 780 sum += row(i,j1,j2).normSq(scale); 781 if (k>0) --k; else ++j1; 782 if (j2<N) ++j2; 783 else if (j1==N) break; 784 } 785 } else if (iscm()) { 786 ptrdiff_t i1=0; 787 ptrdiff_t i2=nlo()+1; 788 ptrdiff_t k=nhi(); 789 for(ptrdiff_t j=0;j<N;++j) { 790 sum += col(j,i1,i2).normSq(scale); 791 if (k>0) --k; else ++i1; 792 if (i2<M) ++i2; 793 else if (i1==M) break; 794 } 795 } else { 796 for(ptrdiff_t i=-nlo();i<=nhi();++i) sum += diag(i).normSq(scale); 797 } 798 } 799 return sum; 800 } 801 802 template <class T> NonLapMaxAbsElement(const GenBandMatrix<T> & m)803 static RT NonLapMaxAbsElement(const GenBandMatrix<T>& m) 804 { 805 RT max = 0; 806 const ptrdiff_t M = m.colsize(); 807 const ptrdiff_t N = m.rowsize(); 808 if (M > 0 && N > 0) { 809 if (m.isrm()) { 810 ptrdiff_t j1=0; 811 ptrdiff_t j2=m.nhi()+1; 812 ptrdiff_t k=m.nlo(); 813 for(ptrdiff_t i=0;i<M;++i) { 814 RT temp = m.row(i,j1,j2).maxAbsElement(); 815 if (temp > max) max = temp; 816 if (k>0) --k; else ++j1; 817 if (j2<N) ++j2; 818 else if (j1==N) break; 819 } 820 } else if (m.iscm()) { 821 ptrdiff_t i1=0; 822 ptrdiff_t i2=m.nlo()+1; 823 ptrdiff_t k=m.nhi(); 824 for(ptrdiff_t j=0;j<N;++j) { 825 RT temp = m.col(j,i1,i2).maxAbsElement(); 826 if (temp > max) max = temp; 827 if (k>0) --k; else ++i1; 828 if (i2<M) ++i2; 829 else if (i1==M) break; 830 } 831 } else { 832 for(ptrdiff_t i=-m.nlo();i<=m.nhi();++i) { 833 RT temp = m.diag(i).maxAbsElement(); 834 if (temp > max) max = temp; 835 } 836 } 837 } 838 return max; 839 } 840 841 template <class T> NonLapMaxAbs2Element(const GenBandMatrix<T> & m)842 static RT NonLapMaxAbs2Element(const GenBandMatrix<T>& m) 843 { 844 RT max = 0; 845 const ptrdiff_t M = m.colsize(); 846 const ptrdiff_t N = m.rowsize(); 847 if (M > 0 && N > 0) { 848 if (m.isrm()) { 849 ptrdiff_t j1=0; 850 ptrdiff_t j2=m.nhi()+1; 851 ptrdiff_t k=m.nlo(); 852 for(ptrdiff_t i=0;i<M;++i) { 853 RT temp = m.row(i,j1,j2).maxAbs2Element(); 854 if (temp > max) max = temp; 855 if (k>0) --k; else ++j1; 856 if (j2<N) ++j2; 857 else if (j1==N) break; 858 } 859 } else if (m.iscm()) { 860 ptrdiff_t i1=0; 861 ptrdiff_t i2=m.nlo()+1; 862 ptrdiff_t k=m.nhi(); 863 for(ptrdiff_t j=0;j<N;++j) { 864 RT temp = m.col(j,i1,i2).maxAbs2Element(); 865 if (temp > max) max = temp; 866 if (k>0) --k; else ++i1; 867 if (i2<M) ++i2; 868 else if (i1==M) break; 869 } 870 } else { 871 for(ptrdiff_t i=-m.nlo();i<=m.nhi();++i) { 872 RT temp = m.diag(i).maxAbs2Element(); 873 if (temp > max) max = temp; 874 } 875 } 876 } 877 return max; 878 } 879 880 // 1-norm = max_j (sum_i |a_ij|) 881 template <class T> NonLapNorm1(const GenBandMatrix<T> & m)882 static RT NonLapNorm1(const GenBandMatrix<T>& m) 883 { 884 const ptrdiff_t M = m.colsize(); 885 const ptrdiff_t N = m.rowsize(); 886 RT max = 0; 887 if (M > 0 && N > 0) { 888 ptrdiff_t i1=0; 889 ptrdiff_t i2=m.nlo()+1; 890 ptrdiff_t k=m.nhi(); 891 for(ptrdiff_t j=0;j<N;++j) { 892 RT temp = m.col(j,i1,i2).norm1(); 893 if (temp > max) max = temp; 894 if (k>0) --k; else ++i1; 895 if (i2<M) ++i2; 896 else if (i1==M) break; 897 } 898 } 899 return max; 900 } 901 902 template <class T> NonLapNormF(const GenBandMatrix<T> & m)903 static RT NonLapNormF(const GenBandMatrix<T>& m) 904 { 905 const RT eps = TMV_Epsilon<T>(); 906 907 RT mmax = m.maxAbs2Element(); 908 if (mmax == RT(0)) return RT(0); 909 else if (TMV_Underflow(mmax * mmax)) { 910 // Then we need to rescale, since underflow has caused 911 // rounding errors. 912 // Epsilon is a pure power of 2, so no rounding errors from 913 // rescaling. 914 const RT inveps = RT(1)/eps; 915 RT scale = inveps; 916 mmax *= scale; 917 const RT eps2 = eps*eps; 918 while (mmax < eps2) { scale *= inveps; mmax *= inveps; } 919 return TMV_SQRT(m.normSq(scale))/scale; 920 } else if (RT(1) / mmax == RT(0)) { 921 // Then mmax is already inf, so no hope of making it more accurate. 922 return mmax; 923 } else if (RT(1) / (mmax*mmax) == RT(0)) { 924 // Then we have overflow, so we need to rescale: 925 const RT inveps = RT(1)/eps; 926 RT scale = eps; 927 mmax *= scale; 928 while (mmax > inveps) { scale *= eps; mmax *= eps; } 929 return TMV_SQRT(m.normSq(scale))/scale; 930 } else { 931 return TMV_SQRT(m.normSq()); 932 } 933 } 934 935 template <class T> NonLapNormInf(const GenBandMatrix<T> & m)936 static inline RT NonLapNormInf(const GenBandMatrix<T>& m) 937 { return NonLapNorm1(m.transpose()); } 938 939 #ifdef XLAP 940 template <class T> LapNorm(const char c,const GenBandMatrix<T> & m)941 static RT LapNorm(const char c, const GenBandMatrix<T>& m) 942 { 943 switch(c) { 944 case 'M' : return NonLapMaxAbsElement(m); 945 case '1' : return NonLapNorm1(m); 946 case 'F' : return NonLapNormF(m); 947 case 'I' : return NonLapNormInf(m); 948 default : TMVAssert(TMV_FALSE); 949 } 950 return RT(0); 951 } 952 #ifdef INST_DOUBLE 953 template <> LapNorm(const char c,const GenBandMatrix<double> & m)954 double LapNorm(const char c, const GenBandMatrix<double>& m) 955 { 956 TMVAssert(m.iscm() || (m.isdm() && m.nlo()==1 && m.nhi()==1)); 957 TMVAssert(m.isSquare()); 958 char cc = c; 959 int n = m.colsize(); 960 double norm; 961 if (BlasIsCM(m)) { 962 int kl = m.nlo(); 963 int ku = m.nhi(); 964 int lda = m.stepj()+1; 965 #ifndef LAPNOWORK 966 int lwork = c=='I' ? n : 0; 967 AlignedArray<double> work(lwork); 968 VectorViewOf(work.get(),lwork).setZero(); 969 #endif 970 norm = LAPNAME(dlangb) ( 971 LAPCM LAPV(cc),LAPV(n),LAPV(kl),LAPV(ku), 972 LAPP(m.cptr()-m.nhi()),LAPV(lda) LAPWK(work.get())); 973 } else { 974 norm = LAPNAME(dlangt) ( 975 LAPCM LAPV(cc),LAPV(n), 976 LAPP(m.cptr()+m.stepi()),LAPP(m.cptr()), 977 LAPP(m.cptr()+m.stepj())); 978 } 979 return norm; 980 } 981 template <> LapNorm(const char c,const GenBandMatrix<std::complex<double>> & m)982 double LapNorm(const char c, const GenBandMatrix<std::complex<double> >& m) 983 { 984 TMVAssert(m.iscm() || (m.isdm() && m.nlo()==1 && m.nhi()==1)); 985 TMVAssert(m.isSquare()); 986 char cc = c; 987 int n = m.colsize(); 988 double norm; 989 if (BlasIsCM(m)) { 990 int kl = m.nlo(); 991 int ku = m.nhi(); 992 int lda = m.stepj()+1; 993 #ifndef LAPNOWORK 994 int lwork = c=='I' ? n : 0; 995 AlignedArray<double> work(lwork); 996 VectorViewOf(work.get(),lwork).setZero(); 997 #endif 998 norm = LAPNAME(zlangb) ( 999 LAPCM LAPV(cc),LAPV(n),LAPV(kl),LAPV(ku), 1000 LAPP(m.cptr()-m.nhi()),LAPV(lda) LAPWK(work.get())); 1001 } else { 1002 norm = LAPNAME(zlangt) ( 1003 LAPCM LAPV(cc),LAPV(n), 1004 LAPP(m.cptr()+m.stepi()),LAPP(m.cptr()), 1005 LAPP(m.cptr()+m.stepj())); 1006 } 1007 return norm; 1008 } 1009 #endif 1010 #ifdef INST_FLOAT 1011 template <> LapNorm(const char c,const GenBandMatrix<float> & m)1012 float LapNorm(const char c, const GenBandMatrix<float>& m) 1013 { 1014 TMVAssert(m.iscm() || (m.isdm() && m.nlo()==1 && m.nhi()==1)); 1015 TMVAssert(m.isSquare()); 1016 char cc = c; 1017 int n = m.colsize(); 1018 float norm; 1019 if (BlasIsCM(m)) { 1020 int kl = m.nlo(); 1021 int ku = m.nhi(); 1022 int lda = m.stepj()+1; 1023 #ifndef LAPNOWORK 1024 int lwork = c=='I' ? n : 0; 1025 AlignedArray<float> work(lwork); 1026 VectorViewOf(work.get(),lwork).setZero(); 1027 #endif 1028 norm = LAPNAME(slangb) ( 1029 LAPCM LAPV(cc),LAPV(n),LAPV(kl),LAPV(ku), 1030 LAPP(m.cptr()-m.nhi()),LAPV(lda) LAPWK(work.get())); 1031 } else { 1032 norm = LAPNAME(slangt) ( 1033 LAPCM LAPV(cc),LAPV(n), 1034 LAPP(m.cptr()+m.stepi()),LAPP(m.cptr()), 1035 LAPP(m.cptr()+m.stepj())); 1036 } 1037 return norm; 1038 } 1039 template <> LapNorm(const char c,const GenBandMatrix<std::complex<float>> & m)1040 float LapNorm(const char c, const GenBandMatrix<std::complex<float> >& m) 1041 { 1042 TMVAssert(m.iscm() || (m.isdm() && m.nlo()==1 && m.nhi()==1)); 1043 TMVAssert(m.isSquare()); 1044 char cc = c; 1045 int n = m.colsize(); 1046 float norm; 1047 if (BlasIsCM(m)) { 1048 int kl = m.nlo(); 1049 int ku = m.nhi(); 1050 int lda = m.stepj()+1; 1051 #ifndef LAPNOWORK 1052 int lwork = c=='I' ? n : 0; 1053 AlignedArray<float> work(lwork); 1054 VectorViewOf(work.get(),lwork).setZero(); 1055 #endif 1056 norm = LAPNAME(clangb) ( 1057 LAPCM LAPV(cc),LAPV(n),LAPV(kl),LAPV(ku), 1058 LAPP(m.cptr()-m.nhi()),LAPV(lda) LAPWK(work.get())); 1059 } else { 1060 norm = LAPNAME(clangt) ( 1061 LAPCM LAPV(cc),LAPV(n), 1062 LAPP(m.cptr()+m.stepi()),LAPP(m.cptr()), 1063 LAPP(m.cptr()+m.stepj())); 1064 } 1065 return norm; 1066 } 1067 #endif 1068 #endif // XLAP 1069 1070 #ifdef INST_INT NonLapNormF(const GenBandMatrix<int> &)1071 static int NonLapNormF(const GenBandMatrix<int>& ) 1072 { TMVAssert(TMV_FALSE); return 0; } NonLapNormF(const GenBandMatrix<std::complex<int>> &)1073 static int NonLapNormF(const GenBandMatrix<std::complex<int> >& ) 1074 { TMVAssert(TMV_FALSE); return 0; } NonLapNorm1(const GenBandMatrix<std::complex<int>> &)1075 static int NonLapNorm1(const GenBandMatrix<std::complex<int> >& ) 1076 { TMVAssert(TMV_FALSE); return 0; } NonLapMaxAbsElement(const GenBandMatrix<std::complex<int>> &)1077 static int NonLapMaxAbsElement(const GenBandMatrix<std::complex<int> >& ) 1078 { TMVAssert(TMV_FALSE); return 0; } 1079 #endif 1080 1081 template <class T> maxAbsElement() const1082 RT GenBandMatrix<T>::maxAbsElement() const 1083 { 1084 #ifdef XLAP 1085 if (BlasIsRM(*this) && this->isSquare()) 1086 return LapNorm('M',transpose()); 1087 else if (BlasIsCM(*this) && this->isSquare()) return LapNorm('M',*this); 1088 else if (isdm() && this->isSquare() && nlo()==1 && nhi()==1) 1089 return LapNorm('M',*this); 1090 else 1091 #endif 1092 return NonLapMaxAbsElement(*this); 1093 } 1094 template <class T> maxAbs2Element() const1095 RT GenBandMatrix<T>::maxAbs2Element() const 1096 { 1097 #ifdef XLAP 1098 if (Traits<T>::iscomplex) return NonLapMaxAbs2Element(*this); 1099 else if (BlasIsRM(*this) && this->isSquare()) 1100 return LapNorm('M',transpose()); 1101 else if (BlasIsCM(*this) && this->isSquare()) return LapNorm('M',*this); 1102 else if (isdm() && this->isSquare() && nlo()==1 && nhi()==1) 1103 return LapNorm('M',*this); 1104 else 1105 #endif 1106 return NonLapMaxAbs2Element(*this); 1107 } 1108 template <class T> norm1() const1109 RT GenBandMatrix<T>::norm1() const 1110 { 1111 #ifdef XLAP 1112 if (BlasIsRM(*this) && this->isSquare()) return LapNorm('I',transpose()); 1113 else if (BlasIsCM(*this) && this->isSquare()) return LapNorm('1',*this); 1114 else if (isdm() && this->isSquare() && nlo()==1 && nhi()==1) 1115 return LapNorm('1',*this); 1116 else 1117 #endif 1118 return NonLapNorm1(*this); 1119 } 1120 template <class T> normF() const1121 RT GenBandMatrix<T>::normF() const 1122 { 1123 #ifdef XLAP 1124 if (BlasIsRM(*this) && this->isSquare()) return LapNorm('F',transpose()); 1125 else if (BlasIsCM(*this) && this->isSquare()) return LapNorm('F',*this); 1126 else if (isdm() && this->isSquare() && nlo()==1 && nhi()==1) 1127 return LapNorm('F',*this); 1128 else 1129 #endif 1130 return NonLapNormF(*this); 1131 } 1132 1133 template <class T> DoNorm2(const GenBandMatrix<T> & m)1134 static RT DoNorm2(const GenBandMatrix<T>& m) 1135 { 1136 if (m.colsize() < m.rowsize()) return DoNorm2(m.transpose()); 1137 if (m.rowsize() == 0) return RT(0); 1138 DiagMatrix<RT> S(m.rowsize()); 1139 SV_Decompose(m,S.view()); 1140 return S(0); 1141 } 1142 1143 template <class T> DoCondition(const GenBandMatrix<T> & m)1144 static RT DoCondition(const GenBandMatrix<T>& m) 1145 { 1146 if (m.colsize() < m.rowsize()) return DoCondition(m.transpose()); 1147 if (m.rowsize() == 0) return RT(1); 1148 DiagMatrix<RT> S(m.rowsize()); 1149 SV_Decompose(m,S.view()); 1150 return S(0)/S(S.size()-1); 1151 } 1152 1153 #ifdef INST_INT DoNorm2(const GenBandMatrix<int> &)1154 static int DoNorm2(const GenBandMatrix<int>& ) 1155 { TMVAssert(TMV_FALSE); return 0; } DoCondition(const GenBandMatrix<int> &)1156 static int DoCondition(const GenBandMatrix<int>& ) 1157 { TMVAssert(TMV_FALSE); return 0; } DoNorm2(const GenBandMatrix<std::complex<int>> &)1158 static int DoNorm2(const GenBandMatrix<std::complex<int> >& ) 1159 { TMVAssert(TMV_FALSE); return 0; } DoCondition(const GenBandMatrix<std::complex<int>> &)1160 static int DoCondition(const GenBandMatrix<std::complex<int> >& ) 1161 { TMVAssert(TMV_FALSE); return 0; } 1162 #endif 1163 1164 template <class T> doNorm2() const1165 RT GenBandMatrix<T>::doNorm2() const 1166 { return tmv::DoNorm2(*this); } 1167 1168 template <class T> doCondition() const1169 RT GenBandMatrix<T>::doCondition() const 1170 { return tmv::DoCondition(*this); } 1171 1172 // 1173 // Modifying Functions 1174 // 1175 1176 template <class T, int A> clip(RT thresh)1177 BandMatrixView<T,A>& BandMatrixView<T,A>::clip(RT thresh) 1178 { 1179 if (this->canLinearize()) linearView().clip(thresh); 1180 else { 1181 const ptrdiff_t M = colsize(); 1182 const ptrdiff_t N = rowsize(); 1183 if (M > 0 && N > 0) { 1184 if (isrm()) { 1185 ptrdiff_t j1=0; 1186 ptrdiff_t j2=nhi()+1; 1187 ptrdiff_t k=nlo(); 1188 for(ptrdiff_t i=0;i<M;++i) { 1189 row(i,j1,j2).clip(thresh); 1190 if (k>0) --k; else ++j1; 1191 if (j2<N) ++j2; 1192 else if (j1==N) break; 1193 } 1194 } else if (iscm()) { 1195 ptrdiff_t i1=0; 1196 ptrdiff_t i2=nlo()+1; 1197 ptrdiff_t k=nhi(); 1198 for(ptrdiff_t j=0;j<N;++j) { 1199 col(j,i1,i2).clip(thresh); 1200 if (k>0) --k; else ++i1; 1201 if (i2<M) ++i2; 1202 else if (i1==M) break; 1203 } 1204 } else { 1205 for(ptrdiff_t i=-nlo();i<=nhi();++i) diag(i).clip(thresh); 1206 } 1207 } 1208 } 1209 return *this; 1210 } 1211 1212 template <class T, int A> setZero()1213 BandMatrixView<T,A>& BandMatrixView<T,A>::setZero() 1214 { 1215 if (this->canLinearize()) linearView().setZero(); 1216 else { 1217 const ptrdiff_t M = colsize(); 1218 const ptrdiff_t N = rowsize(); 1219 if (M > 0 && N > 0) { 1220 if (isrm()) { 1221 ptrdiff_t j1=0; 1222 ptrdiff_t j2=nhi()+1; 1223 ptrdiff_t k=nlo(); 1224 for(ptrdiff_t i=0;i<M;++i) { 1225 row(i,j1,j2).setZero(); 1226 if (k>0) --k; else ++j1; 1227 if (j2<N) ++j2; 1228 else if (j1==N) break; 1229 } 1230 } else if (iscm()) { 1231 ptrdiff_t i1=0; 1232 ptrdiff_t i2=nlo()+1; 1233 ptrdiff_t k=nhi(); 1234 for(ptrdiff_t j=0;j<N;++j) { 1235 col(j,i1,i2).setZero(); 1236 if (k>0) --k; else ++i1; 1237 if (i2<M) ++i2; 1238 else if (i1==M) break; 1239 } 1240 } else { 1241 for(ptrdiff_t i=-nlo();i<=nhi();++i) diag(i).setZero(); 1242 } 1243 } 1244 } 1245 return *this; 1246 } 1247 1248 template <class T, int A> setAllTo(const T & x)1249 BandMatrixView<T,A>& BandMatrixView<T,A>::setAllTo(const T& x) 1250 { 1251 if (this->canLinearize()) linearView().setAllTo(x); 1252 else { 1253 const ptrdiff_t M = colsize(); 1254 const ptrdiff_t N = rowsize(); 1255 if (M > 0 && N > 0) { 1256 if (isrm()) { 1257 ptrdiff_t j1=0; 1258 ptrdiff_t j2=nhi()+1; 1259 ptrdiff_t k=nlo(); 1260 for(ptrdiff_t i=0;i<M;++i) { 1261 row(i,j1,j2).setAllTo(x); 1262 if (k>0) --k; else ++j1; 1263 if (j2<N) ++j2; 1264 else if (j1==N) break; 1265 } 1266 } else if (iscm()) { 1267 ptrdiff_t i1=0; 1268 ptrdiff_t i2=nlo()+1; 1269 ptrdiff_t k=nhi(); 1270 for(ptrdiff_t j=0;j<N;++j) { 1271 col(j,i1,i2).setAllTo(x); 1272 if (k>0) --k; else ++i1; 1273 if (i2<M) ++i2; 1274 else if (i1==M) break; 1275 } 1276 } else { 1277 for(ptrdiff_t i=-nlo();i<=nhi();++i) diag(i).setAllTo(x); 1278 } 1279 } 1280 } 1281 return *this; 1282 } 1283 1284 template <class T, int A> addToAll(const T & x)1285 BandMatrixView<T,A>& BandMatrixView<T,A>::addToAll(const T& x) 1286 { 1287 if (this->canLinearize()) linearView().addToAll(x); 1288 else { 1289 const ptrdiff_t M = colsize(); 1290 const ptrdiff_t N = rowsize(); 1291 if (M > 0 && N > 0) { 1292 if (isrm()) { 1293 ptrdiff_t j1=0; 1294 ptrdiff_t j2=nhi()+1; 1295 ptrdiff_t k=nlo(); 1296 for(ptrdiff_t i=0;i<M;++i) { 1297 row(i,j1,j2).addToAll(x); 1298 if (k>0) --k; else ++j1; 1299 if (j2<N) ++j2; 1300 else if (j1==N) break; 1301 } 1302 } else if (iscm()) { 1303 ptrdiff_t i1=0; 1304 ptrdiff_t i2=nlo()+1; 1305 ptrdiff_t k=nhi(); 1306 for(ptrdiff_t j=0;j<N;++j) { 1307 col(j,i1,i2).addToAll(x); 1308 if (k>0) --k; else ++i1; 1309 if (i2<M) ++i2; 1310 else if (i1==M) break; 1311 } 1312 } else { 1313 for(ptrdiff_t i=-nlo();i<=nhi();++i) diag(i).addToAll(x); 1314 } 1315 } 1316 } 1317 return *this; 1318 } 1319 1320 template <class T, int A> doTransposeSelf()1321 void BandMatrixView<T,A>::doTransposeSelf() 1322 { 1323 TMVAssert(colsize() == rowsize()); 1324 TMVAssert(nlo() == nhi()); 1325 for(ptrdiff_t i=1;i<=nhi();++i) Swap(diag(-i),diag(i)); 1326 } 1327 1328 template <class T, int A> conjugateSelf()1329 BandMatrixView<T,A>& BandMatrixView<T,A>::conjugateSelf() 1330 { 1331 if (this->canLinearize()) linearView().conjugateSelf(); 1332 else { 1333 const ptrdiff_t M = colsize(); 1334 const ptrdiff_t N = rowsize(); 1335 if (isComplex(T()) && M > 0 && N > 0) { 1336 if (isrm()) { 1337 ptrdiff_t j1=0; 1338 ptrdiff_t j2=nhi()+1; 1339 ptrdiff_t k=nlo(); 1340 for(ptrdiff_t i=0;i<M;++i) { 1341 row(i,j1,j2).conjugateSelf(); 1342 if (k>0) --k; else ++j1; 1343 if (j2<N) ++j2; 1344 else if (j1==N) break; 1345 } 1346 } else if (iscm()) { 1347 ptrdiff_t i1=0; 1348 ptrdiff_t i2=nlo()+1; 1349 ptrdiff_t k=nhi(); 1350 for(ptrdiff_t j=0;j<N;++j) { 1351 col(j,i1,i2).conjugateSelf(); 1352 if (k>0) --k; else ++i1; 1353 if (i2<M) ++i2; 1354 else if (i1==M) break; 1355 } 1356 } else { 1357 for(ptrdiff_t i=-nlo();i<=nhi();++i) diag(i).conjugateSelf(); 1358 } 1359 } 1360 } 1361 return *this; 1362 } 1363 1364 1365 // 1366 // Special Constructors 1367 // 1368 1369 template <class T> UpperBiDiagMatrix(const GenVector<T> & v1,const GenVector<T> & v2)1370 BandMatrix<T,DiagMajor> UpperBiDiagMatrix( 1371 const GenVector<T>& v1, const GenVector<T>& v2) 1372 { 1373 if (v1.size() == v2.size()) { 1374 BandMatrix<T,DiagMajor> temp(v1.size(),v1.size()+1,0,1); 1375 temp.diag(0) = v1; 1376 temp.diag(1) = v2; 1377 return temp; 1378 } else { 1379 TMVAssert2(v2.size() == v1.size()-1); 1380 BandMatrix<T,DiagMajor> temp(v1.size(),v1.size(),0,1); 1381 temp.diag(0) = v1; 1382 temp.diag(1) = v2; 1383 return temp; 1384 } 1385 } 1386 1387 template <class T> LowerBiDiagMatrix(const GenVector<T> & v1,const GenVector<T> & v2)1388 BandMatrix<T,DiagMajor> LowerBiDiagMatrix( 1389 const GenVector<T>& v1, const GenVector<T>& v2) 1390 { 1391 if (v1.size() == v2.size()) { 1392 BandMatrix<T,DiagMajor> temp(v2.size()+1,v2.size(),1,0); 1393 temp.diag(-1) = v1; 1394 temp.diag(0) = v2; 1395 return temp; 1396 } else { 1397 TMVAssert2(v1.size() == v2.size()-1); 1398 BandMatrix<T,DiagMajor> temp(v2.size(),v2.size(),1,0); 1399 temp.diag(-1) = v1; 1400 temp.diag(0) = v2; 1401 return temp; 1402 } 1403 } 1404 1405 template <class T> TriDiagMatrix(const GenVector<T> & v1,const GenVector<T> & v2,const GenVector<T> & v3)1406 BandMatrix<T,DiagMajor> TriDiagMatrix( 1407 const GenVector<T>& v1, const GenVector<T>& v2, 1408 const GenVector<T>& v3) 1409 { 1410 if (v1.size() == v2.size()) { 1411 TMVAssert2(v3.size() == v2.size()-1); 1412 BandMatrix<T,DiagMajor> temp(v2.size()+1,v2.size(),1,1); 1413 temp.diag(-1) = v1; 1414 temp.diag(0) = v2; 1415 temp.diag(1) = v3; 1416 return temp; 1417 } else if (v2.size() == v3.size()) { 1418 TMVAssert2(v1.size() == v2.size()-1); 1419 BandMatrix<T,DiagMajor> temp(v2.size(),v2.size()+1,1,1); 1420 temp.diag(-1) = v1; 1421 temp.diag(0) = v2; 1422 temp.diag(1) = v3; 1423 return temp; 1424 } else { 1425 TMVAssert2(v1.size() == v2.size()-1); 1426 TMVAssert2(v3.size() == v2.size()-1); 1427 BandMatrix<T,DiagMajor> temp(v2.size(),v2.size(),1,1); 1428 temp.diag(-1) = v1; 1429 temp.diag(0) = v2; 1430 temp.diag(1) = v3; 1431 return temp; 1432 } 1433 } 1434 1435 1436 // 1437 // Swap 1438 // 1439 1440 template <class T> Swap(BandMatrixView<T> m1,BandMatrixView<T> m2)1441 void Swap(BandMatrixView<T> m1, BandMatrixView<T> m2) 1442 { 1443 TMVAssert2(m1.colsize() == m2.colsize()); 1444 TMVAssert2(m1.rowsize() == m2.rowsize()); 1445 TMVAssert2(m1.nlo() == m2.nlo()); 1446 TMVAssert2(m1.nhi() == m2.nhi()); 1447 const ptrdiff_t M = m1.colsize(); 1448 const ptrdiff_t N = m1.rowsize(); 1449 if (M > 0 && N > 0) { 1450 if (m1.isrm() && m2.isrm()) { 1451 ptrdiff_t j1=0; 1452 ptrdiff_t j2=m1.nhi()+1; 1453 ptrdiff_t k=m1.nlo(); 1454 for(ptrdiff_t i=0;i<M;++i) { 1455 Swap(m1.row(i,j1,j2),m2.row(i,j1,j2)); 1456 if (k>0) --k; else ++j1; 1457 if (j2<N) ++j2; 1458 else if (j1==N) break; 1459 } 1460 } else if (m1.iscm() && m2.iscm()) { 1461 ptrdiff_t i1=0; 1462 ptrdiff_t i2=m1.nlo()+1; 1463 ptrdiff_t k=m1.nhi(); 1464 for(ptrdiff_t j=0;j<N;++j) { 1465 Swap(m1.col(j,i1,i2),m2.col(j,i1,i2)); 1466 if (k>0) --k; else ++i1; 1467 if (i2<M) ++i2; 1468 else if (i1==M) break; 1469 } 1470 } else { 1471 for(ptrdiff_t i=-m1.nlo();i<=m1.nhi();++i) 1472 Swap(m1.diag(i),m2.diag(i)); 1473 } 1474 } 1475 } 1476 1477 // 1478 // m1 == m2 1479 // 1480 1481 template <class T1, class T2> operator ==(const GenBandMatrix<T1> & m1,const GenBandMatrix<T2> & m2)1482 bool operator==(const GenBandMatrix<T1>& m1, const GenBandMatrix<T2>& m2) 1483 { 1484 if (m1.colsize() != m2.colsize()) return false; 1485 else if (m1.rowsize() != m2.rowsize()) return false; 1486 else if (m1.isSameAs(m2)) return true; 1487 else { 1488 ptrdiff_t lo = TMV_MIN(m1.nlo(),m2.nlo()); 1489 ptrdiff_t hi = TMV_MIN(m1.nhi(),m2.nhi()); 1490 for(ptrdiff_t i=-lo;i<=hi;++i) 1491 if (m1.diag(i) != m2.diag(i)) return false; 1492 for(ptrdiff_t i=-m1.nlo();i<-lo;++i) 1493 if (m1.diag(i).maxAbs2Element() != T1(0)) return false; 1494 for(ptrdiff_t i=-m2.nlo();i<-lo;++i) 1495 if (m2.diag(i).maxAbs2Element() != T2(0)) return false; 1496 for(ptrdiff_t i=hi+1;i<=m1.nhi();++i) 1497 if (m1.diag(i).maxAbs2Element() != T1(0)) return false; 1498 for(ptrdiff_t i=hi+1;i<=m2.nhi();++i) 1499 if (m2.diag(i).maxAbs2Element() != T2(0)) return false; 1500 return true; 1501 } 1502 } 1503 1504 template <class T1, class T2> operator ==(const GenBandMatrix<T1> & m1,const GenMatrix<T2> & m2)1505 bool operator==(const GenBandMatrix<T1>& m1, const GenMatrix<T2>& m2) 1506 { 1507 if (m1.colsize() != m2.colsize()) return false; 1508 else if (m1.rowsize() != m2.rowsize()) return false; 1509 else { 1510 ConstBandMatrixView<T2> m2b = 1511 BandMatrixViewOf(m2,m2.colsize()-1,m2.rowsize()-1); 1512 if ( m1.diagRange(-m1.nlo(),m1.nhi()+1) != 1513 m2b.diagRange(-m1.nlo(),m1.nhi()+1)) 1514 return false; 1515 if ( m1.nhi()+1 < m1.rowsize() && 1516 m2b.diagRange(m1.nhi()+1,m1.rowsize()).maxAbs2Element() != T2(0)) 1517 return false; 1518 if ( m1.nlo()+1 < m1.colsize() && 1519 m2b.diagRange(-m1.colsize()+1,-m1.nlo()).maxAbs2Element() != T2(0)) 1520 return false; 1521 return true; 1522 } 1523 } 1524 1525 1526 // 1527 // I/O 1528 // 1529 1530 template <class T> write(const TMV_Writer & writer) const1531 void GenBandMatrix<T>::write(const TMV_Writer& writer) const 1532 { 1533 const ptrdiff_t M = colsize(); 1534 const ptrdiff_t N = rowsize(); 1535 ptrdiff_t j1=0; 1536 ptrdiff_t j2=nhi()+1; 1537 1538 writer.begin(); 1539 writer.writeCode("B"); 1540 writer.writeSize(M); 1541 writer.writeSize(N); 1542 writer.writeFullSize(nlo()); 1543 writer.writeFullSize(nhi()); 1544 writer.writeStart(); 1545 1546 for(ptrdiff_t i=0;i<M;++i) { 1547 writer.writeLParen(); 1548 if (!writer.isCompact()) { 1549 for(ptrdiff_t j=0;j<j1;++j) { 1550 writer.writeValue(T(0)); 1551 if (j<N-1) writer.writeSpace(); 1552 } 1553 } 1554 for(ptrdiff_t j=j1;j<j2;++j) { 1555 if (j > j1) writer.writeSpace(); 1556 writer.writeValue(cref(i,j)); 1557 } 1558 if (!writer.isCompact()) { 1559 for(ptrdiff_t j=j2;j<N;++j) { 1560 writer.writeSpace(); 1561 writer.writeValue(T(0)); 1562 } 1563 } 1564 writer.writeRParen(); 1565 if (i < M-1) writer.writeRowEnd(); 1566 if (j2 < N) ++j2; 1567 if (i >= nlo() && j1 < N) ++j1; 1568 } 1569 writer.writeFinal(); 1570 writer.end(); 1571 } 1572 1573 #ifndef NOTHROW 1574 template <class T> 1575 class BandMatrixReadError : public ReadError 1576 { 1577 public : 1578 BandMatrix<T> m; 1579 ptrdiff_t i,j; 1580 std::string exp,got; 1581 ptrdiff_t cs,rs; 1582 ptrdiff_t lo,hi; 1583 T v1; 1584 bool is, iseof, isbad; 1585 BandMatrixReadError(std::istream & _is)1586 BandMatrixReadError(std::istream& _is) throw() : 1587 ReadError("BandMatrix."), 1588 i(0), j(0), cs(0), rs(0), lo(0), hi(0), v1(0), 1589 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} BandMatrixReadError(std::istream & _is,const std::string & _e,const std::string & _g)1590 BandMatrixReadError( 1591 std::istream& _is, 1592 const std::string& _e, const std::string& _g) throw() : 1593 ReadError("BandMatrix."), 1594 i(0), j(0), exp(_e), got(_g), cs(0), rs(0), lo(0), hi(0), v1(0), 1595 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 1596 BandMatrixReadError(const GenBandMatrix<T> & _m,std::istream & _is,ptrdiff_t _cs,ptrdiff_t _rs,ptrdiff_t _lo,ptrdiff_t _hi)1597 BandMatrixReadError( 1598 const GenBandMatrix<T>& _m, std::istream& _is, 1599 ptrdiff_t _cs, ptrdiff_t _rs, ptrdiff_t _lo, ptrdiff_t _hi) throw() : 1600 ReadError("BandMatrix."), 1601 m(_m), i(0), j(0), cs(_cs), rs(_rs), lo(_lo), hi(_hi), v1(0), 1602 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} BandMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenBandMatrix<T> & _m,std::istream & _is,const std::string & _e,const std::string & _g)1603 BandMatrixReadError( 1604 ptrdiff_t _i, ptrdiff_t _j, const GenBandMatrix<T>& _m, 1605 std::istream& _is, 1606 const std::string& _e, const std::string& _g) throw() : 1607 ReadError("BandMatrix."), 1608 m(_m), i(_i), j(_j), exp(_e), got(_g), 1609 cs(m.colsize()), rs(m.rowsize()), lo(m.nlo()), hi(m.nhi()), v1(0), 1610 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} BandMatrixReadError(ptrdiff_t _i,ptrdiff_t _j,const GenBandMatrix<T> & _m,std::istream & _is,T _v1=0)1611 BandMatrixReadError( 1612 ptrdiff_t _i, ptrdiff_t _j, const GenBandMatrix<T>& _m, 1613 std::istream& _is, T _v1=0) throw() : 1614 ReadError("BandMatrix."), 1615 m(_m), i(_i), j(_j), 1616 cs(m.colsize()), rs(m.rowsize()), lo(m.nlo()), hi(m.nhi()), v1(_v1), 1617 is(_is), iseof(_is.eof()), isbad(_is.bad()) {} 1618 BandMatrixReadError(const BandMatrixReadError<T> & rhs)1619 BandMatrixReadError(const BandMatrixReadError<T>& rhs) : 1620 m(rhs.m), i(rhs.i), j(rhs.j), exp(rhs.exp), got(rhs.got), 1621 cs(rhs.cs), rs(rhs.rs), lo(rhs.lo), hi(rhs.hi), v1(rhs.v1), 1622 is(rhs.is), iseof(rhs.iseof), isbad(rhs.isbad) {} ~BandMatrixReadError()1623 ~BandMatrixReadError() throw() {} 1624 write(std::ostream & os) const1625 void write(std::ostream& os) const throw() 1626 { 1627 os<<"TMV Read Error: Reading istream input for BandMatrix\n"; 1628 if (exp != got) { 1629 os<<"Wrong format: expected '"<<exp<<"', got '"<<got<<"'.\n"; 1630 } 1631 if (cs != m.colsize()) { 1632 os<<"Wrong colsize: expected "<<m.colsize()<< 1633 ", got "<<cs<<".\n"; 1634 } 1635 if (rs != m.rowsize()) { 1636 os<<"Wrong rowsize: expected "<<m.rowsize()<< 1637 ", got "<<rs<<".\n"; 1638 } 1639 if (lo != m.nlo()) { 1640 os<<"Wrong nlo: expected "<<m.nlo()<<", got "<<lo<<".\n"; 1641 } 1642 if (hi != m.nhi()) { 1643 os<<"Wrong nhi: expected "<<m.nhi()<<", got "<<hi<<".\n"; 1644 } 1645 if (!is) { 1646 if (iseof) { 1647 os<<"Input stream reached end-of-file prematurely.\n"; 1648 } else if (isbad) { 1649 os<<"Input stream is corrupted.\n"; 1650 } else { 1651 os<<"Input stream cannot read next character.\n"; 1652 } 1653 } 1654 if (v1 != T(0)) { 1655 os<<"Invalid input. Expected 0, got "<<v1<<".\n"; 1656 } 1657 if (m.colsize() > 0 || m.rowsize() > 0) { 1658 os<<"The portion of the BandMatrix which was successfully " 1659 "read is: \n"; 1660 const ptrdiff_t N = m.rowsize(); 1661 for(ptrdiff_t ii=0;ii<i;++ii) { 1662 os<<"( "; 1663 for(ptrdiff_t jj=0;jj<N;++jj) os<<' '<<m.cref(ii,jj)<<' '; 1664 os<<" )\n"; 1665 1666 } 1667 os<<"( "; 1668 for(ptrdiff_t jj=0;jj<j;++jj) os<<' '<<m.cref(i,jj)<<' '; 1669 os<<" )\n"; 1670 } 1671 } 1672 }; 1673 #endif 1674 1675 template <class T> FinishRead(const TMV_Reader & reader,BandMatrixView<T> m)1676 static void FinishRead(const TMV_Reader& reader, BandMatrixView<T> m) 1677 { 1678 const ptrdiff_t M = m.colsize(); 1679 const ptrdiff_t N = m.rowsize(); 1680 std::string exp,got; 1681 T temp; 1682 ptrdiff_t j1=0; 1683 ptrdiff_t j2=m.nhi()+1; 1684 1685 if (!reader.readStart(exp,got)) { 1686 #ifdef NOTHROW 1687 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1688 exit(1); 1689 #else 1690 throw BandMatrixReadError<T>(0,0,m,reader.getis(),exp,got); 1691 #endif 1692 } 1693 for(ptrdiff_t i=0;i<M;++i) { 1694 if (!reader.readLParen(exp,got)) { 1695 #ifdef NOTHROW 1696 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1697 exit(1); 1698 #else 1699 throw BandMatrixReadError<T>(i,0,m,reader.getis(),exp,got); 1700 #endif 1701 } 1702 if (!reader.isCompact()) { 1703 for(ptrdiff_t j=0;j<j1;++j) { 1704 if (!reader.readValue(temp)) { 1705 #ifdef NOTHROW 1706 std::cerr<<"BandMatrix Read Error: reading value\n"; 1707 exit(1); 1708 #else 1709 throw BandMatrixReadError<T>(i,j,m,reader.getis()); 1710 #endif 1711 } 1712 if (temp != T(0)) { 1713 #ifdef NOTHROW 1714 std::cerr<<"BandMatrix Read Error: "<<temp<<" != 0\n"; 1715 exit(1); 1716 #else 1717 throw BandMatrixReadError<T>(i,j,m,reader.getis(),temp); 1718 #endif 1719 } 1720 if (!reader.readSpace(exp,got)) { 1721 #ifdef NOTHROW 1722 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1723 exit(1); 1724 #else 1725 throw BandMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1726 #endif 1727 } 1728 } 1729 } 1730 for(ptrdiff_t j=j1;j<j2;++j) { 1731 if (j>j1 && !reader.readSpace(exp,got)) { 1732 #ifdef NOTHROW 1733 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1734 exit(1); 1735 #else 1736 throw BandMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1737 #endif 1738 } 1739 if (!reader.readValue(temp)) { 1740 #ifdef NOTHROW 1741 std::cerr<<"BandMatrix Read Error: reading value\n"; 1742 exit(1); 1743 #else 1744 throw BandMatrixReadError<T>(i,j,m,reader.getis()); 1745 #endif 1746 } 1747 m.ref(i,j) = temp; 1748 } 1749 if (!reader.isCompact()) { 1750 for(ptrdiff_t j=j2;j<N;++j) { 1751 if (!reader.readSpace(exp,got)) { 1752 #ifdef NOTHROW 1753 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1754 exit(1); 1755 #else 1756 throw BandMatrixReadError<T>(i,j,m,reader.getis(),exp,got); 1757 #endif 1758 } 1759 if (!reader.readValue(temp)) { 1760 #ifdef NOTHROW 1761 std::cerr<<"BandMatrix Read Error: reading value\n"; 1762 exit(1); 1763 #else 1764 throw BandMatrixReadError<T>(i,j,m,reader.getis()); 1765 #endif 1766 } 1767 if (temp != T(0)) { 1768 #ifdef NOTHROW 1769 std::cerr<<"BandMatrix Read Error: "<<temp<<" != 0\n"; 1770 exit(1); 1771 #else 1772 throw BandMatrixReadError<T>(i,j,m,reader.getis(),temp); 1773 #endif 1774 } 1775 } 1776 } 1777 if (!reader.readRParen(exp,got)) { 1778 #ifdef NOTHROW 1779 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1780 exit(1); 1781 #else 1782 throw BandMatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1783 #endif 1784 } 1785 if (i < M-1 && !reader.readRowEnd(exp,got)) { 1786 #ifdef NOTHROW 1787 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1788 exit(1); 1789 #else 1790 throw BandMatrixReadError<T>(i,N,m,reader.getis(),exp,got); 1791 #endif 1792 } 1793 if (j2 < N) ++j2; 1794 if (i >= m.nlo() && j1 < N) ++j1; 1795 } 1796 if (!reader.readFinal(exp,got)) { 1797 #ifdef NOTHROW 1798 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1799 exit(1); 1800 #else 1801 throw BandMatrixReadError<T>(N,0,m,reader.getis(),exp,got); 1802 #endif 1803 } 1804 } 1805 1806 template <class T, int A> read(const TMV_Reader & reader)1807 void BandMatrix<T,A>::read(const TMV_Reader& reader) 1808 { 1809 std::string exp,got; 1810 if (!reader.readCode("B",exp,got)) { 1811 #ifdef NOTHROW 1812 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1813 exit(1); 1814 #else 1815 throw BandMatrixReadError<T>(reader.getis(),exp,got); 1816 #endif 1817 } 1818 ptrdiff_t cs=colsize(), rs=rowsize(), lo=nlo(), hi=nhi(); 1819 if (!reader.readSize(cs,exp,got) || 1820 !reader.readSize(rs,exp,got) || 1821 !reader.readFullSize(lo,exp,got) || 1822 !reader.readFullSize(hi,exp,got)) { 1823 #ifdef NOTHROW 1824 std::cerr<<"BandMatrix Read Error: reading size\n"; 1825 exit(1); 1826 #else 1827 throw BandMatrixReadError<T>(reader.getis(),exp,got); 1828 #endif 1829 } 1830 if (cs != colsize() || rs != rowsize() || lo != nlo() || hi != nhi()) { 1831 resize(cs,rs,lo,hi); 1832 } 1833 BandMatrixView<T> v = view(); 1834 FinishRead(reader,v); 1835 } 1836 1837 template <class T, int A> read(const TMV_Reader & reader)1838 void BandMatrixView<T,A>::read(const TMV_Reader& reader) 1839 { 1840 std::string exp,got; 1841 if (!reader.readCode("B",exp,got)) { 1842 #ifdef NOTHROW 1843 std::cerr<<"BandMatrix Read Error: "<<got<<" != "<<exp<<std::endl; 1844 exit(1); 1845 #else 1846 throw BandMatrixReadError<T>(reader.getis(),exp,got); 1847 #endif 1848 } 1849 ptrdiff_t cs=colsize(), rs=rowsize(), lo=nlo(), hi=nhi(); 1850 if (!reader.readSize(cs,exp,got) || 1851 !reader.readSize(rs,exp,got) || 1852 !reader.readFullSize(lo,exp,got) || 1853 !reader.readFullSize(hi,exp,got)) { 1854 #ifdef NOTHROW 1855 std::cerr<<"BandMatrix Read Error: reading size\n"; 1856 exit(1); 1857 #else 1858 throw BandMatrixReadError<T>(reader.getis(),exp,got); 1859 #endif 1860 } 1861 if (cs != colsize() || rs != rowsize() || lo != nlo() || hi != nhi()) { 1862 #ifdef NOTHROW 1863 std::cerr<<"BandMatrix Read Error: Wrong size\n"; 1864 exit(1); 1865 #else 1866 throw BandMatrixReadError<T>(*this,reader.getis(),cs,rs,lo,hi); 1867 #endif 1868 } 1869 FinishRead(reader,*this); 1870 } 1871 1872 #undef RT 1873 1874 #define InstFile "TMV_BandMatrix.inst" 1875 #include "TMV_Inst.h" 1876 #undef InstFile 1877 1878 } // namespace tmv 1879 1880 1881