1 /*
  2  *   Copyright (c) 2007, Michael Lehn
  3  *
  4  *   All rights reserved.
  5  *
  6  *   Redistribution and use in source and binary forms, with or without
  7  *   modification, are permitted provided that the following conditions
  8  *   are met:
  9  *
 10  *   1) Redistributions of source code must retain the above copyright
 11  *      notice, this list of conditions and the following disclaimer.
 12  *   2) Redistributions in binary form must reproduce the above copyright
 13  *      notice, this list of conditions and the following disclaimer in
 14  *      the documentation and/or other materials provided with the
 15  *      distribution.
 16  *   3) Neither the name of the FLENS development group nor the names of
 17  *      its contributors may be used to endorse or promote products derived
 18  *      from this software without specific prior written permission.
 19  *
 20  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 21  *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 22  *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 23  *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 24  *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 25  *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 26  *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 27  *   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 28  *   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 29  *   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 30  *   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 31  */
 32 
 33 #ifndef FLENS_STORAGE_FULLSTORAGE_FULLSTORAGE_TCC
 34 #define FLENS_STORAGE_FULLSTORAGE_FULLSTORAGE_TCC 1
 35 
 36 #include <cxxblas/level1extensions/gecopy.h>
 37 #include <flens/storage/fullstorage/trapezoidalfill.h>
 38 #include <flens/typedefs.h>
 39 
 40 namespace flens {
 41 
 42 //= Constructors
 43 
 44 template <typename T, StorageOrder Order, typename I, typename A>
 45 FullStorage<T, Order, I, A>::FullStorage()
 46     :  _data(0),
 47        _numRows(0), _numCols(0),
 48        _firstRow(I::defaultIndexBase), _firstCol(I::defaultIndexBase)
 49 {
 50 }
 51 
 52 template <typename T, StorageOrder Order, typename I, typename A>
 53 FullStorage<T, Order, I, A>::FullStorage(IndexType numRows, IndexType numCols,
 54                                          IndexType firstRow, IndexType firstCol,
 55                                          const ElementType &value,
 56                                          const Allocator &allocator)
 57     : _data(0), _allocator(allocator),
 58       _numRows(numRows), _numCols(numCols),
 59       _firstRow(firstRow), _firstCol(firstCol)
 60 {
 61     ASSERT(_numRows>=0);
 62     ASSERT(_numCols>=0);
 63 
 64     _allocate(value);
 65 }
 66 
 67 template <typename T, StorageOrder Order, typename I, typename A>
 68 FullStorage<T, Order, I, A>::FullStorage(const FullStorage &rhs)
 69     : _data(0), _allocator(rhs.allocator()),
 70       _numRows(rhs.numRows()), _numCols(rhs.numCols()),
 71       _firstRow(rhs.firstRow()), _firstCol(rhs.firstCol())
 72 {
 73     _allocate(ElementType());
 74     cxxblas::gecopy(Order,
 75                     NoTrans, _numRows, _numCols,
 76                     rhs.data(), rhs.leadingDimension(),
 77                     data(), leadingDimension());
 78 }
 79 
 80 template <typename T, StorageOrder Order, typename I, typename A>
 81 template <typename RHS>
 82 FullStorage<T, Order, I, A>::FullStorage(const RHS &rhs)
 83     : _data(0), _allocator(rhs.allocator()),
 84       _numRows(rhs.numRows()), _numCols(rhs.numCols()),
 85       _firstRow(rhs.firstRow()), _firstCol(rhs.firstCol())
 86 {
 87     _allocate(ElementType());
 88     cxxblas::gecopy(Order,
 89                     NoTrans, _numRows, _numCols,
 90                     rhs.data(), rhs.leadingDimension(),
 91                     data(), leadingDimension());
 92 }
 93 
 94 template <typename T, StorageOrder Order, typename I, typename A>
 95 FullStorage<T, Order, I, A>::~FullStorage()
 96 {
 97     _release();
 98 }
 99 
100 //-- operators -----------------------------------------------------------------
101 
102 template <typename T, StorageOrder Order, typename I, typename A>
103 const typename FullStorage<T, Order, I, A>::ElementType &
104 FullStorage<T, Order, I, A>::operator()(IndexType row, IndexType col) const
105 {
106     ASSERT(row>=_firstRow);
107     ASSERT(row<_firstRow+_numRows);
108     ASSERT(col>=_firstCol);
109     ASSERT(col<_firstCol+_numCols);
110 
111     if (Order==ColMajor) {
112         return _data[col*_numRows+row];
113     }
114     return _data[row*_numCols+col];
115 }
116 
117 template <typename T, StorageOrder Order, typename I, typename A>
118 typename FullStorage<T, Order, I, A>::ElementType &
119 FullStorage<T, Order, I, A>::operator()(IndexType row, IndexType col)
120 {
121     ASSERT(row>=_firstRow);
122     ASSERT(row<_firstRow+_numRows);
123     ASSERT(col>=_firstCol);
124     ASSERT(col<_firstCol+_numCols);
125 
126     if (Order==ColMajor) {
127         return _data[col*_numRows+row];
128     }
129     return _data[row*_numCols+col];
130 }
131 
132 //-- Methods -------------------------------------------------------------------
133 template <typename T, StorageOrder Order, typename I, typename A>
134 typename FullStorage<T, Order, I, A>::IndexType
135 FullStorage<T, Order, I, A>::firstRow() const
136 {
137     return _firstRow;
138 }
139 
140 template <typename T, StorageOrder Order, typename I, typename A>
141 typename FullStorage<T, Order, I, A>::IndexType
142 FullStorage<T, Order, I, A>::lastRow() const
143 {
144     return _firstRow+_numRows-1;
145 }
146 
147 template <typename T, StorageOrder Order, typename I, typename A>
148 typename FullStorage<T, Order, I, A>::IndexType
149 FullStorage<T, Order, I, A>::firstCol() const
150 {
151     return _firstCol;
152 }
153 
154 template <typename T, StorageOrder Order, typename I, typename A>
155 typename FullStorage<T, Order, I, A>::IndexType
156 FullStorage<T, Order, I, A>::lastCol() const
157 {
158     return _firstCol+_numCols-1;
159 }
160 
161 template <typename T, StorageOrder Order, typename I, typename A>
162 typename FullStorage<T, Order, I, A>::IndexType
163 FullStorage<T, Order, I, A>::numRows() const
164 {
165     return _numRows;
166 }
167 
168 template <typename T, StorageOrder Order, typename I, typename A>
169 typename FullStorage<T, Order, I, A>::IndexType
170 FullStorage<T, Order, I, A>::numCols() const
171 {
172     return _numCols;
173 }
174 
175 template <typename T, StorageOrder Order, typename I, typename A>
176 typename FullStorage<T, Order, I, A>::IndexType
177 FullStorage<T, Order, I, A>::leadingDimension() const
178 {
179     return (Order==ColMajor) ? std::max(_numRows, IndexType(1))
180                              : std::max(_numCols, IndexType(1));
181 }
182 
183 template <typename T, StorageOrder Order, typename I, typename A>
184 typename FullStorage<T, Order, I, A>::IndexType
185 FullStorage<T, Order, I, A>::strideRow() const
186 {
187     return (Order==ColMajor) ? 1
188                              : leadingDimension();
189 }
190 
191 template <typename T, StorageOrder Order, typename I, typename A>
192 typename FullStorage<T, Order, I, A>::IndexType
193 FullStorage<T, Order, I, A>::strideCol() const
194 {
195     return (Order==ColMajor) ? leadingDimension()
196                              : 1;
197 }
198 
199 template <typename T, StorageOrder Order, typename I, typename A>
200 const typename FullStorage<T, Order, I, A>::ElementType *
201 FullStorage<T, Order, I, A>::data() const
202 {
203 #   ifndef NDEBUG
204     if ((numRows()==0) || numCols()==0) {
205         return 0;
206     }
207 #   endif
208 
209     ASSERT(_data);
210     return &(this->operator()(_firstRow, _firstCol));
211 }
212 
213 template <typename T, StorageOrder Order, typename I, typename A>
214 typename FullStorage<T, Order, I, A>::ElementType *
215 FullStorage<T, Order, I, A>::data()
216 {
217 #   ifndef NDEBUG
218     if ((numRows()==0) || numCols()==0) {
219         return 0;
220     }
221 #   endif
222 
223     ASSERT(_data);
224     return &(this->operator()(_firstRow, _firstCol));
225 }
226 
227 template <typename T, StorageOrder Order, typename I, typename A>
228 const typename FullStorage<T, Order, I, A>::Allocator &
229 FullStorage<T, Order, I, A>::allocator() const
230 {
231     return _allocator;
232 }
233 
234 template <typename T, StorageOrder Order, typename I, typename A>
235 bool
236 FullStorage<T, Order, I, A>::resize(IndexType numRows, IndexType numCols,
237                                     IndexType firstRow, IndexType firstCol,
238                                     const ElementType &value)
239 {
240     if ((_numRows!=numRows) || (_numCols!=numCols)) {
241         _release();
242         _numRows = numRows;
243         _numCols = numCols;
244         _firstRow = firstRow;
245         _firstCol = firstCol;
246         _allocate(value);
247         return true;
248     }
249     changeIndexBase(firstRow, firstCol);
250     return false;
251 }
252 
253 template <typename T, StorageOrder Order, typename I, typename A>
254 bool
255 FullStorage<T, Order, I, A>::resize(const Range<IndexType> &rows,
256                                     const Range<IndexType> &cols,
257                                     const ElementType &value)
258 {
259     if ((_numRows!=rows.length()) || (_numCols!=cols.length())) {
260         _release();
261         _numRows = rows.length();
262         _numCols = cols.length();
263         _firstRow = rows.firstIndex();
264         _firstCol = cols.firstIndex();
265         _allocate(value);
266         return true;
267     }
268     changeIndexBase(rows.firstIndex(), cols.firstIndex());
269     return false;
270 }
271 
272 template <typename T, StorageOrder Order, typename I, typename A>
273 template <typename FS>
274 bool
275 FullStorage<T, Order, I, A>::resize(const FS &rhs, const ElementType &value)
276 {
277     return resize(rhs.numRows(), rhs.numCols(),
278                   rhs.firstRow(), rhs.firstCol(),
279                   value);
280 }
281 
282 template <typename T, StorageOrder Order, typename I, typename A>
283 bool
284 FullStorage<T, Order, I, A>::fill(const ElementType &value)
285 {
286     ASSERT(_data);
287     std::fill_n(data(), numRows()*numCols(), value);
288     return true;
289 }
290 
291 template <typename T, StorageOrder Order, typename I, typename A>
292 bool
293 FullStorage<T, Order, I, A>::fill(StorageUpLo  upLo,
294                                   const ElementType &value)
295 {
296     ASSERT(_data);
297 
298     trapezoidalFill(order, upLo, value,
299                     numRows(), numCols(),
300                     data(), leadingDimension());
301     return true;
302 }
303 
304 
305 template <typename T, StorageOrder Order, typename I, typename A>
306 void
307 FullStorage<T, Order, I, A>::changeIndexBase(IndexType firstRow,
308                                              IndexType firstCol)
309 {
310     if (Order==RowMajor) {
311         _data = data() - (firstRow*leadingDimension() + firstCol);
312     }
313     if (Order==ColMajor) {
314         _data = data() - (firstCol*leadingDimension() + firstRow);
315     }
316     _firstRow = firstRow;
317     _firstCol = firstCol;
318 }
319 
320 // view of fullstorage scheme as an array
321 template <typename T, StorageOrder Order, typename I, typename A>
322 const typename FullStorage<T, Order, I, A>::ConstArrayView
323 FullStorage<T, Order, I, A>::arrayView(IndexType firstViewIndex) const
324 {
325 #   ifndef NDEBUG
326     if (order==ColMajor) {
327         ASSERT(numRows()==leadingDimension());
328     } else {
329         ASSERT(numCols()==leadingDimension());
330     }
331 #   endif
332 
333     return ConstArrayView(numCols()*numRows(),
334                           data(),
335                           IndexType(1),
336                           firstViewIndex,
337                           allocator());
338 }
339 
340 template <typename T, StorageOrder Order, typename I, typename A>
341 typename FullStorage<T, Order, I, A>::ArrayView
342 FullStorage<T, Order, I, A>::arrayView(IndexType firstViewIndex)
343 {
344 #   ifndef NDEBUG
345     if (order==ColMajor) {
346         ASSERT(numRows()==leadingDimension());
347     } else {
348         ASSERT(numCols()==leadingDimension());
349     }
350 #   endif
351 
352     return ArrayView(numCols()*numRows(),
353                      data(),
354                      IndexType(1),
355                      firstViewIndex,
356                      allocator());
357 }
358 
359 // view of rectangular part
360 template <typename T, StorageOrder Order, typename I, typename A>
361 const typename FullStorage<T, Order, I, A>::ConstView
362 FullStorage<T, Order, I, A>::view(IndexType fromRow, IndexType fromCol,
363                                   IndexType toRow, IndexType toCol,
364                                   IndexType firstViewRow,
365                                   IndexType firstViewCol) const
366 {
367     const IndexType numRows = toRow-fromRow+1;
368     const IndexType numCols = toCol-fromCol+1;
369 
370 #   ifndef NDEBUG
371     // prevent an out-of-bound assertion in case a view is empty anyway
372     if ((numRows==0) || (numCols==0)) {
373         return ConstView(numRows, numCols, 0, leadingDimension(),
374                          firstViewRow, firstViewCol, allocator());
375     }
376 #   endif
377 
378     ASSERT(fromRow>=firstRow());
379     ASSERT(fromRow<=toRow);
380     ASSERT(toRow<=lastRow());
381 
382     ASSERT(fromCol>=firstCol());
383     ASSERT(fromCol<=toCol);
384     ASSERT(toCol<=lastCol());
385 
386     return ConstView(numRows,                         // # rows
387                      numCols,                         // # cols
388                      &(operator()(fromRow, fromCol)), // data
389                      leadingDimension(),              // leading dimension
390                      firstViewRow,                    // firstRow
391                      firstViewCol,                    // firstCol
392                      allocator());                    // allocator
393 }
394 
395 template <typename T, StorageOrder Order, typename I, typename A>
396 typename FullStorage<T, Order, I, A>::View
397 FullStorage<T, Order, I, A>::view(IndexType fromRow, IndexType fromCol,
398                                   IndexType toRow, IndexType toCol,
399                                   IndexType firstViewRow,
400                                   IndexType firstViewCol)
401 {
402     const IndexType numRows = toRow-fromRow+1;
403     const IndexType numCols = toCol-fromCol+1;
404 
405 #   ifndef NDEBUG
406     // prevent an out-of-bound assertion in case a view is empty anyway
407     if ((numRows==0) || (numCols==0)) {
408         return View(numRows, numCols, 0, leadingDimension(),
409                     firstViewRow, firstViewCol, allocator());
410     }
411 #   endif
412 
413     ASSERT(fromRow>=firstRow());
414     ASSERT(fromRow<=toRow);
415     ASSERT(toRow<=lastRow());
416 
417     ASSERT(fromCol>=firstCol());
418     ASSERT(fromCol<=toCol);
419     ASSERT(toCol<=lastCol());
420 
421     return View(numRows,                         // # rows
422                 numCols,                         // # cols
423                 &(operator()(fromRow, fromCol)), // data
424                 leadingDimension(),              // leading dimension
425                 firstViewRow,                    // firstRow
426                 firstViewCol,                    // firstCol
427                 allocator());                    // allocator
428 }
429 
430 // view of single row
431 template <typename T, StorageOrder Order, typename I, typename A>
432 const typename FullStorage<T, Order, I, A>::ConstArrayView
433 FullStorage<T, Order, I, A>::viewRow(IndexType row,
434                                      IndexType firstViewIndex) const
435 {
436 #   ifndef NDEBUG
437     // prevent an out-of-bound assertion in case a view is empty anyway
438     if (numCols()==0) {
439         return ConstArrayView(numCols(), 0, strideCol(),
440                               firstViewIndex, allocator());
441     }
442 #   endif
443 
444     ASSERT(row>=firstRow());
445     ASSERT(row<=lastRow());
446 
447     return ConstArrayView(numCols(),
448                           &(operator()(row, _firstCol)),
449                           strideCol(),
450                           firstViewIndex,
451                           allocator());
452 }
453 
454 template <typename T, StorageOrder Order, typename I, typename A>
455 typename FullStorage<T, Order, I, A>::ArrayView
456 FullStorage<T, Order, I, A>::viewRow(IndexType row,
457                                      IndexType firstViewIndex)
458 {
459 #   ifndef NDEBUG
460     // prevent an out-of-bound assertion in case a view is empty anyway
461     if (numCols()==0) {
462         return ArrayView(numCols(), 0, strideCol(),
463                          firstViewIndex, allocator());
464     }
465 #   endif
466 
467     ASSERT(row>=firstRow());
468     ASSERT(row<=lastRow());
469 
470     return ArrayView(numCols(),
471                      &(this->operator()(row, _firstCol)),
472                      strideCol(),
473                      firstViewIndex,
474                      allocator());
475 }
476 
477 template <typename T, StorageOrder Order, typename I, typename A>
478 const typename FullStorage<T, Order, I, A>::ConstArrayView
479 FullStorage<T, Order, I, A>::viewRow(IndexType row,
480                                      IndexType firstCol, IndexType lastCol,
481                                      IndexType firstViewIndex) const
482 {
483     const IndexType length = lastCol-firstCol+1;
484 
485 #   ifndef NDEBUG
486     // prevent an out-of-bound assertion in case a view is empty anyway
487     if (length==0) {
488         return ConstArrayView(length, 0, strideCol(),
489                               firstViewIndex, allocator());
490     }
491 #   endif
492 
493     ASSERT(row>=firstRow());
494     ASSERT(row<=lastRow());
495 
496     return ConstArrayView(length,
497                           &(operator()(row, firstCol)),
498                           strideCol(),
499                           firstViewIndex,
500                           allocator());
501 }
502 
503 template <typename T, StorageOrder Order, typename I, typename A>
504 typename FullStorage<T, Order, I, A>::ArrayView
505 FullStorage<T, Order, I, A>::viewRow(IndexType row,
506                                      IndexType firstCol, IndexType lastCol,
507                                      IndexType firstViewIndex)
508 {
509     const IndexType length = lastCol-firstCol+1;
510 
511 #   ifndef NDEBUG
512     // prevent an out-of-bound assertion in case a view is empty anyway
513     if (length==0) {
514         return ArrayView(length, 0, strideCol(),
515                          firstViewIndex, allocator());
516     }
517 #   endif
518 
519     ASSERT(row>=firstRow());
520     ASSERT(row<=lastRow());
521 
522     return ArrayView(length,
523                      &(operator()(row, firstCol)),
524                      strideCol(),
525                      firstViewIndex,
526                      allocator());
527 }
528 
529 // view of single column
530 template <typename T, StorageOrder Order, typename I, typename A>
531 const typename FullStorage<T, Order, I, A>::ConstArrayView
532 FullStorage<T, Order, I, A>::viewCol(IndexType col,
533                                      IndexType firstViewIndex) const
534 {
535 #   ifndef NDEBUG
536     // prevent an out-of-bound assertion in case a view is empty anyway
537     if (numRows()==0) {
538         return ConstArrayView(numRows(), 0, strideRow(),
539                               firstViewIndex, allocator());
540     }
541 #   endif
542 
543     ASSERT(col>=firstCol());
544     ASSERT(col<=lastCol());
545 
546     return ConstArrayView(numRows(),
547                           &(this->operator()(_firstRow, col)),
548                           strideRow(),
549                           firstViewIndex,
550                           allocator());
551 }
552 
553 template <typename T, StorageOrder Order, typename I, typename A>
554 typename FullStorage<T, Order, I, A>::ArrayView
555 FullStorage<T, Order, I, A>::viewCol(IndexType col,
556                                      IndexType firstViewIndex)
557 {
558 #   ifndef NDEBUG
559     // prevent an out-of-bound assertion in case a view is empty anyway
560     if (numRows()==0) {
561         return ArrayView(numRows(), 0, strideRow(),
562                          firstViewIndex, allocator());
563     }
564 #   endif
565 
566     ASSERT(col>=firstCol());
567     ASSERT(col<=lastCol());
568 
569     return ArrayView(numRows(),
570                      &(operator()(_firstRow, col)),
571                      strideRow(),
572                      firstViewIndex,
573                      allocator());
574 }
575 
576 template <typename T, StorageOrder Order, typename I, typename A>
577 const typename FullStorage<T, Order, I, A>::ConstArrayView
578 FullStorage<T, Order, I, A>::viewCol(IndexType firstRow, IndexType lastRow,
579                                      IndexType col,
580                                      IndexType firstViewIndex) const
581 {
582     const IndexType length = lastRow-firstRow+1;
583 
584 #   ifndef NDEBUG
585     // prevent an out-of-bound assertion in case a view is empty anyway
586     if (length==0) {
587         return ConstArrayView(length, 0, strideRow(),
588                               firstViewIndex, allocator());
589     }
590 #   endif
591 
592     ASSERT(col>=firstCol());
593     ASSERT(col<=lastCol());
594 
595     return ConstArrayView(length,
596                           &(this->operator()(firstRow, col)),
597                           strideRow(),
598                           firstViewIndex,
599                           allocator());
600 }
601 
602 template <typename T, StorageOrder Order, typename I, typename A>
603 typename FullStorage<T, Order, I, A>::ArrayView
604 FullStorage<T, Order, I, A>::viewCol(IndexType firstRow, IndexType lastRow,
605                                      IndexType col,
606                                      IndexType firstViewIndex)
607 {
608     const IndexType length = lastRow-firstRow+1;
609 
610 #   ifndef NDEBUG
611     // prevent an out-of-bound assertion in case a view is empty anyway
612     if (length==0) {
613         return ArrayView(length, 0, strideRow(),
614                          firstViewIndex, allocator());
615     }
616 #   endif
617 
618     ASSERT(col>=firstCol());
619     ASSERT(col<=lastCol());
620 
621     return ArrayView(length,
622                      &(this->operator()(firstRow, col)),
623                      strideRow(),
624                      firstViewIndex,
625                      allocator());
626 }
627 
628 // view of d-th diagonal
629 template <typename T, StorageOrder Order, typename I, typename A>
630 const typename FullStorage<T, Order, I, A>::ConstArrayView
631 FullStorage<T, Order, I, A>::viewDiag(IndexType d,
632                                       IndexType firstViewIndex) const
633 {
634     IndexType _row = (d>0) ? 0 : -d;
635     IndexType _col = (d>0) ? d :  0;
636 
637     IndexType row = firstRow() + _row;
638     IndexType col = firstCol() + _col;
639 
640     return ConstArrayView(std::min(numRows()-_row, numCols()-_col),
641                           &(this->operator()(row,col)),
642                           leadingDimension()+1,
643                           firstViewIndex,
644                           allocator());
645 }
646 
647 template <typename T, StorageOrder Order, typename I, typename A>
648 typename FullStorage<T, Order, I, A>::ArrayView
649 FullStorage<T, Order, I, A>::viewDiag(IndexType d,
650                                       IndexType firstViewIndex)
651 {
652     IndexType _row = (d>0) ? 0 : -d;
653     IndexType _col = (d>0) ? d :  0;
654 
655     IndexType row = firstRow() + _row;
656     IndexType col = firstCol() + _col;
657 
658     return ArrayView(std::min(numRows()-_row, numCols()-_col),
659                      &(this->operator()(row,col)),
660                      leadingDimension()+1,
661                      firstViewIndex,
662                      allocator());
663 }
664 
665 //-- Private Methods -----------------------------------------------------------
666 
667 template <typename T, StorageOrder Order, typename I, typename A>
668 void
669 FullStorage<T, Order, I, A>::_setIndexBase(IndexType firstRow,
670                                            IndexType firstCol)
671 {
672     // assume: _data points to allocated memory
673     
674     if (Order==RowMajor) {
675         _data -= firstRow*_numCols + firstCol;
676     }
677     if (Order==ColMajor) {
678         _data -= firstCol*_numRows + firstRow;
679     }
680     _firstRow = firstRow;
681     _firstCol = firstCol;
682 }
683 
684 template <typename T, StorageOrder Order, typename I, typename A>
685 void
686 FullStorage<T, Order, I, A>::_raw_allocate()
687 {
688     ASSERT(!_data);
689     ASSERT(_numRows>0);
690     ASSERT(_numCols>0);
691 
692     _data = _allocator.allocate(_numRows*_numCols);
693 #ifndef NDEBUG
694     ElementType *p = _data;
695 #endif
696     _setIndexBase(_firstRow, _firstCol);
697     ASSERT(_data);
698 #ifndef NDEBUG
699     ASSERT(p==data());
700 #endif
701 }
702 
703 template <typename T, StorageOrder Order, typename I, typename A>
704 void
705 FullStorage<T, Order, I, A>::_allocate(const ElementType &value)
706 {
707     const IndexType numElements = numRows()*numCols();
708     
709     if (numElements==0) {
710         return;
711     }
712 
713     _raw_allocate();
714     T *p = data();
715     for (IndexType i=0; i<numElements; ++i) {
716         _allocator.construct(p++, value);
717     }
718 }
719 
720 template <typename T, StorageOrder Order, typename I, typename A>
721 void
722 FullStorage<T, Order, I, A>::_release()
723 {
724     if (_data) {
725         T *p = data();
726         for (IndexType i=0; i<numRows()*numCols(); ++i) {
727             _allocator.destroy(p++);
728         }
729         _allocator.deallocate(data(), numRows()*numCols());
730         _data = 0;
731     }
732     ASSERT(_data==0);
733 }
734 
735 // namespace flens
736 
737 #endif // FLENS_STORAGE_FULLSTORAGE_FULLSTORAGE_TCC