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_MATRIXTYPES_GENERAL_IMPL_GEMATRIX_TCC
 34 #define FLENS_MATRIXTYPES_GENERAL_IMPL_GEMATRIX_TCC 1
 35 
 36 #include <flens/blas/blas.h>
 37 #include <flens/typedefs.h>
 38 
 39 #include <flens/matrixtypes/general/impl/ge/constelementclosure.tcc>
 40 #include <flens/matrixtypes/general/impl/ge/elementclosure.tcc>
 41 #include <flens/matrixtypes/general/impl/ge/initializer.tcc>
 42 
 43 namespace flens {
 44 
 45 template <typename FS>
 46 GeMatrix<FS>::GeMatrix()
 47 {
 48 }
 49 
 50 template <typename FS>
 51 GeMatrix<FS>::GeMatrix(IndexType numRows, IndexType numCols)
 52     : _engine(numRows, numCols)
 53 {
 54     ASSERT(numRows>=0);
 55     ASSERT(numCols>=0);
 56 }
 57 
 58 template <typename FS>
 59 GeMatrix<FS>::GeMatrix(IndexType numRows, IndexType numCols,
 60                        IndexType firstRow, IndexType firstCol)
 61     : _engine(numRows, numCols, firstRow, firstCol)
 62 {
 63     ASSERT(numRows>=0);
 64     ASSERT(numCols>=0);
 65 }
 66 
 67 template <typename FS>
 68 GeMatrix<FS>::GeMatrix(const Range<IndexType> &rowRange,
 69                        const Range<IndexType> &colRange)
 70     : _engine(rowRange.numTicks(), colRange.numTicks(),
 71               rowRange.firstIndex(), colRange.firstIndex())
 72 {
 73 }
 74 
 75 template <typename FS>
 76 GeMatrix<FS>::GeMatrix(const Engine &engine)
 77     : _engine(engine)
 78 {
 79 }
 80 
 81 template <typename FS>
 82 GeMatrix<FS>::GeMatrix(const GeMatrix &rhs)
 83     : GeneralMatrix<GeMatrix>(),
 84       _engine(rhs._engine)
 85 {
 86 }
 87 
 88 template <typename FS>
 89 template <typename RHS>
 90 GeMatrix<FS>::GeMatrix(const GeMatrix<RHS> &rhs)
 91     : _engine(rhs.engine())
 92 {
 93 }
 94 
 95 template <typename FS>
 96 template <typename RHS>
 97 GeMatrix<FS>::GeMatrix(GeMatrix<RHS> &rhs)
 98     : _engine(rhs.engine())
 99 {
100 }
101 
102 template <typename FS>
103 template <typename RHS>
104 GeMatrix<FS>::GeMatrix(const Matrix<RHS> &rhs)
105 {
106     assign(rhs, *this);
107 }
108 
109 template <typename FS>
110 template <typename VECTOR>
111 GeMatrix<FS>::GeMatrix(IndexType numRows, IndexType numCols, VECTOR &&rhs)
112     : _engine(numRows, numCols, rhs.engine(), (FS::order==RowMajor) ? numCols
113                                                                     : numRows)
114 {
115 }
116 
117 template <typename FS>
118 template <typename VECTOR>
119 GeMatrix<FS>::GeMatrix(IndexType numRows, IndexType numCols,
120                        VECTOR &&rhs,
121                        IndexType leadingDimension)
122     : _engine(numRows, numCols, rhs.engine(), leadingDimension)
123 {
124 }
125 
126 // -- operators ----------------------------------------------------------------
127 
128 template <typename FS>
129 typename GeMatrix<FS>::Initializer
130 GeMatrix<FS>::operator=(const ElementType &value)
131 {
132     engine().fill(value);
133     return Initializer(*this, firstRow(), firstCol());
134 }
135 
136 template <typename FS>
137 GeMatrix<FS> &
138 GeMatrix<FS>::operator=(const GeMatrix<FS> &rhs)
139 {
140     if (this!=&rhs) {
141         assign(rhs, *this);
142     }
143     return *this;
144 }
145 
146 template <typename FS>
147 template <typename RHS>
148 GeMatrix<FS> &
149 GeMatrix<FS>::operator=(const Matrix<RHS> &rhs)
150 {
151     assign(rhs, *this);
152     return *this;
153 }
154 
155 template <typename FS>
156 template <typename RHS>
157 GeMatrix<FS> &
158 GeMatrix<FS>::operator+=(const Matrix<RHS> &rhs)
159 {
160     plusAssign(rhs, *this);
161     return *this;
162 }
163 
164 template <typename FS>
165 template <typename RHS>
166 GeMatrix<FS> &
167 GeMatrix<FS>::operator-=(const Matrix<RHS> &rhs)
168 {
169     minusAssign(rhs, *this);
170     return *this;
171 }
172 
173 template <typename FS>
174 GeMatrix<FS> &
175 GeMatrix<FS>::operator+=(const ElementType &alpha)
176 {
177     const Underscore<IndexType> _;
178 
179     if (order()==ColMajor) {
180         for (IndexType j=firstCol(); j<=lastCol(); ++j) {
181             (*this)(_,j) += alpha;
182         }
183     } else {
184         for (IndexType i=firstRow(); i<=lastRow(); ++i) {
185             (*this)(i,_) += alpha;
186         }
187     }
188 }
189 
190 template <typename FS>
191 GeMatrix<FS> &
192 GeMatrix<FS>::operator-=(const ElementType &alpha)
193 {
194     const Underscore<IndexType> _;
195 
196     if (order()==ColMajor) {
197         for (IndexType j=firstCol(); j<=lastCol(); ++j) {
198             (*this)(_,j) -= alpha;
199         }
200     } else {
201         for (IndexType i=firstRow(); i<=lastRow(); ++i) {
202             (*this)(i,_) -= alpha;
203         }
204     }
205 }
206 
207 template <typename FS>
208 GeMatrix<FS> &
209 GeMatrix<FS>::operator*=(const ElementType &alpha)
210 {
211     blas::scal(alpha, *this);
212     return *this;
213 }
214 
215 template <typename FS>
216 GeMatrix<FS> &
217 GeMatrix<FS>::operator/=(const ElementType &alpha)
218 {
219     blas::rscal(alpha, *this);
220     return *this;
221 }
222 
223 template <typename FS>
224 inline
225 const typename GeMatrix<FS>::ElementType &
226 GeMatrix<FS>::operator()(IndexType row, IndexType col) const
227 {
228     return _engine(row, col);
229 }
230 
231 template <typename FS>
232 inline
233 typename GeMatrix<FS>::ElementType &
234 GeMatrix<FS>::operator()(IndexType row, IndexType col)
235 {
236     return _engine(row, col);
237 }
238 
239 template <typename FS>
240 template <typename S>
241 const gematrix::ConstElementClosure<GeMatrix<FS>, typename Scalar<S>::Impl>
242 GeMatrix<FS>::operator()(const Scalar<S> &row, const Scalar<S> &col) const
243 {
244     typedef typename Scalar<S>::Impl ScalarImpl;
245     typedef gematrix::ConstElementClosure<GeMatrix, ScalarImpl>  CEC;
246     return CEC(*this, row.impl(), col.impl());
247 
248 }
249 
250 template <typename FS>
251 const typename GeMatrix<FS>::ConstElementClosure
252 GeMatrix<FS>::operator()(const IndexVariable &row,
253                          const IndexVariable &col) const
254 {
255     return ConstElementClosure(*this, row, col);
256 }
257 
258 template <typename FS>
259 typename GeMatrix<FS>::ElementClosure
260 GeMatrix<FS>::operator()(IndexVariable &row, IndexVariable &col)
261 {
262     return ElementClosure(*this, row, col);
263 }
264 
265 // -- methods ------------------------------------------------------------------
266 
267 template <typename FS>
268 typename GeMatrix<FS>::IndexType
269 GeMatrix<FS>::numRows() const
270 {
271     return _engine.numRows();
272 }
273 
274 template <typename FS>
275 typename GeMatrix<FS>::IndexType
276 GeMatrix<FS>::numCols() const
277 {
278     return _engine.numCols();
279 }
280 
281 template <typename FS>
282 typename GeMatrix<FS>::IndexType
283 GeMatrix<FS>::firstRow() const
284 {
285     return _engine.firstRow();
286 }
287 
288 template <typename FS>
289 typename GeMatrix<FS>::IndexType
290 GeMatrix<FS>::lastRow() const
291 {
292     return _engine.lastRow();
293 }
294 
295 template <typename FS>
296 typename GeMatrix<FS>::IndexType
297 GeMatrix<FS>::firstCol() const
298 {
299     return _engine.firstCol();
300 }
301 
302 template <typename FS>
303 typename GeMatrix<FS>::IndexType
304 GeMatrix<FS>::lastCol() const
305 {
306     return _engine.lastCol();
307 }
308 
309 template <typename FS>
310 Range<typename GeMatrix<FS>::IndexType>
311 GeMatrix<FS>::rows() const
312 {
313     return Range<IndexType>(_engine.firstRow(),_engine.lastRow());
314 }
315 
316 template <typename FS>
317 Range<typename GeMatrix<FS>::IndexType>
318 GeMatrix<FS>::cols() const
319 {
320     return Range<IndexType>(_engine.firstCol(),_engine.lastCol());
321 }
322 
323 template <typename FS>
324 const typename GeMatrix<FS>::ElementType *
325 GeMatrix<FS>::data() const
326 {
327     return _engine.data();
328 }
329 
330 template <typename FS>
331 typename GeMatrix<FS>::ElementType *
332 GeMatrix<FS>::data()
333 {
334     return _engine.data();
335 }
336 
337 template <typename FS>
338 typename GeMatrix<FS>::IndexType
339 GeMatrix<FS>::leadingDimension() const
340 {
341     return _engine.leadingDimension();
342 }
343 
344 template <typename FS>
345 StorageOrder
346 GeMatrix<FS>::order() const
347 {
348     return _engine.order;
349 }
350 
351 template <typename FS>
352 template <typename RHS>
353 bool
354 GeMatrix<FS>::resize(const GeMatrix<RHS> &rhs,
355                      const ElementType &value)
356 {
357     return _engine.resize(rhs.engine(), value);
358 }
359 
360 template <typename FS>
361 bool
362 GeMatrix<FS>::resize(IndexType numRows, IndexType numCols,
363                      IndexType firstRowIndex,
364                      IndexType firstColIndex,
365                      const ElementType &value)
366 {
367     return _engine.resize(numRows, numCols,
368                           firstRowIndex, firstColIndex,
369                           value);
370 }
371 
372 template <typename FS>
373 bool
374 GeMatrix<FS>::fill(const ElementType &value)
375 {
376     return _engine.fill(value);
377 }
378 
379 template <typename FS>
380 void
381 GeMatrix<FS>::changeIndexBase(IndexType firstRowIndex, IndexType firstColIndex)
382 {
383     _engine.changeIndexBase(firstRowIndex, firstColIndex);
384 }
385 
386 
387 // -- views --------------------------------------------------------------------
388 // vectorize matrix
389 template <typename FS>
390 const typename GeMatrix<FS>::ConstVectorView
391 GeMatrix<FS>::vectorView() const
392 {
393     return _engine.arrayView();
394 }
395 
396 template <typename FS>
397 typename GeMatrix<FS>::VectorView
398 GeMatrix<FS>::vectorView()
399 {
400     return _engine.arrayView();
401 }
402 
403 // vectorize matrix and select range
404 template <typename FS>
405 const typename GeMatrix<FS>::ConstVectorView
406 GeMatrix<FS>::vectorView(IndexType from, IndexType to) const
407 {
408     return _engine.arrayView().view(from,to);
409 }
410 
411 template <typename FS>
412 typename GeMatrix<FS>::VectorView
413 GeMatrix<FS>::vectorView(IndexType from, IndexType to)
414 {
415     return _engine.arrayView().view(from,to);
416 }
417 
418 // diag views
419 template <typename FS>
420 const typename GeMatrix<FS>::ConstVectorView
421 GeMatrix<FS>::diag(IndexType d) const
422 {
423     return _engine.viewDiag(d);
424 }
425 
426 template <typename FS>
427 typename GeMatrix<FS>::VectorView
428 GeMatrix<FS>::diag(IndexType d)
429 {
430     return _engine.viewDiag(d);
431 }
432 
433 // triangular views
434 template <typename FS>
435 const typename GeMatrix<FS>::ConstTriangularView
436 GeMatrix<FS>::upper() const
437 {
438     return ConstTriangularView(engine(), Upper);
439 }
440 
441 template <typename FS>
442 typename GeMatrix<FS>::TriangularView
443 GeMatrix<FS>::upper()
444 {
445     return TriangularView(engine(), Upper);
446 }
447 
448 template <typename FS>
449 const typename GeMatrix<FS>::ConstTriangularView
450 GeMatrix<FS>::upperUnit() const
451 {
452     return ConstTriangularView(engine(), Upper, Unit);
453 }
454 
455 template <typename FS>
456 typename GeMatrix<FS>::TriangularView
457 GeMatrix<FS>::upperUnit()
458 {
459     return TriangularView(engine(), Upper, Unit);
460 }
461 
462 template <typename FS>
463 const typename GeMatrix<FS>::ConstTriangularView
464 GeMatrix<FS>::strictUpper() const
465 {
466     const Underscore<IndexType> _;
467     const IndexType n = numCols();
468 
469     return operator()(_,_(2,n)).upper();
470 }
471 
472 template <typename FS>
473 typename GeMatrix<FS>::TriangularView
474 GeMatrix<FS>::strictUpper()
475 {
476     const Underscore<IndexType> _;
477     const IndexType n = numCols();
478 
479     return operator()(_,_(2,n)).upper();
480 }
481 
482 template <typename FS>
483 const typename GeMatrix<FS>::ConstTriangularView
484 GeMatrix<FS>::lower() const
485 {
486     return ConstTriangularView(engine(), Lower);
487 }
488 
489 template <typename FS>
490 typename GeMatrix<FS>::TriangularView
491 GeMatrix<FS>::lower()
492 {
493     return TriangularView(engine(), Lower);
494 }
495 
496 template <typename FS>
497 const typename GeMatrix<FS>::ConstTriangularView
498 GeMatrix<FS>::lowerUnit() const
499 {
500     return ConstTriangularView(engine(), Lower, Unit);
501 }
502 
503 template <typename FS>
504 typename GeMatrix<FS>::TriangularView
505 GeMatrix<FS>::lowerUnit()
506 {
507     return TriangularView(engine(), Lower, Unit);
508 }
509 
510 template <typename FS>
511 const typename GeMatrix<FS>::ConstTriangularView
512 GeMatrix<FS>::strictLower() const
513 {
514     const Underscore<IndexType> _;
515     const IndexType m = numRows();
516 
517     return operator()(_(2,m),_).lower();
518 }
519 
520 template <typename FS>
521 typename GeMatrix<FS>::TriangularView
522 GeMatrix<FS>::strictLower()
523 {
524     const Underscore<IndexType> _;
525     const IndexType m = numRows();
526 
527     return operator()(_(2,m),_).lower();
528 }
529 
530 // rectangular views
531 template <typename FS>
532 const typename GeMatrix<FS>::ConstView
533 GeMatrix<FS>::operator()(const Range<IndexType> &rows,
534                          const Range<IndexType> &cols) const
535 {
536     ASSERT(rows.stride()==IndexType(1));
537     ASSERT(cols.stride()==IndexType(1));
538     return engine().view(rows.firstIndex(), cols.firstIndex(),
539                          rows.lastIndex(), cols.lastIndex());
540 }
541 
542 template <typename FS>
543 typename GeMatrix<FS>::View
544 GeMatrix<FS>::operator()(const Range<IndexType> &rows,
545                          const Range<IndexType> &cols)
546 {
547     ASSERT(rows.stride()==IndexType(1));
548     ASSERT(cols.stride()==IndexType(1));
549     return engine().view(rows.firstIndex(), cols.firstIndex(),
550                          rows.lastIndex(), cols.lastIndex());
551 }
552 
553 template <typename FS>
554 template <typename RHS>
555 const typename GeMatrix<FS>::ConstView
556 GeMatrix<FS>::operator()(const GeMatrix<RHS> &A) const
557 {
558     return engine().view(A.firstRow(), A.firstCol(),
559                          A.lastRow(), A.lastCol());
560 }
561 
562 template <typename FS>
563 template <typename RHS>
564 typename GeMatrix<FS>::View
565 GeMatrix<FS>::operator()(const GeMatrix<RHS> &A)
566 {
567     return engine().view(A.firstRow(), A.firstCol(),
568                          A.lastRow(), A.lastCol());
569 }
570 
571 // rectangular views (all rows selected)
572 template <typename FS>
573 const typename GeMatrix<FS>::ConstView
574 GeMatrix<FS>::operator()(const Underscore<IndexType> &,
575                          const Range<IndexType> &cols) const
576 {
577     ASSERT(cols.stride()==IndexType(1));
578     return engine().view(firstRow(), cols.firstIndex(),
579                          lastRow(), cols.lastIndex());
580 }
581 
582 template <typename FS>
583 typename GeMatrix<FS>::View
584 GeMatrix<FS>::operator()(const Underscore<IndexType> &,
585                          const Range<IndexType> &cols)
586 {
587     ASSERT(cols.stride()==IndexType(1));
588     return engine().view(firstRow(), cols.firstIndex(),
589                          lastRow(), cols.lastIndex());
590 }
591 
592 // rectangular views (all columns selected)
593 template <typename FS>
594 const typename GeMatrix<FS>::ConstView
595 GeMatrix<FS>::operator()(const Range<IndexType> &rows,
596                          const Underscore<IndexType> &) const
597 {
598     ASSERT(rows.stride()==IndexType(1));
599     return engine().view(rows.firstIndex(), firstCol(),
600                          rows.lastIndex(), lastCol());
601 }
602 
603 template <typename FS>
604 typename GeMatrix<FS>::View
605 GeMatrix<FS>::operator()(const Range<IndexType> &rows,
606                          const Underscore<IndexType> &)
607 {
608     ASSERT(rows.stride()==IndexType(1));
609     return engine().view(rows.firstIndex(), firstCol(),
610                          rows.lastIndex(), lastCol());
611 }
612 
613 // row view (vector view)
614 template <typename FS>
615 const typename GeMatrix<FS>::ConstVectorView
616 GeMatrix<FS>::operator()(IndexType row, const Underscore<IndexType> &) const
617 {
618     return engine().viewRow(row);
619 }
620 
621 template <typename FS>
622 typename GeMatrix<FS>::VectorView
623 GeMatrix<FS>::operator()(IndexType row, const Underscore<IndexType> &)
624 {
625     return engine().viewRow(row);
626 }
627 
628 template <typename FS>
629 const typename GeMatrix<FS>::ConstVectorView
630 GeMatrix<FS>::operator()(IndexType row, const Range<IndexType> &cols) const
631 {
632     return engine().viewRow(row, cols.firstIndex(), cols.lastIndex());
633 }
634 
635 template <typename FS>
636 typename GeMatrix<FS>::VectorView
637 GeMatrix<FS>::operator()(IndexType row, const Range<IndexType> &cols)
638 {
639     return engine().viewRow(row, cols.firstIndex(), cols.lastIndex());
640 }
641 
642 // column view (vector view)
643 template <typename FS>
644 const typename GeMatrix<FS>::ConstVectorView
645 GeMatrix<FS>::operator()(const Underscore<IndexType> &, IndexType col) const
646 {
647     return engine().viewCol(col);
648 }
649 
650 template <typename FS>
651 typename GeMatrix<FS>::VectorView
652 GeMatrix<FS>::operator()(const Underscore<IndexType> &, IndexType col)
653 {
654     return engine().viewCol(col);
655 }
656 
657 template <typename FS>
658 const typename GeMatrix<FS>::ConstVectorView
659 GeMatrix<FS>::operator()(const Range<IndexType> &rows, IndexType col) const
660 {
661     return engine().viewCol(rows.firstIndex(), rows.lastIndex(), col);
662 }
663 
664 template <typename FS>
665 typename GeMatrix<FS>::VectorView
666 GeMatrix<FS>::operator()(const Range<IndexType> &rows, IndexType col)
667 {
668     return engine().viewCol(rows.firstIndex(), rows.lastIndex(), col);
669 }
670 
671 // -- implementation -----------------------------------------------------------
672 
673 template <typename FS>
674 const typename GeMatrix<FS>::Engine &
675 GeMatrix<FS>::engine() const
676 {
677     return _engine;
678 }
679 
680 template <typename FS>
681 typename GeMatrix<FS>::Engine &
682 GeMatrix<FS>::engine()
683 {
684     return _engine;
685 }
686 
687 // namespace flens
688 
689 #endif // FLENS_MATRIXTYPES_GENERAL_IMPL_GEMATRIX_TCC