1/* 2 * Copyright (C) 2004 Pascal Giorgi, Clément Pernet 3 * 2013, 2014 the LinBox group 4 * 5 * Written by : 6 * Pascal Giorgi <pascal.giorgi@ens-lyon.fr> 7 * Clément Pernet <clement.pernet@imag.fr> 8 * Brice Boyer (briceboyer) <boyer.brice@gmail.com> 9 * 10 * ========LICENCE======== 11 * This file is part of the library LinBox. 12 * 13 * LinBox is free software: you can redistribute it and/or modify 14 * it under the terms of the GNU Lesser General Public 15 * License as published by the Free Software Foundation; either 16 * version 2.1 of the License, or (at your option) any later version. 17 * 18 * This library is distributed in the hope that it will be useful, 19 * but WITHOUT ANY WARRANTY; without even the implied warranty of 20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 21 * Lesser General Public License for more details. 22 * 23 * You should have received a copy of the GNU Lesser General Public 24 * License along with this library; if not, write to the Free Software 25 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 26 * ========LICENCE======== 27 */ 28 29/*!@internal 30 * @file matrix/densematrix/blas-matrix.inl 31 * @ingroup densematrix 32 * A \c BlasMatrix<\c _Field > represents a matrix as an array of 33 * <code>_Field</code>s. 34 */ 35 36#ifndef __LINBOX_densematrix_blas_submatrix_INL 37#define __LINBOX_densematrix_blas_submatrix_INL 38 39///////////////// 40// PRIVATE // 41///////////////// 42 43namespace LinBox 44{ 45} 46 47////////////////// 48// CONSTRUCTORS // 49////////////////// 50 51namespace LinBox 52{ 53#if 0 54 template < class _Matrix > 55 BlasSubmatrix<_Matrix>::BlasSubmatrix () : 56 _Mat(NULL),_row(0),_col(0),_r0(0),_c0(0),_stride(0),_off(0) 57 { 58 } 59#endif 60 61 template < class _Matrix > 62 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::matrixType &Mat, 63 size_t row, 64 size_t col, 65 size_t Rowdim, 66 size_t Coldim) : 67 _Mat (Mat), 68 _row (Rowdim), _col(Coldim), 69 _r0(row),_c0(col), 70 _stride(Mat.coldim()),_off(row*_stride+col) 71 ,_AD(Mat.field()) 72 ,_VD(Mat.field()) 73 { 74 // std::cout << "sub cstor 0 called" << std::endl; 75 linbox_check ( _r0 <= Mat.rowdim() ); // allow for NULL matrix 76 linbox_check ( _c0 <= Mat.coldim() ); 77 // linbox_check ( _row <= Mat.rowdim() ); 78 linbox_check ( _off <= Mat.rowdim()*Mat.coldim() ); 79 linbox_check ( _col <= Mat.coldim() ); 80 } 81 82 template < class _Matrix > 83 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::constMatrixType &Mat, 84 size_t row, 85 size_t col, 86 size_t Rowdim, 87 size_t Coldim) : 88 _Mat (Mat), 89 _row (Rowdim), _col(Coldim), 90 _r0(row),_c0(col), 91 _stride(Mat.coldim()),_off(row*_stride+col) 92 ,_AD(Mat.field()) 93 ,_VD(Mat.field()) 94 { 95 // std::cout << "sub const cstor 0 called" << std::endl; 96 linbox_check ( _r0 <= Mat.rowdim() ); // allow for NULL matrix 97 linbox_check ( _c0 <= Mat.coldim() ); 98 // linbox_check ( _row <= Mat.rowdim() ); 99 linbox_check ( _off <= Mat.rowdim()*Mat.coldim() ); 100 linbox_check ( _col <= Mat.coldim() ); 101 } 102 103 104 template < class _Matrix > 105 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::matrixType &Mat) : 106 _Mat (Mat), 107 _row(Mat.rowdim()), _col(Mat.coldim()), 108 _r0(0), _c0(0), 109 _stride(Mat.coldim()),_off(0) 110 ,_AD(Mat.field()) 111 ,_VD(Mat.field()) 112 { 113 // std::cout << "sub cstor 2 called" << std::endl; 114 } 115 116 template < class _Matrix > 117 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::constMatrixType &Mat) : 118 _Mat (Mat), 119 _row(Mat.rowdim()), _col(Mat.coldim()), 120 _r0(0), _c0(0), 121 _stride(Mat.coldim()),_off(0) 122 ,_AD(Mat.field()) 123 ,_VD(Mat.field()) 124 { 125 // std::cout << "sub const cstor 2 called" << std::endl; 126 } 127 128 template < class _Matrix > 129 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::constSubMatrixType &SM, 130 size_t row, 131 size_t col, 132 size_t Rowdim, 133 size_t Coldim) : 134 _Mat (SM._Mat), 135 _row ( Rowdim ), _col ( Coldim ) , 136 _r0 ( SM._r0 + row ), _c0 ( SM._c0 + col ), 137 _stride(SM._stride),_off(SM._off+(row*_stride+col)) 138 ,_AD(SM.field()) 139 ,_VD(SM.field()) 140 { 141 // std::cout << "sub const cstor 3 called" << std::endl; 142 linbox_check ( _r0 <= SM._Mat.rowdim() ); // allow for NULL matrix 143 linbox_check ( _c0 <= SM._stride ); 144 linbox_check ( _row <= SM._Mat.rowdim() ); 145 linbox_check ( _col <= SM._stride ); 146 linbox_check ( _off <= SM._Mat.rowdim()*(SM._Mat.coldim()+1) ); 147 } 148 149 150 template < class _Matrix > 151 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::subMatrixType &SM, 152 size_t row, 153 size_t col, 154 size_t Rowdim, 155 size_t Coldim) : 156 _Mat (SM._Mat), 157 _row ( Rowdim ), _col ( Coldim ) , 158 _r0 ( SM._r0 + row ), _c0 ( SM._c0 + col ), 159 _stride(SM._stride),_off(SM._off+(row*_stride+col)) 160 ,_AD(SM.field()) 161 ,_VD(SM.field()) 162 { 163 // std::cout << "sub cstor 3 called" << std::endl; 164 linbox_check ( _r0 <= SM._Mat.rowdim() ); // allow for NULL matrix 165 linbox_check ( _c0 <= SM._stride ); 166 linbox_check ( _row <= SM._Mat.rowdim() ); 167 linbox_check ( _col <= SM._stride ); 168 linbox_check ( _off <= SM._Mat.rowdim()*(SM._Mat.coldim()+1) ); 169 } 170 171 template < class _Matrix > 172 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::subMatrixType &SM) : 173 _Mat (SM._Mat), 174 _row ( SM._row ), _col ( SM._col ) , 175 _r0 ( SM._r0 ), _c0 (SM._c0 ), 176 _stride(SM._stride),_off(SM._off) 177 ,_AD(SM.field()) 178 ,_VD(SM.field()) 179 { 180 // std::cout << "sub cstor 4 called" << std::endl; 181 linbox_check ( _r0 <= SM._Mat.rowdim() ); // allow for NULL matrix 182 linbox_check ( _c0 <= SM._stride ); 183 linbox_check ( _row <= SM._Mat.rowdim() ); 184 linbox_check ( _col <= SM._stride ); 185 linbox_check ( _off <= SM._Mat.rowdim()*(SM._Mat.coldim()+1) ); 186 } 187 188 template < class _Matrix > 189 BlasSubmatrix<_Matrix>::BlasSubmatrix (typename BlasSubmatrix<_Matrix>::constSubMatrixType &SM) : 190 _Mat (SM._Mat), 191 _row ( SM._row ), _col ( SM._col ) , 192 _r0 ( SM._r0 ), _c0 (SM._c0 ), 193 _stride(SM._stride),_off(SM._off) 194 ,_AD(SM.field()) 195 ,_VD(SM.field()) 196 { 197 // std::cout << "sub const cstor 4 called" << std::endl; 198 linbox_check ( _r0 <= SM._Mat.rowdim() ); // allow for NULL matrix 199 linbox_check ( _c0 <= SM._stride ); 200 linbox_check ( _row <= SM._Mat.rowdim() ); 201 linbox_check ( _col <= SM._stride ); 202 linbox_check ( _off <= SM._Mat.rowdim()*(SM._Mat.coldim()+1) ); 203 } 204 205 // shallow copy 206 template < class _Matrix > 207 BlasSubmatrix<_Matrix>& BlasSubmatrix<_Matrix>::operator=(const BlasSubmatrix<_Matrix> &SM) 208 { 209 if ( &SM == this) 210 return *this ; 211 212 _Mat = SM._Mat ; 213 _r0 = SM._r0 ; 214 _row = SM._row; 215 _c0 = SM._c0 ; 216 _col = SM._col; 217 _stride = SM._stride; 218 _off = SM._off ; 219 _AD = SM._AD ; 220 _VD = SM._VD ; 221 222 return *this; 223 } 224 225 // function for repurposing Submatrices. 226 template < class _Matrix > 227 BlasSubmatrix<_Matrix>& BlasSubmatrix<_Matrix>::submatrix(typename BlasSubmatrix<_Matrix>::constSelf_t &SM, 228 size_t row, size_t col, size_t Rowdim, size_t Coldim) 229 { 230 _Mat = SM._Mat; 231 _row = Rowdim; _col = Coldim; 232 _r0 = SM._r0 + row; _c0 = SM._c0 + col; 233 _stride = SM._stride; _off = SM._off+(row*_stride+col); 234 _AD = SM._AD; 235 _VD = SM._VD; 236 237 return *this; 238 } 239 240 // deep copy 241 template < class _Matrix > 242 template<class Matrix> 243 BlasSubmatrix<_Matrix>& BlasSubmatrix<_Matrix>::copy( const Matrix & B) 244 { 245 for (size_t i = 0 ; i < rowdim() ; ++i) 246 for (size_t j = 0 ; j < coldim() ; ++j) { 247 setEntry(i,j,B.getEntry(i,j)); 248 } 249 return *this; 250 } 251 252 template < class _Matrix > 253 BlasSubmatrix<_Matrix> &BlasSubmatrix<_Matrix>::swap( typename BlasSubmatrix<_Matrix>::Self_t & B) 254 { 255 Element temp; _Mat.field().init(temp); 256 Element hold; _Mat.field().init(hold); 257 for (size_t i = 0 ; i < rowdim() ; ++i) 258 for (size_t j = 0 ; j < coldim() ; ++j) { 259 getEntry(hold,i,j); 260 setEntry(i,j,B.getEntry(temp,i,j)); 261 B.setEntry(i,j,hold); 262 } 263 return *this; 264 } 265 266 template < class _Matrix > 267 BlasSubmatrix<_Matrix> &BlasSubmatrix<_Matrix>::zero() 268 { 269 for (size_t i = 0 ; i < rowdim() ; ++i) 270 for (size_t j = 0 ; j < coldim() ; ++j) { 271 setEntry(i,j,_Mat.field().zero); 272 } 273 return *this; 274 } 275 276 template < class _Matrix > 277 void BlasSubmatrix<_Matrix>::random() 278 { 279 typename Field::RandIter r(_Mat.field()); 280 Element temp; _Mat.field().init(temp); 281 for (size_t i = 0 ; i < rowdim() ; ++i) 282 for (size_t j = 0 ; j < coldim() ; ++j) { 283 setEntry(i,j,r.random(temp)); 284 } 285 return ; 286 } 287 288} // LinBox 289 290 291////////////////// 292// DIMENSIONS // 293////////////////// 294 295namespace LinBox 296{ 297 298 template < class _Matrix > 299 size_t BlasSubmatrix<_Matrix>::rowdim() const 300 { 301 return _row ; 302 } 303 304 template < class _Matrix > 305 size_t BlasSubmatrix<_Matrix>::coldim() const 306 { 307 return _col ; 308 } 309 310 template < class _Matrix > 311 size_t BlasSubmatrix<_Matrix>::getStride() const 312 { 313 return _stride; 314 } 315 316 317} // LinBox 318 319////////////////// 320// ELEMENTS // 321////////////////// 322 323namespace LinBox 324{ 325 326 327 template < class _Matrix > 328 typename BlasSubmatrix<_Matrix>::pointer 329 BlasSubmatrix<_Matrix>::getPointer() 330 { 331 return _Mat.getPointer()+_off; 332 } 333 334 template < class _Matrix > 335 typename BlasSubmatrix<_Matrix>::const_pointer 336 BlasSubmatrix<_Matrix>::getPointer() const 337 { 338 return _Mat.getPointer()+_off; 339 } 340 341 template < class _Matrix > 342 typename BlasSubmatrix<_Matrix>::pointer 343 BlasSubmatrix<_Matrix>::getWritePointer() 344 { 345 return (_Mat.getWritePointer()+_off); 346 } 347 348 template < class _Matrix > 349 const typename LinBox::BlasSubmatrix<_Matrix>::Element & BlasSubmatrix<_Matrix>::setEntry (size_t i, size_t j, const Element &a_ij) 350 { 351 return _Mat.setEntry (_r0 + i, _c0 + j, a_ij); 352 } 353 354 template < class _Matrix > 355 typename LinBox::BlasSubmatrix<_Matrix>::Element &BlasSubmatrix<_Matrix>::refEntry (size_t i, size_t j) 356 { 357 return _Mat.refEntry ( _r0+i, _c0+j ); 358 } 359 360 template < class _Matrix > 361 const typename LinBox::BlasSubmatrix<_Matrix>::Element &BlasSubmatrix<_Matrix>::getEntry (size_t i, size_t j) const 362 { 363 return _Mat.getEntry ( _r0+i , _c0+j ); 364 } 365 366 template < class _Matrix > 367typename LinBox::BlasSubmatrix<_Matrix>::Element &BlasSubmatrix<_Matrix>::getEntry (Element &x, size_t i, size_t j) const 368 { 369 return _Mat.getEntry (x, _r0+i , _c0+j ); 370 } 371 372} 373 374/////////////////// 375// TRANSPOSE &AL // 376/////////////////// 377 378namespace LinBox 379{ 380} 381 382/////////////////// 383// ITERATORS // 384/////////////////// 385 386namespace LinBox 387{ 388 389 390 /*! Raw Iterators. 391 * @ingroup iterators 392 * 393 * The raw iterator is a method for accessing all entries in the matrix 394 * in some unspecified order. This can be used, e.g. to reduce all 395 * matrix entries modulo a prime before passing the matrix into an 396 * algorithm. 397 */ 398 template < class _Matrix > 399 class BlasSubmatrix<_Matrix>::Iterator { 400 public: 401 Iterator (){} 402 403 /*! @internal 404 * @brief NO DOC 405 */ 406 Iterator (const typename BlasMatrix<typename _Matrix::Field, typename _Matrix::Rep>::Iterator& cur, 407 const size_t c_dim, 408 const size_t stride, 409 const size_t c_idx) : 410 _cur (cur), _c_dim (c_dim), _stride(stride), _c_idx (c_idx) 411 {} 412 413 /*! @internal 414 * @brief copy operator. 415 * @param r Iterator to copy. 416 */ 417 Iterator& operator = (const Iterator& r) 418 { 419 _cur = r._cur; 420 _c_dim = r._c_dim; 421 _stride = r._stride; 422 _c_idx = r._c_idx; 423 return *this; 424 } 425 426 /*! @internal 427 * increment. 428 * ?? 429 */ 430 Iterator& operator ++() 431 { 432 if (_c_idx < _c_dim - 1){ 433 ++_cur; ++_c_idx; 434 } 435 else { 436 _cur = _cur + _stride - _c_dim + 1; 437 _c_idx = 0; 438 } 439 440 return *this; 441 } 442 443 /*! @internal 444 * increment. 445 * ?? 446 */ 447 Iterator& operator++ (int) 448 { 449 return this->operator++ (); 450 } 451 452 453 /*! @internal 454 * @brief operator !=. 455 * @param r Iterator to test inequaltity from. 456 */ 457 bool operator != (const Iterator& r) const 458 { 459 return (_cur != r._cur || _c_dim != r._c_dim) || (_stride != r._stride) || (_c_idx != r._c_idx); 460 } 461 462 //! @internal operator *. 463 typename _Matrix::Element& operator * () 464 { 465 return *_cur; 466 } 467 468 //! @internal operator *. 469 const typename _Matrix::Element& operator * () const 470 { 471 return *_cur; 472 } 473 474 protected: 475 typename BlasMatrix<typename _Matrix::Field,typename _Matrix::Rep>::Iterator _cur; 476 size_t _c_dim; 477 size_t _stride; 478 size_t _c_idx; 479 }; 480 481 /*! Raw Iterators (const version). 482 * @ingroup iterators 483 * The raw iterator is a method for accessing all entries in the matrix 484 * in some unspecified order. This can be used, e.g. to reduce all 485 * matrix entries modulo a prime before passing the matrix into an 486 * algorithm. 487 */ 488 template < class _Matrix > 489 class BlasSubmatrix<_Matrix>::ConstIterator { 490 public: 491 //! @internal Null constructor 492 ConstIterator (){} 493 494 495 /*! @internal 496 * @brief NO DOC 497 */ 498 ConstIterator (const typename BlasMatrix<typename _Matrix::Field, typename _Matrix::Rep>::ConstIterator& cur, 499 const size_t c_dim, 500 const size_t stride, 501 const size_t c_idx) : 502 _cur (cur), _c_dim (c_dim), _stride(stride), _c_idx (c_idx) 503 {} 504 505 /*! @internal 506 * @brief copy operator. 507 * @param r Iterator to copy. 508 */ 509 ConstIterator& operator = (const ConstIterator& r) 510 { 511 _cur = r._cur; 512 _c_dim = r._c_dim; 513 _stride = r._stride; 514 _c_idx = r._c_idx; 515 return *this; 516 } 517 518 /*! @internal 519 * increment. 520 * ?? 521 */ 522 ConstIterator& operator ++() 523 { 524 if (_c_idx < _c_dim - 1){ 525 ++_cur; ++_c_idx; 526 } 527 else { 528 linbox_check(_stride > _c_dim); 529 _cur = _cur + (ptrdiff_t)(_stride - _c_dim + 1); 530 _c_idx = 0; 531 } 532 533 return *this; 534 } 535 536 /*! @internal 537 * increment. 538 * ?? 539 */ 540 ConstIterator& operator++ (int) 541 { 542 return this->operator++ (); 543 } 544 545 /*! @internal 546 * @brief operator !=. 547 * @param r Iterator to test inequaltity from. 548 */ 549 bool operator != (const ConstIterator& r) const 550 { 551 return (_cur != r._cur) || (_c_dim != r._c_dim) || (_stride != r._stride) || (_c_idx != r._c_idx); 552 } 553 554 //! @internal operator *. 555 const typename BlasSubmatrix<_Matrix>::Element& operator * () const 556 { 557 return *_cur; 558 } 559 560 protected: 561 typename BlasMatrix<typename _Matrix::Field, typename _Matrix::Rep>::ConstIterator _cur; 562 size_t _c_dim; 563 size_t _stride; 564 size_t _c_idx; 565 }; 566 567#if 0 568 template < class _Matrix > 569 class BlasSubmatrix<_Matrix>::ConstIterator { 570 public: 571 ConstIterator (){} 572 573 ConstIterator ( const typename BlasMatrix< _Field>::ConstIterator& cur, 574 size_t cont_len, 575 size_t gap_len) : 576 _beg (beg), _cur (cur), _cont_len (cont_len), _gap_len (gap_len) 577 {} 578 579 ConstIterator& operator = (const Iterator& r) 580 { 581 _cur = r._cur; 582 _beg = r._beg; 583 _cont_len = r._cont_len; 584 _gap_len = r._gap_len; 585 return *this; 586 } 587 588 ConstIterator& operator = (const ConstIterator& r) 589 { 590 _cur = r._cur; 591 _beg = r._beg; 592 _cont_len = r._cont_len; 593 _gap_len = r._gap_len; 594 return *this; 595 } 596 597 ConstIterator& operator++() 598 { 599 if (((_cur - _beg + 1) % _cont_len) != 0) 600 ++_cur; 601 else { 602 _cur = _cur + _gap_len + 1; 603 _beg = _beg + _gap_len + _cont_len; 604 } 605 return *this; 606 } 607 608 ConstIterator operator++(int) 609 { 610 ConstIterator tmp = *this; 611 this->operator++(); 612 return tmp; 613 } 614 615 bool operator != (const ConstIterator& r) const 616 { 617 return (_cur != r._cur) || (_beg != r._beg) || (_cont_len != r._cont_len) || (_gap_len != r._gap_len); 618 } 619 620 const typename _Field::Element& operator*() 621 { return *_cur; } 622 623 _Field& operator*() 624 { return *_cur; } 625 626 const _Field& operator*() const 627 { return *_cur; } 628 629 protected: 630 typename BlasMatrix< _Field>::ConstIterator _beg; 631 typename BlasMatrix< _Field>::ConstIterator _cur; 632 size_t _cont_len; 633 size_t _gap_len; 634 }; 635#endif 636 637 template < class _Matrix > 638 typename BlasSubmatrix<_Matrix>::Iterator BlasSubmatrix<_Matrix>::Begin () 639 { 640 return Iterator (_Mat.Begin () + (ptrdiff_t)( _off ), 641 _col, _stride, 0); 642 } 643 644 template < class _Matrix > 645 typename BlasSubmatrix<_Matrix>::Iterator BlasSubmatrix<_Matrix>::End () 646 { 647 return Iterator (_Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + _off ), 648 _col, _stride, 0); 649 } 650 651 template < class _Matrix > 652 typename BlasSubmatrix<_Matrix>::ConstIterator BlasSubmatrix<_Matrix>::Begin () const 653 { 654 return ConstIterator (_Mat.Begin () +(ptrdiff_t) ( _off ), 655 _col, _stride, 0); 656 } 657 658 template < class _Matrix > 659 typename BlasSubmatrix<_Matrix>::ConstIterator BlasSubmatrix<_Matrix>::End () const 660 { 661 return ConstIterator (_Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + _off ), 662 _col, _stride, 0); 663 } 664 665#if 0 666 template < class _Matrix > 667 typename BlasSubmatrix<_Matrix>::ConstIterator BlasSubmatrix<_Matrix>::Begin () const 668 { 669 return ConstIterator (_Mat.Begin () +(ptrdiff_t) ( _off ), 670 _Mat.Begin () +(ptrdiff_t) ( _off ), 671 _col, _stride - _col); 672 } 673 674 template < class _Matrix > 675 typename BlasSubmatrix<_Matrix>::ConstIterator BlasSubmatrix<_Matrix>::End () const 676 { 677 return ConstIterator (_Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + _off ), 678 _Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + _off ), 679 _col, _stride - _col); 680 } 681#endif 682 683 /*! Raw Indexed Iterator. 684 * @ingroup iterators 685 * 686 * Like the raw iterator, the indexed iterator is a method for 687 * accessing all entries in the matrix in some unspecified order. 688 * At each position of the the indexed iterator, it also provides 689 * the row and column indices of the currently referenced entry. 690 * This is provided through it's \c rowIndex() and \c colIndex() functions. 691 */ 692 template < class _Matrix > 693 class BlasSubmatrix<_Matrix>::IndexedIterator { 694 public: 695 IndexedIterator (){} 696 697 IndexedIterator (const typename BlasMatrix<typename _Matrix::Field,typename _Matrix::Rep>::Iterator& cur, 698 size_t c_dim, 699 size_t stride, 700 size_t r_idx, 701 size_t c_idx) : 702 _cur (cur), _c_dim (c_dim), _stride (stride), _r_idx(r_idx), _c_idx (c_idx) 703 {} 704 705 IndexedIterator& operator = (const IndexedIterator& r) 706 { 707 _cur = r._cur; 708 _stride = r._stride; 709 _c_dim = r._c_dim; 710 _r_idx = r._r_idx; 711 _c_idx = r._c_idx; 712 return *this; 713 } 714 715 IndexedIterator& operator++() 716 { 717 if (_c_idx < _c_dim - 1){ 718 ++_c_idx; 719 ++_cur; 720 } 721 else { 722 _cur = _cur + _stride - _c_dim + 1; 723 _c_idx = 0; 724 ++_r_idx; 725 } 726 return *this; 727 } 728 729 IndexedIterator& operator--() 730 { 731 if (_c_idx > 0){ 732 --_c_idx; 733 --_cur; 734 } 735 else { 736 _cur = _cur - _stride + _c_dim -1; 737 _c_idx = 0; 738 --_r_idx; 739 } 740 return *this; 741 } 742 743 IndexedIterator operator++(int) 744 { 745 IndexedIterator tmp = *this; 746 this->operator++(); 747 return tmp; 748 } 749 750 IndexedIterator operator--(int) 751 { 752 IndexedIterator tmp = *this; 753 this->operator--(); 754 return tmp; 755 } 756 757 bool operator != (const IndexedIterator& r) const 758 { 759 return ((_c_idx != r._c_idx) || (_r_idx != r._r_idx) ||(_stride != r._stride) || (_c_dim != r._c_dim) ); 760 } 761 762 const typename _Matrix::Field& operator*() const {return *_cur;} 763 764 typename _Matrix::Element& operator*() {return *_cur;} 765 766 size_t rowIndex () const { return _r_idx; } 767 768 size_t colIndex () const { return _c_idx; } 769 770 const typename _Matrix::Element& value () const {return *_cur;} 771 772 protected: 773 typename BlasMatrix<typename _Matrix::Field, typename _Matrix::Rep>::Iterator _cur; 774 size_t _stride; 775 size_t _c_dim; 776 size_t _r_idx; 777 size_t _c_idx; 778 }; 779 780 template < class _Matrix > 781 typename BlasSubmatrix<_Matrix>::IndexedIterator BlasSubmatrix<_Matrix>::IndexedBegin () 782 { 783 return IndexedIterator (_Mat.Begin () +(ptrdiff_t) ( (_off) ), 784 _col , _stride, 0, 0); 785 } 786 787 template < class _Matrix > 788 typename BlasSubmatrix<_Matrix>::IndexedIterator BlasSubmatrix<_Matrix>::IndexedEnd () 789 { 790 return IndexedIterator (_Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + (_col+_off) ), 791 _col, _stride, _row-1, _col-1); 792 } 793 794 /*! Raw Indexed Iterator (const version). 795 * @ingroup iterators 796 * 797 * Like the raw iterator, the indexed iterator is a method for 798 * accessing all entries in the matrix in some unspecified order. 799 * At each position of the the indexed iterator, it also provides 800 * the row and column indices of the currently referenced entry. 801 * This is provided through it's \c rowIndex() and \c colIndex() functions. 802 */ 803 template < class _Matrix > 804 class BlasSubmatrix<_Matrix>::ConstIndexedIterator { 805 public: 806 ConstIndexedIterator (){} 807 808 ConstIndexedIterator (const typename BlasMatrix<typename _Matrix::Field, typename _Matrix::Rep>::ConstIterator& cur, 809 size_t c_dim, 810 size_t stride, 811 size_t r_idx, 812 size_t c_idx) : 813 _cur (cur), _stride (stride), _c_dim (c_dim), _r_idx(r_idx), _c_idx (c_idx) 814 {} 815 816 ConstIndexedIterator& operator = (const IndexedIterator& r) 817 { 818 _cur = r._cur; 819 _stride = r._stride; 820 _c_dim = r._c_dim; 821 _r_idx = r._r_idx; 822 _c_idx = r._c_idx; 823 return *this; 824 } 825 826 ConstIndexedIterator& operator = (const ConstIndexedIterator& r) 827 { 828 _cur = r._cur; 829 _stride = r._stride; 830 _c_dim = r._c_dim; 831 _r_idx = r._r_idx; 832 _c_idx = r._c_idx; 833 return *this; 834 } 835 836 ConstIndexedIterator& operator++() 837 { 838 if (_c_idx < _c_dim - 1){ 839 ++_c_idx; 840 ++_cur; 841 } 842 else { 843 _cur = _cur + _stride - _c_dim +1; 844 _c_idx = 0; 845 ++_r_idx; 846 } 847 return *this; 848 } 849 850 IndexedIterator& operator--() 851 { 852 if (_c_idx > 0){ 853 --_c_idx; 854 --_cur; 855 } 856 else { 857 _cur = _cur - _stride + _c_dim -1; 858 _c_idx = 0; 859 --_r_idx; 860 } 861 return *this; 862 } 863 864 ConstIndexedIterator operator++(int) 865 { 866 ConstIndexedIterator tmp = *this; 867 this->operator++(); 868 return tmp; 869 } 870 871 ConstIndexedIterator operator--(int) 872 { 873 ConstIndexedIterator tmp = *this; 874 this->operator--(); 875 return tmp; 876 } 877 878 size_t rowIndex () const 879 { 880 return _r_idx; 881 } 882 883 size_t colIndex () const 884 { 885 return _c_idx; 886 } 887 888 bool operator != (const ConstIndexedIterator& r) const 889 { 890 return ((_c_idx != r._c_idx) || (_r_idx != r._r_idx) ||(_stride != r._stride) || (_c_dim != r._c_dim) ); 891 } 892 893 const typename _Matrix::Element& operator*() const 894 { 895 return *_cur; 896 } 897 898 899 friend std::ostream& operator<<(std::ostream& out, const ConstIndexedIterator m) 900 { 901 return out /* << m._cur << ' ' */ 902 << m._stride << ' ' 903 << m._c_dim << ' ' 904 << m._r_idx << ' ' 905 << m._c_idx; 906 } 907 908 const typename _Matrix::Element & value() const 909 { 910 return this->operator*(); 911 912 } 913 914 protected: 915 typename BlasMatrix<typename _Matrix::Field, typename _Matrix::Rep>::ConstIterator _cur; 916 size_t _stride; 917 size_t _c_dim; 918 size_t _r_idx; 919 size_t _c_idx; 920 }; 921 922 template < class _Matrix > 923 typename BlasSubmatrix<_Matrix>::ConstIndexedIterator BlasSubmatrix<_Matrix>::IndexedBegin () const 924 { 925 return ConstIndexedIterator (_Mat.Begin () +(ptrdiff_t) ( _off ), 926 _row, _stride, 0, 0); 927 } 928 929 template < class _Matrix > 930 typename BlasSubmatrix<_Matrix>::ConstIndexedIterator BlasSubmatrix<_Matrix>::IndexedEnd () const 931 { 932 return ConstIndexedIterator (_Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + (_off+_col) ), 933 _col, _stride, _row-1, _col-1); 934 } 935 936 //////// 937 template < class _Matrix > 938 typename BlasSubmatrix<_Matrix>::RowIterator BlasSubmatrix<_Matrix>::rowBegin () 939 { 940 return RowIterator (_Mat.Begin () +(ptrdiff_t) ( _off ), 941 _col, _stride); 942 } 943 944 template < class _Matrix > 945 typename BlasSubmatrix<_Matrix>::RowIterator BlasSubmatrix<_Matrix>::rowEnd () 946 { 947 return RowIterator (_Mat.Begin () +(ptrdiff_t) ( (_row) * _stride + _off ), 948 _col, _stride); 949 } 950 951 template < class _Matrix > 952 typename BlasSubmatrix<_Matrix>::ConstRowIterator BlasSubmatrix<_Matrix>::rowBegin () const 953 { 954 return ConstRowIterator (_Mat.Begin () + (ptrdiff_t)( _off ), 955 _col, _stride); 956 } 957 958 template < class _Matrix > 959 typename BlasSubmatrix<_Matrix>::ConstRowIterator BlasSubmatrix<_Matrix>::rowEnd () const 960 { 961 return ConstRowIterator (_Mat.Begin () + (ptrdiff_t)( (_row) * _stride + _off ), 962 _col, _stride); 963 } 964 965 template < class _Matrix > 966 typename BlasSubmatrix<_Matrix>::ColIterator BlasSubmatrix<_Matrix>::colBegin () 967 { 968 return ColIterator (_Mat.Begin () + (ptrdiff_t)( _off ), 969 _stride, _row); 970 } 971 972 template < class _Matrix > 973 typename BlasSubmatrix<_Matrix>::ColIterator BlasSubmatrix<_Matrix>::colEnd () 974 { 975 return ColIterator (_Mat.Begin () + (ptrdiff_t)( (_col) + _off ), 976 _stride, _row); 977 } 978 979 template < class _Matrix > 980 typename BlasSubmatrix<_Matrix>::ConstColIterator BlasSubmatrix<_Matrix>::colBegin () const 981 { 982 return ConstColIterator (_Mat.Begin () + (ptrdiff_t)( _off ), 983 _stride, _row); 984 } 985 986 template < class _Matrix > 987 typename BlasSubmatrix<_Matrix>::ConstColIterator BlasSubmatrix<_Matrix>::colEnd () const 988 { 989 return ConstColIterator (_Mat.Begin () + (ptrdiff_t)( (_col) + _off ), 990 _stride, _row); 991 } 992 993 /* operators */ 994 template < class _Matrix > 995 typename BlasSubmatrix<_Matrix>::Row BlasSubmatrix<_Matrix>::operator[] (size_t i) 996 { 997 return Row (_Mat.Begin () +(ptrdiff_t) (_r0+i) * _stride, _Mat.Begin () + (ptrdiff_t)((_r0+i) * _stride + _stride) ); 998 } 999 1000 template < class _Matrix > 1001 typename BlasSubmatrix<_Matrix>::ConstRow BlasSubmatrix<_Matrix>::operator[] (size_t i) const 1002 { 1003 return Row (_Mat.Begin () + (ptrdiff_t)(_r0+i) * _stride, _Mat.Begin () + (ptrdiff_t)((_r0+i) * _stride + _stride) ); 1004 } 1005 1006} // LinBox 1007 1008/////////////////// 1009// I/O // 1010/////////////////// 1011 1012namespace LinBox 1013{ 1014 1015 template < class _Matrix > 1016 std::istream& BlasSubmatrix<_Matrix>::read (std::istream &file) 1017 { 1018#if 0 1019 Iterator p; 1020 int m,n; 1021 char c; 1022 file>>m>>n>>c; 1023 1024 if (m*n < _row*_col) 1025 cerr<<"NOT ENOUGH ELEMENT TO READ\n"; 1026 else { 1027 for (p = Begin (); p != End (); ++p) { 1028 integer tmp; 1029 file>>tmp;cout<<tmp<<endl; 1030 //file.ignore(1); 1031 _Mat.field().read (file, *p); 1032 } 1033 } 1034#endif 1035 1036 1037 Iterator p; 1038 int m=0,n=0; 1039 char c='\0'; 1040 file>>m>>n>>c; 1041 // std::cout << m << 'x' << n << ':' << c << std::endl; 1042 1043 // this is bogus!! -bds 1044 _row = m; _col = n; 1045 1046 // resize(_row,_col); 1047 1048 if ((c != 'M') && (c != 'm')) { 1049 for (p = Begin (); p != End (); ++p) { 1050 //file.ignore(1); 1051 _Mat.field().read (file, *p); 1052 } 1053 1054 } 1055 else { // sparse file format - needs fixing 1056 int i=0, j=0; 1057 while (true) 1058 { 1059 file >> i >> j; 1060 //file.ignore(1); 1061 //if (! file) break; 1062 if (i+j <= 0) break; 1063 // std::cout << i << ',' << j << ':' ; 1064 _Mat.field().read (file, _Mat.refEntry(i-1, j-1)); 1065 } 1066 } 1067 1068 return file; 1069 } 1070 1071 template <class _Matrix> 1072 std::ostream &BlasSubmatrix< _Matrix >::write (std::ostream &os, 1073 Tag::FileFormat f ) const 1074 { 1075 1076 ConstRowIterator p; 1077 switch(f) { 1078 case (Tag::FileFormat::MatrixMarket ) : /* Matrix Market */ 1079 { 1080 writeMMArray(os, *this, "BlasSubmatrix"); 1081 } 1082 break; 1083 case (Tag::FileFormat::Plain) : /* raw output */ 1084 { 1085 integer c; 1086 int wid; 1087 1088 _Mat.field().cardinality (c); 1089 1090 if (c >0) 1091 wid = (int) ceil (log ((double) c) / M_LN10); 1092 else { 1093// integer tmp; 1094// size_t max=0; 1095// ConstIterator it = Begin(); 1096// for (; it != End(); ++it){ 1097// _Mat.field().convert(tmp,*it); 1098// if (tmp.bitsize() > max) 1099// max= tmp.bitsize(); 1100// } 1101// wid= (int) ceil ((double)max / M_LN10)+1; 1102 wid=1000; 1103 1104 } 1105 1106 for (p = rowBegin (); p != rowEnd ();++p) { 1107 typename ConstRow::const_iterator pe; 1108 1109 os << " [ "; 1110 1111 for (pe = p->begin (); pe != p->end (); ++pe) { 1112 os.width (wid); 1113 /*! @warning 1114 * matrix base does not provide this field(), maybe should? 1115 * _Mat.field ().write (os, *pe); 1116 * os << *pe; 1117 * fixed by using extra field 1118 */ 1119 1120 _Mat.field().write (os, *pe); 1121 os << " "; 1122 } 1123 1124 os << "]" << std::endl; 1125 } 1126 } 1127 break; 1128 case (Tag::FileFormat::Maple) : /* maple format */ 1129 { 1130 1131 os << "Matrix( " << rowdim() << ',' << coldim() << ",\n[" ; 1132 for (p = rowBegin (); p != rowEnd (); ) { 1133 typename ConstRow::const_iterator pe; 1134 if (p!=rowBegin()) os << ' '; 1135 os << "[ "; 1136 1137 for (pe = p->begin (); pe != p->end (); ) { 1138 _Mat.field().write (os, *pe); 1139 ++pe ; 1140 if (pe != p->end()) 1141 os << ", "; 1142 } 1143 1144 os << "]" ; 1145 ++p ; 1146 if (p != rowEnd() ) 1147 os << ',' << std::endl;; 1148 1149 } 1150 os << "])" ; 1151 } 1152 break; 1153 case (Tag::FileFormat::HTML) : /* HTML format */ 1154 { 1155 1156 os << "<table border=\"1\">" ; 1157 for (p = rowBegin (); p != rowEnd (); ) { 1158 typename ConstRow::const_iterator pe; 1159 1160 os << "<tr>"; 1161 1162 for (pe = p->begin (); pe != p->end (); ) { 1163 _Mat.field().write (os<< "<td>", *pe)<<"</td>"; 1164 ++pe ; 1165 } 1166 1167 os << "</tr>" << std::endl; 1168 ++p ; 1169 1170 } 1171 os << "</table>" ; 1172 } 1173 break; 1174 case (Tag::FileFormat::LaTeX) : /* LaTex format */ 1175 { 1176 os << "\\begin{pmatrix} " << std::endl; 1177 for (p = rowBegin (); p != rowEnd (); ) { 1178 typename ConstRow::const_iterator pe; 1179 1180 for (pe = p->begin (); pe != p->end (); ) { 1181 _Mat.field().write (os, *pe); 1182 ++pe ; 1183 if (pe != p->end()) 1184 os << "& "; 1185 } 1186 1187 os << "\\\\" << std::endl; 1188 ++p ; 1189 1190 } 1191 os << "\\end{pmatrix}" ; 1192 } 1193 break; 1194 default : /* this is an error */ 1195 { 1196 throw LinBoxError("unknown format to write matrix in"); 1197 } 1198 } 1199 1200 return os; 1201 } 1202 1203 template <class _Matrix> 1204 template<typename _Tp1, class _Rep2> 1205 struct BlasSubmatrix< _Matrix>::rebind { 1206 typedef BlasMatrix<_Tp1,_Rep2> other; 1207 1208 void operator() (other & Ap, const Self_t& A) { 1209 typedef typename BlasSubmatrix<_Matrix>::ConstIterator ConstSelfIterator ; 1210 typedef typename other::Iterator OtherIterator ; 1211 OtherIterator Ap_i; 1212 ConstSelfIterator A_i; 1213 Hom<Field, _Tp1> hom(A. field(), Ap. field()); 1214 for (A_i = A. Begin(), Ap_i = Ap.Begin(); 1215 A_i != A. End(); ++ A_i, ++ Ap_i) 1216 hom.image (*Ap_i, *A_i); 1217 } 1218 }; 1219 1220} // LinBox 1221 1222////////////////// 1223// MISC // 1224////////////////// 1225 1226namespace LinBox 1227{ 1228} // LinBox 1229#endif // __LINBOX_densematrix_blas_submatrix_INL 1230 1231// Local Variables: 1232// mode: C++ 1233// tab-width: 4 1234// indent-tabs-mode: nil 1235// c-basic-offset: 4 1236// End: 1237// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 1238