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_matrix_INL 37#define __LINBOX_densematrix_blas_matrix_INL 38 39///////////////// 40// PRIVATE // 41///////////////// 42 43namespace LinBox 44{ 45 template<class _Field, class _Rep> 46 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const BlasMatrix< _Field, _Rep >& A) 47 { 48#ifndef NDEBUG 49 if (!areFieldEqual(A.field(),field())) { 50 A.field().write(std::cout) << "!=" ; 51 field().write(std::cout) << std::endl; 52 } 53#endif 54 linbox_check( areFieldEqual(A.field(),field())); 55 createBlasMatrix(A,MatrixContainerCategory::BlasContainer()); 56 } 57 58 template<class _Field, class _Rep> 59 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const Element * v) 60 { 61 // Element * iter_v = const_cast<Element*>(v) ; 62 Element * v_end = const_cast<Element*>(v+(_col*_row)) ; 63 // Iterator iter_addr = this->Begin(); 64 Element * iter_addr = _ptr ; 65 for (; v != v_end ; ++v, ++iter_addr) 66 { 67 field().init(*iter_addr); 68 field().assign(*iter_addr,*v); 69 } 70 } 71 72 template<class _Field, class _Rep> 73 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const std::vector<Element> & v) 74 { 75 typename std::vector< Element>::const_iterator iter_value = v.begin(); 76 Iterator iter_addr = this->Begin(); 77 for (;iter_value != v.end(); ++iter_value,++iter_addr) 78 { 79 field().init(*iter_addr); 80 field().assign(*iter_addr,*iter_value); 81 } 82 } 83 84 85 template<class _Field, class _Rep> 86 template <class _Matrix> 87 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const _Matrix& A, 88 MatrixContainerCategory::BlasContainer) 89 { 90 // std::cout << "creator 4" << std::endl; 91 linbox_check( areFieldEqual(A.field(),field())); 92#if 0 93 typename _Matrix::ConstIterator iter_value = A.Begin(); 94 Iterator iter_addr = this->Begin(); 95 for (;iter_value != A.End(); ++iter_value,++iter_addr) 96 { 97 field().init(*iter_addr); 98 field().assign(*iter_addr, *iter_value); 99 } 100#else 101 for (size_t i = 0 ; i < A.rowdim() ; ++i) 102 for (size_t j = 0 ; j < A.coldim() ; ++j) 103 _rep[i*_col+j] = A.getEntry(i,j) ; 104#endif 105 } 106 107 template<class _Field, class _Rep> 108 template <class Matrix> 109 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const Matrix& A, 110 MatrixContainerCategory::Container) 111 { 112 // std::cout << "creator 5" << std::endl; 113 linbox_check( areFieldEqual(A.field(),field())); 114 // const Field & F = A.field(); 115 //!@bug With both iterators, it is Segfaulting !!!! 116 typename Matrix::ConstIndexedIterator iter_index = A.IndexedBegin(); 117 for (;iter_index != A.IndexedEnd(); ++iter_index) 118 setEntry( iter_index.rowIndex(), 119 iter_index.colIndex(), 120 A.getEntry(iter_index.rowIndex(),iter_index.colIndex()) 121 ); 122 } 123 124 template<class _Field, class _Rep> 125 template <class Matrix> 126 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const Matrix& A, 127 MatrixContainerCategory::Blackbox) 128 { 129 // std::cout << "creator 6" << std::endl; 130 linbox_check( areFieldEqual(A.field(),field()) ); 131 132 BlasVector<Field> e(A.field(),A.coldim(), field().zero), tmp(A.field(),A.rowdim()); 133 ColIterator col_p; 134 135 typename BlasMatrix< _Field, _Rep >::Col::iterator elt_p; 136 typename BlasVector<Field>::iterator e_p, tmp_p; 137 138 139 for (col_p = colBegin(), e_p = e.begin(); 140 e_p != e.end(); ++ col_p, ++ e_p) 141 { 142 143 field().assign(*e_p, field().one); 144 145 A.apply (tmp, e); 146 147 for (tmp_p = tmp.begin(), elt_p = col_p -> begin(); 148 tmp_p != tmp.end(); ++ tmp_p, ++ elt_p) 149 150 field().assign(*elt_p, *tmp_p); 151 152 field().assign(*e_p, field().zero); 153 } 154 } 155 156 template<class _Field, class _Rep> 157 template <class _Matrix> 158 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const _Matrix& A, 159 const size_t i0,const size_t j0, 160 const size_t m, const size_t n, 161 MatrixContainerCategory::Container) 162 { 163 linbox_check( areFieldEqual(A.field(),field() ) ); 164 165 typename _Matrix::ConstIterator iter_value = A.Begin(); 166 typename _Matrix::ConstIndexedIterator iter_index = A.IndexedBegin(); 167 168 for (;iter_value != A.End(); ++iter_value,++iter_index){ 169 int i,j; 170 i=(int)iter_index.rowIndex()-(int)i0; 171 j=(int)iter_index.colIndex()-(int)j0; 172 if (( i >= 0) && (j >= 0) && (i< (int)m) && (j < (int)n)) 173 setEntry(i, j, *iter_value); 174 } 175 } 176 177 template<class _Field, class _Rep> 178 template <class Matrix> 179 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const Matrix& A, 180 const size_t i0,const size_t j0, 181 const size_t m, const size_t n, 182 MatrixContainerCategory::BlasContainer) 183 { 184 linbox_check( areFieldEqual(A.field(),field() ) ); 185 186 typename Matrix::ConstIterator iter_value = A.Begin(); 187 typename Matrix::ConstIndexedIterator iter_index = A.IndexedBegin(); 188 189 for (;iter_value != A.End(); ++iter_value,++iter_index){ 190 int i,j; 191 i=(int)iter_index.rowIndex()-(int)i0; 192 j=(int)iter_index.colIndex()-(int)j0; 193 if ( (i>=0) && (j>=0) && (i< (int)m) && (j < (int)n)) 194 setEntry((size_t)i, (size_t)j, *iter_value); 195 } 196 } 197 198 template<class _Field, class _Rep> 199 template <class Matrix> 200 void BlasMatrix< _Field, _Rep >::createBlasMatrix (const Matrix& A, 201 const size_t i0,const size_t j0, 202 const size_t m, const size_t n, 203 MatrixContainerCategory::Blackbox) 204 { 205 linbox_check( areFieldEqual(A.field(),field() ) ); 206 207 208 BlasVector<Field> e(A.field(),A.coldim(), field().zero), tmp(A.field(),A.rowdim()); 209 ColIterator col_p; 210 211 typename BlasMatrix< _Field, _Rep >::Col::iterator elt_p; 212 typename BlasVector<Element>::iterator e_p, tmp_p; 213 214 215 for (col_p = colBegin(), e_p = e.begin()+(ptrdiff_t)j0; 216 e_p != e.begin()+(ptrdiff_t)(j0+n); ++ col_p, ++ e_p) { 217 218 field().assign(*e_p, field().one); 219 220 A.apply (tmp, e); 221 222 for (tmp_p = tmp.begin()+(ptrdiff_t)i0, elt_p = col_p -> begin(); 223 elt_p != col_p.end(); ++ tmp_p, ++ elt_p) { 224 field().assign(*elt_p, *tmp_p); 225 } 226 227 field().assign(*e_p, field().zero); 228 } 229 } 230 231} // LinBox 232 233////////////////// 234// CONSTRUCTORS // 235////////////////// 236 237namespace LinBox 238{ 239 240 241 template < class _Field, class _Rep > 242 BlasMatrix< _Field, _Rep >::BlasMatrix (const _Field &F) : 243 _row(0),_col(0),_rep(0) 244 // ,_use_fflas(false) 245 ,_ptr(NULL) 246 ,_field(&F),_MD(F),_VD(F) 247 // ,_AD(F) 248 { 249 //std::cout << "cstor 1 called" << std::endl; 250 _use_fflas = Protected::checkBlasApply(field(),_col); 251 } 252 253#if 0 254 template < class _Field, class _Rep > 255 BlasMatrix< _Field, _Rep >::BlasMatrix () : 256 _row(0),_col(0),_rep(0),_ptr(NULL), 257 _field(Field()),_MD(_field ),_VD(_field ) 258 {} 259#endif 260 261 template < class _Field, class _Rep > 262 void BlasMatrix< _Field, _Rep >::init(const _Field &F, const size_t & r, const size_t & c) 263 { 264 _field = &F; _row = r; _col = c; 265 _rep.resize(r*c, F.zero); 266 _ptr = _rep.data(); 267 _VD.init(F); _MD.init(F); 268 // _AD.init(F); 269 } 270 271 template < class _Field, class _Rep > 272 BlasMatrix< _Field, _Rep >::BlasMatrix ( const _Field &F, const size_t & m, const size_t & n) : 273 _row(m),_col(n),_rep(_row*_col, F.zero),_ptr(_rep.data()), 274 _field(&F),_MD(F),_VD(F) 275 // ,_AD(F) 276 { 277 // std::cout << "cstor 2 called" << std::endl; 278 _use_fflas = Protected::checkBlasApply(field(),_col); 279 } 280 281 282 template < class _Field, class _Rep > 283 BlasMatrix< _Field, _Rep >::BlasMatrix(MatrixStream<_Field>& ms) : 284 _row(0),_col(0),_rep(0), 285 _field(&(ms.getField())),_MD(field() ),_VD(field() ) 286 // ,_AD(field()) 287 { 288 // std::cout << "cstor 3 called" << std::endl; 289 if( !ms.getArray(_rep) || !ms.getDimensions(_row, _col) ) 290 throw ms.reportError(__FUNCTION__,__LINE__); 291 _ptr = _rep.data(); 292 _use_fflas = Protected::checkBlasApply(field(), _col); 293 } 294 295 template < class _Field, class _Rep > 296 template <class StreamVector> 297 BlasMatrix< _Field, _Rep >::BlasMatrix (const Field &F, VectorStream<StreamVector> &stream) : 298 _row(stream.size ()), _col(stream.dim ()), _rep(_row*_col), _ptr(_rep.data()), 299 _field (&F), _MD (F), _VD(F) 300 // ,_AD(F) 301 { 302 // std::cout << "cstor 4 called" << std::endl; 303 StreamVector tmp(F); 304 typename BlasMatrix<Field,_Rep>::RowIterator p; 305 306 VectorWrapper::ensureDim (tmp, stream.dim ()); 307 308 for (p = BlasMatrix<Field,_Rep>::rowBegin (); p != BlasMatrix<Field,_Rep>::rowEnd (); ++p) { 309 stream >> tmp; 310 _VD.copy (*p, tmp); 311 } 312 _use_fflas = Protected::checkBlasApply(field(), _col); 313 } 314 315 template < class _Field, class _Rep > 316 template <class Matrix> 317 BlasMatrix< _Field, _Rep >::BlasMatrix (const Matrix &A) : 318 _row(A.rowdim()),_col(A.coldim()),_rep(_row*_col),_ptr(_rep.data()), 319 _field(&(A.field())),_MD(field() ),_VD(field() ) 320 // ,_AD(field()) 321 { 322 // std::cout << "cstor 5 called" << std::endl; 323 // makePointer(); 324 _use_fflas = Protected::checkBlasApply(field(), _col); 325 createBlasMatrix(A, typename MatrixContainerTrait<Matrix>::Type()); 326 } 327 328 template < class _Field, class _Rep > 329 template <class Matrix> 330 BlasMatrix< _Field, _Rep >::BlasMatrix (const Matrix& A, 331 const size_t &i0, const size_t &j0, 332 const size_t &m, const size_t &n) : 333 _row(m),_col(n),_rep(_row*_col),_ptr(_rep.data()), 334 _field(&(A.field())),_MD(field() ),_VD(field() ) 335 // ,_AD(field()) 336 { 337 // std::cout << "cstor 6 called" << std::endl; 338 _use_fflas = Protected::checkBlasApply(field(), _col); 339 // makePointer(); 340 createBlasMatrix(A, i0, j0, m, n, 341 typename MatrixContainerTrait<Matrix>::Type()); 342 } 343 344 template < class _Field, class _Rep > 345 template<class _Matrix> 346 BlasMatrix< _Field, _Rep >::BlasMatrix (const _Matrix &A, const _Field &F) : 347 _row(A.rowdim()), _col(A.coldim()),_rep(_row*_col),_ptr(_rep.data()), 348 _field(&F),_MD(field() ),_VD(field() ) 349 // ,_AD(field()) 350 { 351 // std::cout << "cstor 7 called" << std::endl; 352 _use_fflas = Protected::checkBlasApply(field(), _col); 353 typename _Matrix::template rebind<_Field>()(*this,A); 354 } 355 356 template < class _Field, class _Rep > 357 BlasMatrix< _Field, _Rep >::BlasMatrix (const BlasMatrix< _Field, _Rep >& A) : 358 _row(A.rowdim()), _col(A.coldim()),_rep(_row*_col),_ptr(_rep.data()), 359 _field(&(A.field())),_MD(field() ),_VD(field() ) 360 // ,_AD(field()) 361 { 362 // std::cout << "cstor 8 called" << std::endl; 363 _use_fflas = Protected::checkBlasApply(field(), _col); 364 // makePointer(); 365 createBlasMatrix(A); 366 } 367 368 template < class _Field, class _Rep > 369 BlasMatrix< _Field, _Rep >::BlasMatrix (const _Field &F, 370 const std::vector<typename _Field::Element>& v, 371 const size_t & m, const size_t & n) : 372 _row(m), _col(n),_rep(_row*_col),_ptr(_rep.data()), 373 _field(&F),_MD(field() ),_VD(field() ) 374 // ,_AD(field()) 375 { 376 // std::cout << "cstor 9 called" << std::endl; 377 linbox_check(v.size() == m*n); 378 // makePointer(); 379 _use_fflas = Protected::checkBlasApply(field(), _col); 380 createBlasMatrix(v); 381 } 382 383 template < class _Field, class _Rep > 384 BlasMatrix< _Field, _Rep >::BlasMatrix (const _Field &F, 385 const typename _Field::Element * v, 386 const size_t & m, const size_t & n) : 387 _row(m), _col(n),_rep(_row*_col),_ptr(_rep.data()), 388 _field(&F), _MD(field() ),_VD(field() ) 389 // ,_AD(field()) 390 { 391 // std::cout << "cstor 10 called" << std::endl; 392 // makePointer(); 393 _use_fflas = Protected::checkBlasApply(field(), _col); 394 createBlasMatrix(v); 395 } 396 397 template < class _Field, class _Rep > 398 BlasMatrix< _Field, _Rep >::~BlasMatrix () 399 { 400 // if (_ptr) 401 // free(_ptr); 402 } 403 404 template < class _Field, class _Rep > 405 void BlasMatrix< _Field, _Rep >::zero() 406 { 407 typename _Field::Element x; field().init(x); 408 for (size_t i = 0; i < rowdim(); ++i) 409 for (size_t j = 0; j < coldim(); ++j) 410 setEntry(i, j, field().zero); 411 } 412 413} // LinBox 414 415/////////////////// 416// I/O // 417/////////////////// 418 419namespace LinBox 420{ 421 422 template < class _Field, class _Rep > 423 std::istream &BlasMatrix< _Field, _Rep >::read (std::istream &file) 424 { 425 MatrixStream<Field> ms(field(), file); 426 if( !ms.getArray(_rep) || !ms.getDimensions(_row, _col) ) 427 throw ms.reportError(__FUNCTION__,__LINE__); 428 _ptr = _rep.data(); 429 _use_fflas = Protected::checkBlasApply(field(), _col); 430 return file; 431 } 432 433 template < class _Field, class _Rep > 434 BlasMatrix< _Field, _Rep >& BlasMatrix< _Field, _Rep >::operator= (const BlasMatrix< _Field, _Rep >& A) 435 { 436 if ( &A == this) 437 return *this; 438 439 _col = A.coldim(); 440 _row = A.rowdim(); 441 _rep = Rep(_row*_col); 442 _ptr = _rep.data() ; 443 // const_cast<_Field&>(_field ) = A.field(); 444 // changeField( A.field() ); 445 createBlasMatrix(A); 446 447 return *this; 448 } 449 450#if 0 /* loop rebind */ 451 template < class _Field, class _Rep > 452 template<typename _Tp1> 453 struct BlasMatrix< _Field, _Rep >::rebind { 454 typedef BlasMatrix<_Tp1> other; 455 456 void operator() (other & Ap, const Self_t& A, const _Tp1& F) { 457 // typedef typename BlasMatrix< _Field, _Rep >::ConstIndexedIterator ConstIndexedIterator ; 458 // typedef typename BlasMatrix< _Field, _Rep >::ConstIterator ConstIterator ; 459 ConstIterator iter_value = A.Begin(); 460 ConstIndexedIterator iter_index = A.IndexedBegin(); 461 typename _Tp1::Element tmp; 462 for (;iter_value != A.End(); ++iter_value,++iter_index){ 463 Ap.field().init(tmp); 464 Ap.field().assign(tmp, *iter_value); 465 Ap.setEntry(iter_index.rowIndex(), iter_index.colIndex(),tmp); 466 } 467 } 468 }; 469#endif 470 471#if 1 /* HOM */ 472 //! @bug other rep 473 template < class _Field, class _Rep > 474 template<typename _Tp1> 475 struct BlasMatrix< _Field, _Rep >::rebind { 476 typedef BlasMatrix<_Tp1,typename Vector<_Tp1>::Dense> other; 477 478 void operator() (other & Ap, const Self_t& A) { 479 // typedef Self_t::ConstIterator ConstSelfIterator ; 480 typedef typename BlasMatrix< _Field, _Rep >::ConstIterator ConstSelfIterator ; 481 typedef typename other::Iterator OtherIterator ; 482 OtherIterator Ap_i = Ap.Begin(); 483 ConstSelfIterator A_i = A.Begin(); 484 Hom<Field, _Tp1> hom(A. field(), Ap. field()) ; 485 for ( ; A_i != A. End(); ++ A_i, ++ Ap_i) 486 hom.image (*Ap_i, *A_i); 487 } 488 }; 489#endif 490 491 492} // LinBox 493 494////////////////// 495// DIMENSIONS // 496////////////////// 497 498namespace LinBox 499{ 500 501 template < class _Field, class _Rep > 502 size_t BlasMatrix< _Field, _Rep >::rowdim() const 503 { 504 return _row ; 505 } 506 507 template < class _Field, class _Rep > 508 size_t BlasMatrix< _Field, _Rep >::coldim() const 509 { 510 return _col ; 511 } 512 513 template < class _Field, class _Rep > 514 size_t BlasMatrix< _Field, _Rep >::getStride() const 515 { 516 return _col; 517 } 518 519 template < class _Field, class _Rep > 520 size_t& BlasMatrix< _Field, _Rep >::getWriteStride() 521 { 522 return _col; 523 } 524 525 template < class _Field, class _Rep > 526 void BlasMatrix< _Field, _Rep >::resize (const size_t & m, const size_t & n, const Element& val ) 527 { 528#ifndef NDEBUG 529 if (_col > 0 && _col != n) 530 std::cerr << " ***Warning*** you are resizing a matrix, possibly loosing data. " << std::endl; 531#endif 532 _row = m; 533 _col = n; 534 _rep.resize (m * n, val); 535 _ptr = _rep.data(); 536#if 0 537 if (_ptr) { 538 if (m && n) 539 realloc(_ptr,m*n*sizeof(Element)); 540 else { 541 free(_ptr); 542 _ptr=NULL ; 543 } 544 } 545 else 546 makePointer(); 547#endif 548 } 549 550} // LinBox 551 552////////////////// 553// ELEMENTS // 554////////////////// 555 556namespace LinBox 557{ 558 559 560 template < class _Field, class _Rep > 561 typename BlasMatrix< _Field, _Rep >::pointer 562 BlasMatrix< _Field, _Rep >::getPointer() 563 { 564 return _ptr; 565 } 566 567 template < class _Field, class _Rep > 568 typename BlasMatrix< _Field, _Rep >::const_pointer 569 BlasMatrix< _Field, _Rep >::getPointer() const 570 { 571 return (const_pointer)_ptr; 572 } 573 574 template < class _Field, class _Rep > 575 typename BlasMatrix< _Field, _Rep >::const_pointer 576 BlasMatrix< _Field, _Rep >::getConstPointer() const 577 { 578 return (const_pointer)_ptr; 579 } 580 581 582 template < class _Field, class _Rep > 583 typename BlasMatrix< _Field, _Rep >::pointer& 584 BlasMatrix< _Field, _Rep >::getWritePointer() 585 { 586 return _ptr; 587 } 588 589 template < class _Field, class _Rep > 590 const typename _Field::Element & BlasMatrix< _Field, _Rep >::setEntry (size_t i, size_t j, const Element &a_ij) 591 { 592// return _ptr[i * _col + j] = a_ij; 593 _ptr[i * _col + j] = a_ij; 594 return a_ij; 595 } 596 597 template < class _Field, class _Rep > 598 typename _Field::Element & BlasMatrix< _Field, _Rep >::refEntry (size_t i, size_t j) 599 { 600 return _ptr[i * _col + j]; 601 } 602 603 template < class _Field, class _Rep > 604 const typename _Field::Element & BlasMatrix< _Field, _Rep >::getEntry (size_t i, size_t j) const 605 { 606 return _rep[i * _col + j]; 607 } 608 609 template < class _Field, class _Rep > 610 typename _Field::Element & BlasMatrix< _Field, _Rep >::getEntry (Element &x, size_t i, size_t j) const 611 { 612 return x = _rep[i * _col + j]; 613 } 614 615} // LinBox 616 617/////////////////// 618// TRANSPOSE &AL // 619/////////////////// 620 621namespace LinBox 622{ 623 template < class _Field, class _Rep > 624 BlasMatrix< _Field, _Rep > BlasMatrix< _Field, _Rep >::transpose(BlasMatrix< _Field, _Rep > & tM) const 625 { 626 size_t r = this->rowdim() ; 627 size_t c = this->coldim() ; 628 linbox_check(tM.coldim() == r ); 629 linbox_check(tM.rowdim() == c); 630 for (size_t i = 0 ; i < r ; ++i) 631 for (size_t j = 0 ; j < c ; ++j) 632 tM.setEntry(j,i,this->getEntry(i,j)); 633 return tM; 634 } 635 636 namespace Protected 637 { 638 /*! @internal 639 * @brief In-Place Tranpose. 640 * Adapted from the Wikipedia article. 641 * @todo make it for strides and Submatrices 642 * @todo use specialized versions when available (eg dgetmi) 643 * @todo make transpose have an inplace=true default parameter 644 * (execpt maybe when below a threshold). 645 * @param m pointer to the beginning of a row-major matrix vector 646 * @param w rows in the matrix 647 * @param h cols in the matrix 648 */ 649 template<class T> 650 void transposeIP(T *m, size_t h, size_t w) 651 { 652 size_t start; 653 T tmp; 654 655 for (start = 0; start <= w * h - 1; ++start) { 656 size_t next, i ; 657 next = start; 658 i = 0; 659 do { 660 ++i; 661 next = (next % h) * w + next / h; 662 } while (next > start); 663 if (next < start || i == 1) 664 continue; 665 666 tmp = m[next = start]; 667 do { 668 i = (next % h) * w + next / h; 669 m[next] = (i == start) ? tmp : m[i]; 670 next = i; 671 } while (next > start); 672 } 673 } 674 } // Protected 675 676 template < class _Field, class _Rep > 677 template<bool _IP> 678 void BlasMatrix< _Field, _Rep >::transpose() 679 { 680 size_t r = this->rowdim() ; 681 size_t c = this->coldim() ; 682 683 if ( r == c) { 684 for (size_t i = 0 ; i < r ; ++i) 685 for (size_t j = i+1 ; j < c ; ++j) 686 std::swap(this->refEntry(i,j),this->refEntry(j,i)); 687 } 688 else { 689 if (!_IP) { 690 BlasMatrix< _Field, _Rep > MM(*this); 691 std::swap(_row,_col); 692 // iterating row first seems faster. 693#ifdef _BLOCKIT 694 size_t B ; 695 B = 1024; 696 697 for (size_t bi = 0 ; bi < r/B ; ++bi) { 698 for (size_t bj = 0 ; bj < c/B ; ++bj){ 699 for (size_t i = 0 ; i < B ; ++i) 700 for (size_t j = 0 ; j < B ; ++j) 701 this->setEntry(bj*B+j,bi*B+i, 702 MM.getEntry(bi*B+i,bj*B+j)); 703 } 704 } 705 for (size_t i = r/B*B ; i < r ; ++i) 706 for (size_t j = c/B*B ; j < c ; ++j) 707 this->setEntry(j,i, 708 MM.getEntry(i,j)); 709#else 710 for (size_t i = 0 ; i < r ; ++i) 711 for (size_t j = 0 ; j < c ; ++j) 712 this->setEntry(j,i, 713 MM.getEntry(i,j)); 714#endif 715 716 } 717 else { 718 Protected::transposeIP<Element>(_ptr,_row,_col); 719 std::swap(_row,_col); 720 } 721 } 722 } 723 724 template < class _Field, class _Rep > 725 void BlasMatrix< _Field, _Rep >::transpose() 726 { 727 this->transpose<false>(); 728 } 729 730 template < class _Field, class _Rep > 731 void BlasMatrix< _Field, _Rep >::reverseRows() 732 { 733 size_t r = this->rowdim()/2 ; 734 for (size_t i = 0 ; i < r ; ++i) { 735 _VD.swap( this->rowBegin()+i, 736 this->rowBegin()+(r-1-i) ); 737 } 738 739 } 740 741 template < class _Field, class _Rep > 742 void BlasMatrix< _Field, _Rep >::reverseCols() 743 { 744 size_t r = this->rowdim() ; 745 size_t c = this->coldim()/2 ; 746 for (size_t j = 0 ; j < c ; ++j) { 747 for (size_t i = 0 ; i < r ; ++i) { 748 std::swap(this->refEntry(i,j), 749 this->refEntry(i,c-j-1)); 750 751 } 752 } 753 } 754 755 template < class _Field, class _Rep > 756 void BlasMatrix< _Field, _Rep >::reverse() 757 { 758 size_t r = this->rowdim() ; 759 size_t c = this->coldim() ; 760 for (size_t j = 0 ; j < c ; ++j) { 761 for (size_t i = 0 ; i < r ; ++i) { 762 std::swap(this->refEntry(i,j), 763 this->refEntry(r-i-1,c-j-1)); 764 765 } 766 } 767 } 768} // LinBox 769 770/////////////////// 771// ITERATORS // 772/////////////////// 773 774namespace LinBox 775{ 776 777 template < class _Field, class _Rep > 778 class BlasMatrix< _Field, _Rep >::ConstRowIterator { 779 public: 780 ConstRowIterator (const typename Rep::const_iterator& p, size_t len, size_t d) : 781 _row (p, p + (ptrdiff_t)len), _dis (d) 782 {} 783 784 ConstRowIterator () {} 785 786 ConstRowIterator (const ConstRowIterator& colp) : 787 _row (colp._row), _dis (colp._dis) 788 {} 789 790 ConstRowIterator& operator = (const ConstRowIterator& colp) 791 { 792 _row = colp._row; 793 _dis = colp._dis; 794 return *this; 795 } 796 797 ConstRowIterator& operator --() 798 { 799 _row = ConstRow (_row.begin () - (ptrdiff_t)_dis, _row.end () - (ptrdiff_t)_dis); 800 return *this; 801 } 802 803 ConstRowIterator operator-- (int) 804 { 805 ConstRowIterator tmp (*this); 806 --*this; 807 return tmp; 808 } 809 810 811 ConstRowIterator& operator++ () 812 { 813 _row = ConstRow (_row.begin () + (ptrdiff_t)_dis, _row.end () + (ptrdiff_t) _dis); 814 return *this; 815 } 816 817 ConstRowIterator operator++ (int) 818 { 819 ConstRowIterator tmp (*this); 820 ++*this; 821 return tmp; 822 } 823 824 ConstRowIterator operator+ (int i) 825 { 826 return ConstRowIterator (_row.begin () + (ptrdiff_t)((int)_dis * i), _row.size (), _dis); 827 } 828 829 ConstRowIterator& operator += (int i) 830 { 831 _row = ConstRow (_row.begin () + (ptrdiff_t)((int)_dis * i), _row.end () + (ptrdiff_t)((int)_dis * i)); 832 return *this; 833 } 834 835 ConstRow operator[] (int i) const 836 { 837 return ConstRow (_row.begin () + (ptrdiff_t)((int)_dis * i), _row.end () + (ptrdiff_t)((int)_dis * i)); 838 } 839 840 ConstRow* operator-> () 841 { 842 return &_row; 843 } 844 845 ConstRow& operator* () 846 { 847 return _row; 848 } 849 850 bool operator!= (const ConstRowIterator& c) const 851 { 852 return (_row.begin () != c._row.begin ()) || (_row.end () != c._row.end ()) || (_dis != c._dis); 853 } 854 855 private: 856 ConstRow _row; 857 size_t _dis; 858 }; 859 860 template < class _Field, class _Rep > 861 class BlasMatrix< _Field, _Rep >::RowIterator { 862 public: 863 RowIterator (const typename Rep::iterator& p, size_t len, size_t d) : 864 _row (p, p + (ptrdiff_t)len), _dis (d) 865 {} 866 867 RowIterator () {} 868 869 RowIterator (const RowIterator& colp) : 870 _row (colp._row), _dis (colp._dis) 871 {} 872 873 RowIterator& operator = (const RowIterator& colp) 874 { 875 _row = colp._row; 876 _dis = colp._dis; 877 return *this; 878 } 879 880 RowIterator& operator ++ () 881 { 882 _row = Row (_row.begin () + (ptrdiff_t)_dis, _row.end () + (ptrdiff_t)_dis); 883 return *this; 884 } 885 886 RowIterator operator ++ (int) 887 { 888 RowIterator tmp (*this); 889 ++*this; 890 return tmp; 891 } 892 893 RowIterator& operator -- () 894 { 895 _row = Row (_row.begin () - (ptrdiff_t)_dis, _row.end () - (ptrdiff_t)_dis); 896 return *this; 897 } 898 899 RowIterator operator -- (int) 900 { 901 RowIterator tmp (*this); 902 --*this; 903 return tmp; 904 } 905 906 RowIterator operator + (int i) 907 { 908 return RowIterator (_row.begin () + (ptrdiff_t)((int)_dis * i), _row.size (), _dis); 909 } 910 911 RowIterator& operator += (int i) 912 { 913 _row = Row (_row.begin () + (ptrdiff_t)((int)_dis * i), _row.end () + (ptrdiff_t)((int)_dis * i)); 914 return *this; 915 } 916 917 Row operator[] (int i) const 918 { 919 return Row (const_cast<Row&> (_row).begin () + (ptrdiff_t)((int)_dis * i), 920 const_cast<Row&> (_row).end () + (ptrdiff_t)((int)_dis * i)); 921 } 922 923 Row* operator-> () 924 { 925 return &(this->_row); 926 } 927 928 Row& operator* () 929 { 930 return _row; 931 } 932 933 bool operator!= (const RowIterator& c) const 934 { 935 return (_row.begin () != c._row.begin ()) || (_row.end () != c._row.end ()) || (_dis != c._dis); 936 } 937 938 operator ConstRowIterator () 939 { 940 return ConstRowIterator (_row.begin (), _row.size (), _dis); 941 } 942 943 private: 944 Row _row; 945 size_t _dis; 946 }; 947 948 template < class _Field, class _Rep > 949 class BlasMatrix< _Field, _Rep >::ConstColIterator { 950 public: 951 ConstColIterator (typename Rep::const_iterator p, size_t stride, size_t len) : 952 _col (Subiterator<typename Rep::const_iterator> (p, (ptrdiff_t)stride), 953 Subiterator<typename Rep::const_iterator> (p + (ptrdiff_t)(len * stride), (ptrdiff_t)stride)), 954 _stride (stride) 955 {} 956 957 ConstColIterator (const ConstCol& col, size_t stride) : 958 _col (col), 959 _stride (stride) 960 {} 961 962 ConstColIterator () {} 963 964 ConstColIterator (const ConstColIterator& rowp) : 965 _col (rowp._col), 966 _stride (rowp._stride) 967 {} 968 969 ConstColIterator& operator= (const ConstColIterator& rowp) 970 { 971 _col = rowp._col; 972 _stride = rowp._stride; 973 return *this; 974 } 975 976 ConstColIterator& operator++ () 977 { 978 _col = ConstCol (Subiterator<typename Rep::const_iterator> (_col.begin ().operator-> () + 1, (ptrdiff_t)_stride), 979 Subiterator<typename Rep::const_iterator> (_col.end ().operator-> () + 1, (ptrdiff_t)_stride)); 980 return *this; 981 } 982 983 ConstColIterator operator++ (int) 984 { 985 ConstColIterator old(*this); 986 this->operator++ (); 987 return old; 988 } 989 990 ConstColIterator operator + (int i) 991 { 992 return ConstColIterator (_col.begin ().operator-> () + i, _stride, _col.size ()); 993 } 994 995 ConstColIterator& operator += (int i) 996 { 997 _col = ConstCol (Subiterator<typename Rep::const_iterator> (_col.begin ().operator-> () + i, _stride), 998 Subiterator<typename Rep::const_iterator> (_col.end ().operator-> () + i, _stride)); 999 return *this; 1000 } 1001 1002 ConstCol operator[] (int i) const 1003 { 1004 return ConstCol (Subiterator<typename Rep::const_iterator> (_col.begin ().operator-> () + i, _stride), 1005 Subiterator<typename Rep::const_iterator> (_col.end ().operator-> () + i, _stride)); 1006 } 1007 1008 ConstCol* operator-> () 1009 { 1010 return &_col; 1011 } 1012 1013 ConstCol& operator* () 1014 { 1015 return _col; 1016 } 1017 1018 bool operator!= (const ConstColIterator& c) const 1019 { 1020 return (_col.begin () != c._col.begin ()) || (_col.end () != c._col.end ()); 1021 } 1022 1023 private: 1024 ConstCol _col; 1025 size_t _stride; 1026 }; 1027 1028 template < class _Field, class _Rep > 1029 class BlasMatrix< _Field, _Rep >::ColIterator { 1030 public: 1031 ColIterator (typename Rep::iterator p, size_t stride, size_t len) : 1032 _col (Subiterator<typename Rep::iterator> (p, (long)stride), 1033 Subiterator<typename Rep::iterator> (p + (ptrdiff_t)(len * stride),(long) stride)), _stride (stride) 1034 {} 1035 1036 ColIterator () {} 1037 1038 ColIterator (const ColIterator& rowp) : 1039 _col (rowp._col) 1040 {} 1041 1042 ColIterator& operator= (const ColIterator& rowp) 1043 { 1044 _col = rowp._col; 1045 _stride = rowp._stride; 1046 return *this; 1047 } 1048 1049 const ColIterator& operator= (const ColIterator& rowp) const 1050 { 1051 const_cast<ColIterator*> (this)->_col = rowp._col; 1052 return *this; 1053 } 1054 1055 ColIterator& operator++ () 1056 { 1057 _col = Col (Subiterator<typename Rep::iterator> (_col.begin ().operator-> () + 1, (const long)_stride), 1058 Subiterator<typename Rep::iterator> (_col.end ().operator-> () + 1,(const long) _stride)); 1059 return *this; 1060 } 1061 1062 ColIterator operator++ (int) 1063 { 1064 Col tmp (_col); 1065 this->operator++ (); 1066 return tmp; 1067 } 1068 1069 ColIterator operator + (int i) 1070 { 1071 return ColIterator (_col.begin ().operator-> () + i, _stride, _col.size ()); 1072 } 1073 1074 ColIterator& operator += (int i) 1075 { 1076 _col = Col (Subiterator<typename Rep::iterator> (_col.begin ().operator-> () + i, _stride), 1077 Subiterator<typename Rep::iterator> (_col.end ().operator-> () + i, _stride)); 1078 return *this; 1079 } 1080 1081 Col operator[] (int i) const 1082 { 1083 return Col (Subiterator<typename Rep::iterator> (const_cast<Col&> (_col).begin ().operator-> () + i, _stride), 1084 Subiterator<typename Rep::iterator> (const_cast<Col&> (_col).end ().operator-> () + i, _stride)); 1085 } 1086 1087 Col* operator-> () 1088 { 1089 return &(this->_col); 1090 } 1091 1092 Col& operator* () 1093 { 1094 return _col; 1095 } 1096 1097 bool operator!= (const ColIterator& c) const 1098 { 1099 return (_col.begin () != c._col.begin ()) || (_col.end () != c._col.end ()); 1100 } 1101 1102 operator ConstColIterator () 1103 { 1104 return ConstColIterator (reinterpret_cast<ConstCol&> (_col) , _stride); 1105 } 1106 1107 private: 1108 1109 Col _col; 1110 size_t _stride; 1111 }; 1112 1113 /*! Indexed Iterator. 1114 * @ingroup iterators 1115 * @brief NO DOC 1116 */ 1117 template < class _Field, class _Rep > 1118 class BlasMatrix< _Field, _Rep >::IndexedIterator { 1119 size_t _r_index; 1120 size_t _c_index; 1121 size_t _dim; 1122 typename Rep::iterator _begin; 1123 typedef typename _Field::Element value_type; 1124 1125 public: 1126 IndexedIterator (const size_t &dim, 1127 const size_t &r_index, 1128 const size_t &c_index, 1129 const typename Rep::iterator &begin) : 1130 _r_index (r_index), _c_index (c_index), _dim (dim), _begin (begin) 1131 {} 1132 1133 IndexedIterator () : 1134 _r_index (0), _c_index (0), _dim (1), _begin (0) 1135 {} 1136 1137 IndexedIterator (const IndexedIterator& r) : 1138 _r_index (r._r_index), _c_index (r._c_index), _dim (r._dim), _begin (r._begin) 1139 {} 1140 1141 IndexedIterator& operator = (const IndexedIterator &iter) 1142 { 1143 _r_index = iter._r_index; 1144 _c_index = iter._c_index; 1145 _dim = iter._dim; 1146 _begin = iter._begin; 1147 return *this; 1148 } 1149 1150 bool operator == (const IndexedIterator &iter) const 1151 { 1152 return (_r_index == iter._r_index) && 1153 (_c_index == iter._c_index) && 1154 (_dim == iter._dim) && 1155 (_begin==iter._begin); 1156 } 1157 1158 bool operator != (const IndexedIterator& iter) const 1159 { 1160 return (_r_index != iter._r_index) || 1161 (_c_index != iter._c_index) || 1162 (_dim != iter._dim) || 1163 (_begin!=iter._begin); 1164 } 1165 1166 IndexedIterator &operator ++ () 1167 { 1168 ++_c_index; 1169 1170 if (_c_index == _dim) { 1171 _c_index = 0; 1172 ++_r_index; 1173 } 1174 1175 return *this; 1176 } 1177 1178 1179 IndexedIterator operator ++ (int) 1180 { 1181 IndexedIterator tmp = *this; 1182 ++(*this); 1183 return tmp; 1184 } 1185 1186 IndexedIterator &operator -- () 1187 { 1188 if (_c_index) 1189 --_c_index; 1190 else { 1191 --_r_index; 1192 _c_index = _dim - 1; 1193 } 1194 1195 return *this; 1196 } 1197 1198 1199 IndexedIterator operator -- (int) 1200 { 1201 IndexedIterator tmp = *this; 1202 --(*this); 1203 return tmp; 1204 } 1205 1206 value_type &operator * () const 1207 { 1208 return *(_begin +(ptrdiff_t) (_r_index * _dim + _c_index)); 1209 } 1210 1211 1212 value_type * operator -> () const 1213 { 1214 return _begin + (ptrdiff_t)(_r_index * _dim + _c_index); 1215 } 1216 1217 1218 size_t rowIndex () const 1219 { 1220 return _r_index; 1221 } 1222 1223 size_t colIndex () const 1224 { 1225 return _c_index; 1226 } 1227 1228 const value_type &value () const 1229 { 1230 return *(_begin + (ptrdiff_t)(_r_index * _dim + _c_index)); 1231 } 1232 1233 1234 }; 1235 1236 template < class _Field, class _Rep > 1237 class BlasMatrix< _Field, _Rep >::ConstIndexedIterator { 1238 size_t _r_index; 1239 size_t _c_index; 1240 size_t _dim; 1241 typedef typename _Field::Element value_type; 1242 typename Rep::const_iterator _begin; 1243 1244 public: 1245 ConstIndexedIterator (const size_t &my_dim, 1246 const size_t &r_index, 1247 const size_t &c_index, 1248 const typename Rep::const_iterator &begin) : 1249 _r_index (r_index), _c_index (c_index), _dim (my_dim), _begin (begin) 1250 {} 1251 1252 ConstIndexedIterator () : 1253 _r_index (0), _c_index (0), _dim (1), _begin (0) 1254 {} 1255 1256 ConstIndexedIterator (const ConstIndexedIterator& r) : 1257 _r_index (r._r_index), _c_index (r._c_index), _dim (r._dim), _begin (r._begin) 1258 {} 1259 1260 ConstIndexedIterator& operator = (const ConstIndexedIterator &iter) 1261 { 1262 _r_index = iter._r_index; 1263 _c_index = iter._c_index; 1264 _dim = iter._dim; 1265 _begin = iter._begin; 1266 return *this; 1267 } 1268 1269 bool operator == (const ConstIndexedIterator &iter) const 1270 { 1271 return (_r_index == iter._r_index) && 1272 (_c_index == iter._c_index) && 1273 (_dim == iter._dim) && 1274 (_begin==iter._begin); 1275 } 1276 1277 bool operator != (const ConstIndexedIterator& iter) const 1278 { 1279 return (_r_index != iter._r_index) || 1280 (_c_index != iter._c_index) || 1281 (_dim != iter._dim) || 1282 (_begin!=iter._begin); 1283 } 1284 1285 ConstIndexedIterator &operator ++ () 1286 { 1287 ++_c_index; 1288 1289 if (_c_index == _dim) { 1290 _c_index = 0; 1291 ++_r_index; 1292 } 1293 1294 return *this; 1295 } 1296 1297 1298 ConstIndexedIterator operator ++ (int) 1299 { 1300 ConstIndexedIterator tmp = *this; 1301 ++(*this); 1302 return tmp; 1303 } 1304 1305 ConstIndexedIterator &operator -- () 1306 { 1307 if (_c_index) 1308 --_c_index; 1309 else { 1310 --_r_index; 1311 _c_index = _dim - 1; 1312 } 1313 1314 return *this; 1315 } 1316 1317 ConstIndexedIterator operator -- (int) 1318 { 1319 ConstIndexedIterator tmp = *this; 1320 --(*this); 1321 return tmp; 1322 } 1323 1324 const value_type &operator * () const 1325 { 1326 return *(_begin + (ptrdiff_t)(_r_index * _dim + _c_index)); 1327 } 1328 1329 const value_type *operator -> () const 1330 { 1331 return _begin + (ptrdiff_t)(_r_index * _dim + _c_index); 1332 } 1333 1334 size_t rowIndex () const 1335 { 1336 return _r_index; 1337 } 1338 1339 size_t colIndex () const 1340 { 1341 return _c_index; 1342 } 1343 1344 const value_type &value() const 1345 { 1346 return *(_begin + (ptrdiff_t)(_r_index * _dim + _c_index)); 1347 } 1348 }; 1349 1350 /* */ 1351 1352 // Entry access view. Size m*n vector in C (row major) order. 1353 template < class _Field, class _Rep > 1354 typename BlasMatrix< _Field, _Rep >::Iterator BlasMatrix< _Field, _Rep >::Begin () 1355 { 1356 return _rep.begin (); 1357 } 1358 1359 template < class _Field, class _Rep > 1360 typename BlasMatrix< _Field, _Rep >::Iterator BlasMatrix< _Field, _Rep >::End () 1361 { 1362 return _rep.end (); 1363 } 1364 1365 template < class _Field, class _Rep > 1366 typename BlasMatrix< _Field, _Rep >::ConstIterator BlasMatrix< _Field, _Rep >::Begin () const 1367 { 1368 return _rep.begin (); 1369 } 1370 1371 template < class _Field, class _Rep > 1372 typename BlasMatrix< _Field, _Rep >::ConstIterator BlasMatrix< _Field, _Rep >::End () const 1373 { 1374 return _rep.end (); 1375 } 1376 1377 /* Indexed */ 1378 1379 template < class _Field, class _Rep > 1380 typename BlasMatrix< _Field, _Rep >::IndexedIterator BlasMatrix< _Field, _Rep >::IndexedBegin () 1381 { 1382 return IndexedIterator (coldim (), 0, 0, _rep.begin ()); 1383 } 1384 1385 template < class _Field, class _Rep > 1386 typename BlasMatrix< _Field, _Rep >::IndexedIterator BlasMatrix< _Field, _Rep >::IndexedEnd () 1387 { 1388 return IndexedIterator (coldim (), rowdim (), 0, _rep.begin ()); 1389 } 1390 1391 template < class _Field, class _Rep > 1392 typename BlasMatrix< _Field, _Rep >::ConstIndexedIterator BlasMatrix< _Field, _Rep >::IndexedBegin () const 1393 { 1394 return ConstIndexedIterator (coldim (), 0, 0, _rep.begin ()); 1395 } 1396 1397 template < class _Field, class _Rep > 1398 typename BlasMatrix< _Field, _Rep >::ConstIndexedIterator BlasMatrix< _Field, _Rep >::IndexedEnd () const 1399 { 1400 return ConstIndexedIterator (coldim (), rowdim (), 0, _rep.begin ()); 1401 } 1402 1403 /* Row */ 1404 1405 template < class _Field, class _Rep > 1406 typename BlasMatrix< _Field, _Rep >::RowIterator BlasMatrix< _Field, _Rep >::rowBegin () 1407 { 1408 return RowIterator (_rep.begin (), _col, _col); 1409 } 1410 1411 template < class _Field, class _Rep > 1412 typename BlasMatrix< _Field, _Rep >::RowIterator BlasMatrix< _Field, _Rep >::rowEnd () 1413 { 1414 return RowIterator (_rep.end (), _col, _col); 1415 } 1416 1417 template < class _Field, class _Rep > 1418 typename BlasMatrix< _Field, _Rep >::ConstRowIterator BlasMatrix< _Field, _Rep >::rowBegin () const 1419 { 1420 return ConstRowIterator (_rep.begin (), _col, _col); 1421 } 1422 1423 template < class _Field, class _Rep > 1424 typename BlasMatrix< _Field, _Rep >::ConstRowIterator BlasMatrix< _Field, _Rep >::rowEnd () const 1425 { 1426 return ConstRowIterator (_rep.end (), _col, _col); 1427 } 1428 1429 /* Col */ 1430 1431 template < class _Field, class _Rep > 1432 typename BlasMatrix< _Field, _Rep >::ColIterator BlasMatrix< _Field, _Rep >::colBegin () 1433 { 1434 return typename BlasMatrix< _Field, _Rep >::ColIterator (_rep.begin (), _col, _row); 1435 } 1436 1437 template < class _Field, class _Rep > 1438 typename BlasMatrix< _Field, _Rep >::ColIterator BlasMatrix< _Field, _Rep >::colEnd () 1439 { 1440 return typename BlasMatrix< _Field, _Rep >::ColIterator (_rep.begin ()+(ptrdiff_t)_col, _col, _row); 1441 } 1442 1443 template < class _Field, class _Rep > 1444 typename BlasMatrix< _Field, _Rep >::ConstColIterator BlasMatrix< _Field, _Rep >::colBegin () const 1445 { 1446 return typename BlasMatrix< _Field, _Rep >::ConstColIterator (_rep.begin (), _col, _row); 1447 } 1448 1449 template < class _Field, class _Rep > 1450 typename BlasMatrix< _Field, _Rep >::ConstColIterator BlasMatrix< _Field, _Rep >::colEnd () const 1451 { 1452 return typename BlasMatrix< _Field, _Rep >::ConstColIterator (_rep.begin ()+(ptrdiff_t)_col, _col, _row); 1453 } 1454 1455 /* operators */ 1456 template < class _Field, class _Rep > 1457 typename BlasMatrix< _Field, _Rep >::Row BlasMatrix< _Field, _Rep >::operator[] (size_t i) 1458 { 1459 return Row (_rep.begin () +(ptrdiff_t)( i * _col), _rep.begin () + (ptrdiff_t)(i * _col +_col)); 1460 } 1461 1462 template < class _Field, class _Rep > 1463 typename BlasMatrix< _Field, _Rep >::ConstRow BlasMatrix< _Field, _Rep >::operator[] (size_t i) const 1464 { 1465 return Row (_rep.begin () +(ptrdiff_t) (i * _col), _rep.begin () + (ptrdiff_t)( i * _col + _col)); 1466 } 1467 1468} // LinBox 1469 1470////////////////// 1471// MISC // 1472////////////////// 1473 1474namespace LinBox 1475{ 1476 template < class _Field, class _Rep > 1477 template <class Vector> 1478 Vector& BlasMatrix< _Field, _Rep >::columnDensity (Vector &v) const 1479 { 1480 std::fill (v.begin (), v.end (), _row); 1481 return v; 1482 } 1483 1484} // LinBox 1485 1486/////////////////// 1487// BLACK BOX // 1488/////////////////// 1489 1490namespace LinBox 1491{ 1492#if 1 1493 template < class _Field, class _Rep > 1494 template <class Vector1, class Vector2> 1495 Vector1& BlasMatrix< _Field, _Rep >::apply (Vector1& y, const Vector2& x) const 1496 { 1497 // std::cout << "prepare apply1 Matrix" << std::endl; 1498 BlasSubmatrix<constSelf_t> A(*this); 1499 // std::cout << "...................." << std::endl; 1500 A.apply(y,x); 1501 // std::cout << ".......done........." << std::endl; 1502 return y; 1503 } 1504#endif 1505 1506#if 1 1507 template < class _Field, class _Rep > 1508 template<class _Vrep> 1509 BlasVector<_Field,_Vrep>& BlasMatrix< _Field, _Rep >::apply (BlasVector<_Field,_Vrep>& y, const BlasVector<_Field,_Vrep>& x) const 1510 { 1511 // std::cout << "prepare apply2 Matrix" << std::endl; 1512 BlasSubmatrix<constSelf_t> A(*this); 1513 // std::cout << "...................." << std::endl; 1514 A.apply(y,x); 1515 // std::cout << ".......done........." << std::endl; 1516 1517 return y; 1518 } 1519#endif 1520 1521 template < class _Field, class _Rep > 1522 template <class Vector1, class Vector2> 1523 Vector1& BlasMatrix< _Field, _Rep >::applyTranspose (Vector1& y, const Vector2& x) const 1524 { 1525 // std::cout << "prepare applyT Matrix" << std::endl; 1526 BlasSubmatrix<constSelf_t> A(*this); 1527 // std::cout << "...................." << std::endl; 1528 A.applyTranspose(y,x); 1529 // std::cout << ".......done........." << std::endl; 1530 1531 return y; 1532 } 1533 1534 template < class _Field, class _Rep > 1535 const _Field& BlasMatrix< _Field, _Rep >::field() const 1536 { 1537 return *_field; 1538 } 1539 1540#if 0 /* why not ? */ 1541 template < class _Field, class _Rep > 1542 _Field& BlasMatrix< _Field, _Rep >::field() 1543 { 1544 return const_cast<_Field&>(_field ); 1545 } 1546#endif 1547} 1548 1549 1550#endif // __LINBOX_densematrix_blas_matrix_INL 1551 1552// Local Variables: 1553// mode: C++ 1554// tab-width: 4 1555// indent-tabs-mode: nil 1556// c-basic-offset: 4 1557// End: 1558// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 1559