1 /* linbox/matrix/dense-matrix.h 2 * Copyright (C) 2011 B. David Saunders, 3 * See COPYING for license information 4 * 5 * evolved from dense-submatrix and blas-matrix by B. David Saunders <saunders@cis.udel.edu>, BSY 6 */ 7 8 /*! @file matrix/sliced3/dense-submatrix.h 9 * @ingroup matrix 10 * @brief Representation of a submatrix of a dense matrix, not resizeable. 11 * This matrix type conforms to the \c LinBox::DenseMat interface. 12 * \c LinBox::BlasMatrix is an example of DenseSubmatrix. 13 */ 14 15 #ifndef __LINBOX_dense_matrix_H 16 #define __LINBOX_dense_matrix_H 17 18 #include <utility> 19 #include "linbox/linbox-config.h" 20 21 //#include "linbox/util/debug.h" 22 //#include "linbox/matrix/matrix-domain.h" 23 #include "linbox/vector/subvector.h" 24 #include "submat-iterator.h" 25 26 namespace LinBox 27 { 28 29 /** @brief to be used in standard matrix domain 30 31 * Matrix variable declaration, sizing, entry initialization may involve one to 3 steps. 32 * Matrix ops are container ops. (sizing, copying) 33 * Mathematically meaningful operations are to be found only in an associated matrix domain 34 * 35 * A matrix may be allocated or not. A matrix initialized by a submatrix() call is not allocated. 36 * When an allocated matrix goes out of scope or is reinitialized with init(), the memory is released 37 * and all submatrices of it become invalid. 38 * 39 * Allocating: 40 * DenseMat A(2, 3); // allocation of mem for 6 entries at construction 41 * DenseMat B; B.init(10, 10); // default constr and subsequent allocation for 100 entries. 42 * 43 * Allocation of memory plus entry initialization: 44 * // a meaningful value of DenseMat::Entry x is set by a field or matrix domain. 45 * DenseMat A(10, 10, x); 46 * DenseMat B: B.init(10, 10, x); 47 * DenseMat C(A); // allocation at copy construction. A could be a submatrix of another. 48 * A.read(istream) 49 * 50 * Nonallocation sizing: 51 * // assume D declared, A initialized, n = A.coldim(). 52 * D.submatrix(A, 0, 1, 2, n-1); // D is 2 by n-1 in upper right of A. 53 * 54 * Entry initialization (and overwriting) in already sized matrices: 55 * A.setEntry(i, j, x); 56 * A = B; // A and B must have the same shape. 57 58 * Entry read access. OK on const matrices 59 * getEntry, write 60 61 * Under consideration: Require A.clear() on an allocated matrix before any action that would abandon 62 * the allocated mem (init or submatrix). 63 \ingroup matrix 64 */ 65 template<class _Element> 66 class DenseMat { 67 typedef _Element Entry; 68 public: // protected: //TODO scope 69 Entry *_rep; // matrix entries on the heap. 70 bool _alloc; // iff _alloc, I am responsible for _rep allocation/deallocation. 71 size_t _rows; 72 size_t _cols; 73 size_t _stride; // stride from row to row. Entries in a row are contiguous. 74 75 public: rowdim()76 size_t rowdim() const { return _rows; } coldim()77 size_t coldim() const { return _cols; } 78 getEntry(Entry & x,size_t i,size_t j)79 Entry& getEntry(Entry& x, size_t i, size_t j) const { 80 linbox_check(i < _rows && j < _cols); 81 return x = *(_rep + i*_stride + j); 82 } 83 setEntry(size_t i,size_t j,const Entry & x)84 Entry& setEntry(size_t i, size_t j, const Entry& x ) { 85 //linbox_check((i < _rows && j < _cols)); 86 return *(_rep + i*_stride + j) = x; 87 } 88 // no refEntry - not well supportable in some matrix domains. 89 DenseMat()90 DenseMat() : _rep(NULL), _alloc(false), _rows(0), _cols(0), _stride(0) {} 91 92 void init(size_t m = 0, size_t n = 0) { 93 //std::cerr << m << " " << n << " <<<<<<" << std::endl; 94 if (_alloc) delete _rep; // abandon any prior def 95 if (m*n != 0) { 96 _rep = new Entry[m*n]; _alloc = true; 97 } else { 98 _rep = NULL; _alloc = false; 99 } 100 size(m, n); 101 } 102 init(size_t m,size_t n,const Entry & filler)103 void init(size_t m, size_t n, const Entry& filler) { 104 if (_alloc) delete _rep; // abandon any prior def 105 if (m*n != 0) { 106 _rep = new Entry[m*n]; _alloc = true; 107 for (size_t i = 0; i < m*n; ++i) _rep[i] = filler; 108 } else { 109 _rep = NULL; _alloc = false; 110 } 111 size(m, n); 112 } 113 114 /* Copy construction makes this a completely distinct copy of A. 115 * Entries of this are stored contiguously even if A is not contiguous (is a submatrix of another). 116 * 117 * If memcpy is not valid for your entry type, specialize DenseMat for it. 118 */ DenseMat(const DenseMat & A)119 DenseMat(const DenseMat& A) 120 : _rep(new Entry[A._rows*A._cols]), _alloc(true), _rows(A._rows), _cols(A._cols), _stride(A._cols) 121 { 122 //std::cout << "copy construction " << _rep << std::endl; 123 //std::cout << "copy cons " << _rows << " " << _cols << std::endl; 124 *this = A; // copy A's data into _rep 125 } 126 ~DenseMat()127 ~DenseMat() { 128 if (_alloc) delete _rep; 129 } 130 131 // For assignment, the matrices must have same size and not overlap in memory. 132 // This restriction could be loosened... 133 DenseMat& operator=(const DenseMat& B) { 134 linbox_check(_rows == B._rows && _cols == B._cols); 135 if (_cols == _stride && B._cols == B._stride) // both are contiguous 136 memcpy(_rep, B._rep, sizeof(Entry)*_rows*_cols); 137 else 138 for (size_t i = 0; i < _rows; ++i) // copy row by row 139 memcpy(_rep + i*_stride, B._rep + i*B._stride, sizeof(Entry)*_cols); 140 return *this; 141 } 142 143 /** Set this to be an m by n submatrix of A with upper left corner at i,j position of A. 144 * Requires i+m <= A.rowdim(), j+n <= A.coldim(). 145 * 146 * For instance, B.submatrix(A, i, 0, 1, A.coldim()) makes B the i-th row of A. 147 */ submatrix(const DenseMat & A,size_t i,size_t j,size_t m,size_t n)148 void submatrix(const DenseMat & A, size_t i, size_t j, size_t m, size_t n) { 149 linbox_check(i+m <= A._rows && j+n <= A._cols); 150 if (_alloc) delete _rep; // abandon any prior def 151 _rep = A._rep + (i*A._stride + j); 152 _alloc = false; 153 _rows = m; _cols = n; _stride = A._stride; 154 } 155 156 /* SLICED helper functions... maybe shouldn't be public */ 157 158 // size a matrix w/o allocating memory size(size_t m,size_t n)159 void size(size_t m, size_t n){ 160 _rows = m; _cols = _stride = n; 161 } 162 163 template<class vp> did_swapback(vp & swaps,size_t hot,size_t l,size_t r)164 int did_swapback(vp &swaps, size_t hot, size_t l, size_t r){ 165 for(size_t i=l; i<r; ++i){ 166 if(swaps[i].second == hot) 167 return swaps[i].first; 168 } 169 return -1; 170 } 171 172 template<class vp> flocs(vp & tos,vp & swaps)173 void flocs(vp & tos, vp & swaps){ 174 std::vector<size_t> pos; 175 for(size_t i =0; i < coldim(); ++i) pos.push_back(i); 176 177 for(int i = coldim()-1; i>= 0; --i){ 178 //std::cerr << swaps[i] << std::endl; 179 swap(pos[swaps[i].first], pos[swaps[i].second]); 180 } 181 182 for(size_t i = 0; i < coldim(); ++i){ 183 std::pair<size_t, size_t> p(i, pos[i]); 184 tos.push_back(p); 185 } 186 } 187 188 template<class vp> final_loc(vp & swaps,size_t j)189 size_t final_loc(vp &swaps, size_t j){ 190 // first check if we were swapped back prior to being reached 191 int to = did_swapback(swaps, j, 0, j+1); 192 if(to >= 0) return to; 193 // next check if we get swapped back from our immediate neightbor 194 //size_t nbor = swaps[j].second; 195 //to = did_swapback(swaps, nbor, j+1, nbor); 196 //if(to >= 0) return to; 197 // bad luck, we are chained forward, so 198 size_t i = j; 199 while(to < 0){ 200 size_t lc = swaps[i].first; 201 size_t rc = swaps[i].second; 202 to = did_swapback(swaps, rc, lc+1, rc+1); 203 //std::cerr << j << ": " << to << "= ds(swps, " << rc << "," << lc+1 << "," << rc+1 << ")" << std::endl; 204 205 i = rc; 206 if(to < 0 && rc == swaps.size() - 1){ 207 to = swaps.size() - 1; 208 } 209 } 210 return to; 211 } 212 randomColPermutation()213 void randomColPermutation() { 214 typedef std::pair<size_t, size_t> pair; 215 typedef std::vector<pair> vp; 216 vp swaps; 217 vp tos; 218 vp tos2; 219 // create pairs 220 for (size_t j = 0; j < coldim(); ++j){ 221 // Each iteration swap col j with a random col in range [j..n-1]. 222 int k = j + rand()%(coldim()-j); 223 //std::cerr << j << "->" << k << std::endl; 224 pair p(j,k); 225 swaps.push_back(p); 226 //for (size_t i = 0; i < rowdim(); ++i) 227 //swap( _rep[i*_stride + j], _rep[i*_stride + k]); 228 } 229 /* 230 //flocs(tos2, swaps); 231 // find each final location 232 for(size_t j = 0; j < coldim(); ++j){ 233 size_t to = final_loc(swaps, j); 234 pair p(j, to); 235 if(p != tos2[j]){ 236 std::cerr << "HOUSTON: PROBLEM" << std::endl; 237 } 238 else{ 239 std::cout << "YAY "; 240 } 241 tos.push_back(p); 242 } 243 */ 244 flocs(tos, swaps); 245 // tos is the correct mapping now permute 246 Entry *perm_row = new Entry[coldim()]; 247 for(size_t i = 0; i < rowdim(); ++i){ 248 for(vp::iterator ti = tos.begin(); ti!= tos.end(); ++ti) 249 perm_row[(*ti).second] = _rep[i*_stride + (*ti).first]; 250 memcpy(&(_rep[i*_stride]), perm_row, coldim()*sizeof(Entry)); 251 } 252 delete[] perm_row; 253 254 // CHECKING CODE 255 /* 256 std::cerr << "swaps: "; 257 for(vp::iterator i = swaps.begin(); i!= swaps.end(); ++i){ 258 std::cerr << *i << " "; 259 } 260 std::cerr << std::endl; 261 std::cerr << "tos: "; 262 for(vp::iterator i = tos.begin(); i!= tos.end(); ++i){ 263 std::cerr << *i << " "; 264 } 265 std::cerr << std::endl; 266 std::vector<int> orig; 267 for(size_t i = 0; i < coldim(); ++i) orig.push_back((int)i); 268 std::vector<int> out(orig); 269 std::vector<int> out2(orig); 270 for(vp::iterator i = swaps.begin(); i!= swaps.end(); ++i){ 271 std::swap(out[(*i).first], out[(*i).second]); 272 } 273 for(vp::iterator i = tos.begin(); i!= tos.end(); ++i){ 274 out2[(*i).second] = orig[(*i).first]; 275 } 276 for(std::vector<int>::iterator i=out.begin(), i2=out2.begin(); i!=out.end(); ++i, ++i2){ 277 if(*i != *i2){ 278 std::cerr << "LOSER: "; 279 std::cerr << *i << " v " << *i2 << std::endl; 280 } 281 } 282 */ 283 } 284 randomLowerTriangularColTransform()285 void randomLowerTriangularColTransform() { 286 for (size_t j = 0; j < coldim(); ++j){ 287 // Each iteration swap col j with a random col in range [j..n-1]. 288 //int l = 1, k = 0; 289 //do { l *= RAND_MAX; k = k*RAND_MAX + rand(); } while (l < j); 290 //k = k%j; 291 int k = rand()%j; 292 for (size_t i = 0; i < rowdim(); ++i) 293 _rep[i*_stride + j] += _rep[i*_stride + k]; 294 } 295 } 296 /* Iterators */ 297 298 typedef SubMatIterator<_Element> RawIterator; 299 typedef ConstSubMatIterator<_Element> ConstRawIterator; 300 rawBegin()301 RawIterator rawBegin() { 302 return RawIterator(_rep, _cols, _stride); }; rawEnd()303 RawIterator rawEnd() { 304 return RawIterator(_rep + _rows*_stride); } rowBegin(size_t i)305 RawIterator rowBegin(size_t i) { 306 //return rawBegin() + i*_cols; 307 return RawIterator(_rep+i*_stride, _cols, _stride); 308 } rowEnd(size_t i)309 RawIterator rowEnd(size_t i) { 310 return rowBegin(i + 1); 311 } 312 /* 313 RawIterator rowBegin(size_t i) { 314 return rawBegin() + i*_cols; } 315 RawIterator rowEnd(size_t i) { 316 return rowBegin(i + 1); } 317 */ 318 rawBegin()319 ConstRawIterator rawBegin() const { 320 return ConstRawIterator(_rep, _cols, _stride); } rawEnd()321 ConstRawIterator rawEnd() const { 322 return ConstRawIterator(_rep + _rows*_stride); } rowBegin(size_t i)323 ConstRawIterator rowBegin(size_t i) const { 324 return rawBegin() + i*_cols; } rowEnd(size_t i)325 ConstRawIterator rowEnd(size_t i) const { 326 return rowBegin(i + 1); } 327 328 typedef DenseMat<_Element> Self_t; //!< Self type 329 330 #if 0 // TODO factor out 331 /** @name typedef'd Row Iterators. 332 *\brief 333 * The row iterator gives the rows of the 334 * matrix in ascending order. Dereferencing the iterator yields 335 * a row vector in dense format 336 * @{ 337 */ 338 typedef Entry* RowIterator; 339 typedef const Entry* ConstRowIterator; 340 typedef Subvector<Entry*> Row; 341 typedef Subvector<const Entry*> ConstRow; 342 //@} Row Iterators 343 344 /** @name typedef'd Column Iterators. 345 *\brief 346 * The columns iterator gives the columns of the 347 * matrix in ascending order. Dereferencing the iterator yields 348 * a column vector in dense format 349 * @{ 350 */ 351 typedef Subiterator<Entry> ColIterator; 352 typedef Subiterator<const Entry> ConstColIterator; 353 typedef Subvector<ColIterator> Col; 354 typedef Subvector<ConstColIterator> ConstCol; 355 //@} // Column Iterators 356 #endif 357 358 template<typename _Tp1> 359 struct rebind { 360 typedef DenseMat<typename _Tp1::Element> other; 361 }; 362 363 /** Read the matrix from an input stream. 364 * @param file Input stream from which to read 365 * @param field 366 */ 367 template<class Field> 368 std::istream& read (std::istream &file, const Field& field); 369 370 /** Write the matrix to an output stream. 371 * @param os Output stream to which to write 372 * @param field 373 * @param mapleFormat write in Maple(r) format ? 374 */ 375 template<class Field> 376 std::ostream& write (std::ostream &os, const Field& field, 377 bool mapleFormat = false) const; 378 379 /** Write the matrix to an output stream. 380 * This a raw version of \c write(os,F) (no field is given). 381 * @param os Output stream to which to write 382 * @param mapleFormat write in Maple(r) format ? 383 */ 384 std::ostream& write (std::ostream &os, 385 bool mapleFormat = false) const; 386 387 /* TODO factor out row/col iterator 388 RowIterator rowBegin (); 389 RowIterator rowEnd (); 390 ConstRowIterator rowBegin () const; 391 ConstRowIterator rowEnd () const; 392 393 ColIterator colBegin (); 394 ColIterator colEnd (); 395 ConstColIterator colBegin () const; 396 ConstColIterator colEnd () const; 397 */ 398 399 }; 400 /*! Write a matrix to a stream. 401 * The C++ way using <code>operator<<</code> 402 * @param o output stream 403 * @param Mat matrix to write. 404 */ 405 template<class T> 406 std::ostream& operator<< (std::ostream & o, const DenseMat<T> & Mat) 407 { 408 return Mat.write(o); 409 } 410 411 /*! @internal 412 * @brief MatrixTraits 413 */ 414 /* necessary? 415 template <class Element> 416 struct MatrixTraits< DenseMat<Element> > { 417 typedef DenseMat<Element> MatrixType; 418 typedef typename MatrixCategories::RowColMatrixTag MatrixCategory; 419 }; 420 */ 421 422 template <class _Element> 423 template <class Field> read(std::istream & file,const Field & field)424 std::istream& DenseMat<_Element>::read (std::istream &file, const Field& field) 425 { 426 RawIterator p; 427 428 for (p = rawBegin (); p != rawEnd (); ++p) { 429 // each entry is seperated by one space. 430 file.ignore (1); 431 field.read (file, *p); 432 } 433 434 return file; 435 } 436 437 template <class _Element> 438 template <class Field> write(std::ostream & os,const Field & field,bool mapleFormat)439 std::ostream &DenseMat<_Element>::write (std::ostream &os, const Field& field, 440 bool mapleFormat) const 441 { 442 return os; 443 } 444 /* TODO factor out RowIterator 445 ConstRowIterator p; 446 447 // integer c; 448 //int wid; 449 450 // field.cardinality (c); 451 //wid = (int) ceil (log ((double) c) / M_LN10); //BB : not used ! 452 453 typename ConstRow::const_iterator pe; 454 455 if (mapleFormat) os << "["; 456 457 for (p = rowBegin (); p != rowEnd (); ++p) { 458 if (mapleFormat && (p != rowBegin())) 459 os << ','; 460 if (mapleFormat) os << "["; 461 462 for (pe = p->begin (); pe != p->end (); ++pe) { 463 if (mapleFormat && (pe != p->begin())) os << ','; 464 // matrix base does not provide this field(), maybe should? 465 //_M.field ().write (os, *pe); 466 //os << *pe; 467 //fixed by using extra field 468 469 field.write (os, *pe); 470 os << " "; 471 } 472 473 if (!mapleFormat) 474 os << std::endl; 475 else os << ']'; 476 } 477 478 if (mapleFormat) os << ']'; 479 return os; 480 } 481 */ 482 483 template <class _Element> write(std::ostream & os,bool mapleFormat)484 std::ostream &DenseMat<_Element>::write (std::ostream &os, bool mapleFormat) const 485 { 486 return os; 487 } 488 /* 489 ConstRowIterator p; 490 491 492 493 typename ConstRow::const_iterator pe; 494 495 if (mapleFormat) os << "["; 496 497 for (p = rowBegin (); p != rowEnd (); ++p) { 498 if (mapleFormat && (p != rowBegin())) 499 os << ','; 500 if (mapleFormat) os << "["; 501 502 for (pe = p->begin (); pe != p->end (); ++pe) { 503 if (mapleFormat && (pe != p->begin())) os << ','; 504 505 os << *pe; 506 os << " "; 507 } 508 509 if (!mapleFormat) 510 os << std::endl; 511 else os << ']'; 512 } 513 514 if (mapleFormat) os << ']'; 515 return os; 516 } 517 */ 518 519 } // namespace LinBox 520 521 #endif // __LINBOX_dense_matrix_H 522 523 // Local Variables: 524 // mode: C++ 525 // tab-width: 4 526 // indent-tabs-mode: nil 527 // c-basic-offset: 4 528 // End: 529 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 530