matrix_assign.hpp
Go to the documentation of this file.
1 /*!
2  * \brief Kernels for matrix-expression assignments
3  *
4  * \author O. Krause
5  * \date 2013
6  *
7  *
8  * \par Copyright 1995-2015 Shark Development Team
9  *
10  * <BR><HR>
11  * This file is part of Shark.
12  * <http://image.diku.dk/shark/>
13  *
14  * Shark is free software: you can redistribute it and/or modify
15  * it under the terms of the GNU Lesser General Public License as published
16  * by the Free Software Foundation, either version 3 of the License, or
17  * (at your option) any later version.
18  *
19  * Shark is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU Lesser General Public License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public License
25  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26  *
27  */
28 #ifndef SHARK_LINALG_BLAS_KERNELS_MATRIX_ASSIGN_HPP
29 #define SHARK_LINALG_BLAS_KERNELS_MATRIX_ASSIGN_HPP
30 
31 #include "vector_assign.hpp"
32 #include <algorithm>
33 namespace shark {
34 namespace blas {
35 
36 namespace detail{
37 
38 //base version of the matrix_transpose proxy which does not rely on the existence of the assignment functions
39 //we define it here so that we don't need to duplicate all versions of assign for the case of row/column major
40 //first argument
41 template<class M>
42 class internal_transpose_proxy: public matrix_expression<internal_transpose_proxy<M> > {
43 public:
44  typedef typename M::size_type size_type;
45  typedef typename M::difference_type difference_type;
46  typedef typename M::value_type value_type;
47  typedef typename M::scalar_type scalar_type;
48  typedef typename M::const_reference const_reference;
49  typedef typename reference<M>::type reference;
50  typedef typename M::const_pointer const_pointer;
51  typedef typename M::pointer pointer;
52 
53  typedef typename M::index_type index_type;
54  typedef typename M::const_index_pointer const_index_pointer;
55  typedef typename index_pointer<M>::type index_pointer;
56 
57  typedef typename closure<M>::type matrix_closure_type;
58  typedef const internal_transpose_proxy<M> const_closure_type;
59  typedef internal_transpose_proxy<M> closure_type;
60  typedef typename M::orientation::transposed_orientation orientation;
61  typedef typename M::storage_category storage_category;
62  typedef typename M::evaluation_category evaluation_category;
63 
64  // Construction and destruction
65  explicit internal_transpose_proxy(matrix_closure_type m):
66  m_expression(m) {}
67 
68  // Accessors
69  size_type size1() const {
70  return m_expression.size2();
71  }
72  size_type size2() const {
73  return m_expression.size1();
74  }
75 
76  // Expression accessors
77  matrix_closure_type const &expression() const{
78  return m_expression;
79  }
80  matrix_closure_type &expression(){
81  return m_expression;
82  }
83 
84  // Element access
85  reference operator()(size_type i, size_type j)const{
86  return m_expression(j, i);
87  }
88 
89  typedef typename matrix_closure_type::const_column_iterator const_row_iterator;
90  typedef typename matrix_closure_type::column_iterator row_iterator;
91  typedef typename matrix_closure_type::const_row_iterator const_column_iterator;
92  typedef typename matrix_closure_type::row_iterator column_iterator;
93 
94  //iterators
95  const_row_iterator row_begin(std::size_t i) const {
96  return m_expression.column_begin(i);
97  }
98  const_row_iterator row_end(std::size_t i) const {
99  return m_expression.column_end(i);
100  }
101  const_column_iterator column_begin(std::size_t j) const {
102  return m_expression.row_begin(j);
103  }
104  const_column_iterator column_end(std::size_t j) const {
105  return m_expression.row_end(j);
106  }
107 
108  row_iterator row_begin(std::size_t i) {
109  return m_expression.column_begin(i);
110  }
111  row_iterator row_end(std::size_t i) {
112  return m_expression.column_end(i);
113  }
114  column_iterator column_begin(std::size_t j) {
115  return m_expression.row_begin(j);
116  }
117  column_iterator column_end(std::size_t j) {
118  return m_expression.row_end(j);
119  }
120 
121  typedef typename major_iterator<internal_transpose_proxy<M> >::type major_iterator;
122 
123  major_iterator set_element(major_iterator pos, size_type index, value_type value){
124  return m_expression.set_element(pos,index,value);
125  }
126 
127  major_iterator clear_range(major_iterator start, major_iterator end){
128  return m_expression.clear_range(start,end);
129  }
130 
131  major_iterator clear_element(major_iterator elem){
132  return m_expression.clear_element(elem);
133  }
134 
135  void clear(){
136  expression().clear();
137  }
138 
139  void reserve(size_type non_zeros) {
140  m_expression.reserve(non_zeros);
141  }
142 
143  void reserve_row(size_type row, size_type non_zeros) {
144  m_expression.reserve_row(row,non_zeros);
145  }
146  void reserve_column(size_type column, size_type non_zeros) {
147  m_expression.reserve_column(column,non_zeros);
148  }
149 private:
150  matrix_closure_type m_expression;
151 };
152 
153 }
154 
155 //////////////////////////////////////////////////////
156 ////Scalar Assignment to Matrix
157 /////////////////////////////////////////////////////
158 
159 namespace kernels{
160 // Explicitly iterating row major
161 template<template <class T1, class T2> class F, class M>
162 void assign(
164  typename M::value_type t,
165  row_major
166 ){
167  for(std::size_t i = 0; i != m().size1(); ++i){
168  matrix_row<M> rowM(m(),i);
169  assign<F>(rowM,t);
170  }
171 }
172 // Explicitly iterating column major
173 template<template <class T1, class T2> class F, class M>
174 void assign(
176  typename M::value_type t,
177  column_major
178 ){
179  for(std::size_t i = 0; i != m().size2(); ++i){
180  matrix_column<M> columnM(m(),i);
181  assign<F>(columnM,t);
182  }
183 }
184 
185 
186 // Spcial case packed - just calls the first two implementations.
187 template<template <class T1, class T2> class F, class M, class Orientation, class Triangular>
188 void assign(
190  typename M::value_type t,
191  packed<Orientation,Triangular>
192 ){
193  assign<F>(m,t,Orientation());
194 }
195 
196 // Dispatcher
197 template<template <class T1, class T2> class F, class M>
198 void assign(
200  typename M::value_type t
201 ){
202  typedef typename M::orientation orientation;
203  assign<F> (m, t, orientation());
204 }
205 
206 /////////////////////////////////////////////////////////////////
207 //////Matrix Assignment implmenting op=
208 ////////////////////////////////////////////////////////////////
209 
210 //direct assignment without functor
211 //the cases were both arguments have the same orientation can be implemented using assign.
212 template<class M, class E,class TagE, class TagM>
213 void assign(
215  matrix_expression<E> const& e,
216  row_major, row_major,TagE, TagM
217 ) {
218  for(std::size_t i = 0; i != m().size1(); ++i){
219  matrix_row<M> rowM(m(),i);
220  kernels::assign(rowM,row(e,i));
221  }
222 }
223 
224 //remain the versions where both argumnts to not have the same orientation
225 
226 //dense-dense case
227 template<class M, class E>
228 void assign(
230  matrix_expression<E> const& e,
231  row_major, column_major,dense_random_access_iterator_tag, dense_random_access_iterator_tag
232 ) {
233  //compute blockwise and wrelem the transposed block.
234  std::size_t const blockSize = 16;
235  typename M::value_type blockStorage[blockSize][blockSize];
236 
237  typedef typename M::size_type size_type;
238  size_type size1 = m().size1();
239  size_type size2 = m().size2();
240  for (size_type iblock = 0; iblock < size1; iblock += blockSize){
241  for (size_type jblock = 0; jblock < size2; jblock += blockSize){
242  std::size_t blockSizei = std::min(blockSize,size1-iblock);
243  std::size_t blockSizej = std::min(blockSize,size2-jblock);
244 
245  //compute block values
246  for (size_type j = 0; j < blockSizej; ++j){
247  for (size_type i = 0; i < blockSizei; ++i){
248  blockStorage[i][j] = e()(iblock+i,jblock+j);
249  }
250  }
251 
252  //copy block in a different order to m
253  for (size_type i = 0; i < blockSizei; ++i){
254  for (size_type j = 0; j < blockSizej; ++j){
255  m()(iblock+i,jblock+j) = blockStorage[i][j];
256  }
257  }
258  }
259  }
260 }
261 
262 // dense-sparse
263 template<class M, class E>
264 void assign(
266  matrix_expression<E> const& e,
267  row_major, column_major,dense_random_access_iterator_tag, sparse_bidirectional_iterator_tag
268 ) {
269  for(std::size_t i = 0; i != m().size2(); ++i){
270  matrix_column<M> columnM(m(),i);
271  kernels::assign(columnM,column(e,i));
272  }
273 }
274 
275 
276 //sparse-dense
277 template<class M, class E>
278 void assign(
280  matrix_expression<E> const& e,
281  row_major, column_major, sparse_bidirectional_iterator_tag, dense_random_access_iterator_tag
282 ) {
283  for(std::size_t i = 0; i != m().size1(); ++i){
284  matrix_column<M> rowM(m(),i);
285  kernels::assign(rowM,row(e,i));
286  }
287 }
288 
289 //sparse-sparse
290 template<class M, class E>
291 void assign(
293  matrix_expression<E> const& e,
294  row_major, column_major,sparse_bidirectional_iterator_tag,sparse_bidirectional_iterator_tag
295 ) {
296  //apply the transposition of e()
297  //first evaluate e and fill the values into a vector which is then sorted by row_major order
298  //this gives this algorithm a run time of O(eval(e)+k*log(k))
299  //where eval(e) is the time to evaluate and k*log(k) the number of non-zero elements
300  typedef typename M::value_type value_type;
301  typedef typename M::size_type size_type;
302  typedef row_major::sparse_element<value_type> Element;
303  std::vector<Element> elements;
304 
305  size_type size2 = m().size2();
306  size_type size1 = m().size1();
307  for(size_type j = 0; j != size2; ++j){
308  typename E::const_column_iterator pos_e = e().column_begin(j);
309  typename E::const_column_iterator end_e = e().column_end(j);
310  for(; pos_e != end_e; ++pos_e){
311  Element element;
312  element.i = pos_e.index();
313  element.j = j;
314  element.value = *pos_e;
315  elements.push_back(element);
316  }
317  }
318  std::sort(elements.begin(),elements.end());
319  //fill m with the contents
320  m().clear();
321  m().reserve(elements.size());//reserve a bit of space
322  for(size_type current = 0; current != elements.size();){
323  //count elements in row and reserve enough space for it
324  size_type row = elements[current].i;
325  size_type row_end = current;
326  while(row_end != elements.size() && elements[row_end].i == row)
327  ++ row_end;
328  m().reserve_row(row,row_end - current);
329 
330  //copy elements
331  typename M::row_iterator row_pos = m().row_begin(row);
332  for(; current != row_end; ++current){
333  row_pos = m().set_element(row_pos,elements[current].j,elements[current].value);
334  ++row_pos;
335  }
336  }
337 }
338 
339 //packed row_major,row_major
340 template<class M, class E,class Triangular>
341 void assign(
343  matrix_expression<E> const& e,
344  packed<row_major,Triangular>, packed<row_major,Triangular>,
345  packed_random_access_iterator_tag, packed_random_access_iterator_tag
346 ) {
347  typedef typename M::row_iterator MIter;
348  typedef typename E::const_row_iterator EIter;
349 
350  for(std::size_t i = 0; i != m().size1(); ++i){
351  MIter mpos = m().row_begin(i);
352  EIter epos = e().row_begin(i);
353  MIter mend = m().row_end(i);
354  SIZE_CHECK(mpos.index() == epos.index());
355  for(; mpos!=mend; ++mpos,++epos){
356  *mpos = *epos;
357  }
358  }
359 }
360 
361 ////packed row_major,column_major
362 //todo: this is suboptimal as we do strided access!!!!
363 template<class M, class E,class Triangular>
364 void assign(
366  matrix_expression<E> const& e,
367  packed<row_major,Triangular>, packed<column_major,Triangular>,
368  packed_random_access_iterator_tag, packed_random_access_iterator_tag
369 ) {
370  typedef typename M::row_iterator MIter;
371 
372  for(std::size_t i = 0; i != m().size1(); ++i){
373  MIter mpos = m().row_begin(i);
374  MIter mend = m().row_end(i);
375  for(; mpos!=mend; ++mpos){
376  *mpos = e()(i,mpos.index());
377  }
378  }
379 }
380 
381 //general dispatcher: if the second argument has an unknown orientation
382 // it is chosen the same as the first one
383 template<class M, class E, class TagE, class TagM>
384 void assign(
386  matrix_expression<E> const& e,
387  row_major, unknown_orientation ,TagE tagE, TagM tagM
388 ) {
389  assign(m,e,row_major(),row_major(),tagE,tagM);
390 }
391 
392 //general dispatcher: if the first argumeent is column major, we transpose the whole expression
393 //so that it is row-major, this saves us to implment everything twice.
394 template<class M, class E,class EOrientation, class TagE, class TagM>
395 void assign(
397  matrix_expression<E> const& e,
398  column_major, EOrientation,TagE tagE, TagM tagM
399 ) {
400  typedef typename EOrientation::transposed_orientation TEOrientation;
401  detail::internal_transpose_proxy<M> transM(m());
402  detail::internal_transpose_proxy<E const> transE(e());
403  assign(transM,transE,row_major(),TEOrientation(),tagE,tagM);
404 }
405 
406 //first argument is column_major->transpose to row_major
407 template<class M, class E,class EOrientation, class Triangular, class TagM, class TagE>
408 void assign(
410  matrix_expression<E> const& e,
411  packed<column_major,Triangular>, packed<EOrientation,Triangular>,
412  TagM tagM, TagE tagE
413 ) {
414  typedef typename M::orientation::transposed_orientation TMPacked;
415  typedef typename E::orientation::transposed_orientation TEPacked;
416  detail::internal_transpose_proxy<M> transM(m());
417  detail::internal_transpose_proxy<E const> transE(e());
418  assign(transM,transE,TMPacked(),TEPacked(),tagM,tagE);
419 }
420 
421 // Dispatcher
422 template<class M, class E>
424  SIZE_CHECK(m().size1() == e().size1());
425  SIZE_CHECK(m().size2() == e().size2());
426  typedef typename M::orientation MOrientation;
427  typedef typename E::orientation EOrientation;
428  typedef typename major_iterator<M>::type::iterator_category MCategory;
429  typedef typename major_iterator<E>::type::iterator_category ECategory;
430  assign(m, e, MOrientation(),EOrientation(),MCategory(), ECategory());
431 }
432 
433 
434 ///////////////////////////////////////////////////////////////////////////////////////////
435 //////Matrix Assignment With Functor implementing +=,-=...
436 ///////////////////////////////////////////////////////////////////////////////////////////
437 
438 //we only implement the case where the target is row major and than map the rest into these
439 
440 //when both are row-major we can map to vector case
441 //this is not necessarily efficient if m is sparse.
442 template<template <class, class> class F, class M, class E, class TagE, class TagM>
443 void assign(
445  matrix_expression<E> const& e,
446  row_major, row_major,TagM, TagE
447 ) {
448  for(std::size_t i = 0; i != m().size1(); ++i){
449  matrix_row<M> rowM(m(),i);
450  kernels::assign<F>(rowM,row(e,i));
451  }
452 }
453 
454 //we only need to implement the remaining versions for column major second argument
455 
456 //dense-dense case
457 template<template <class, class> class F,class M, class E>
458 void assign(
460  matrix_expression<E> const& e,
461  row_major, column_major,dense_random_access_iterator_tag, dense_random_access_iterator_tag
462 ) {
463  F<typename major_iterator<M>::type::reference, typename E::value_type> f;
464  //compute blockwise and wrelem the transposed block.
465  std::size_t const blockSize = 16;
466  typename M::value_type blockStorage[blockSize][blockSize];
467 
468  typedef typename M::size_type size_type;
469  size_type size1 = m().size1();
470  size_type size2 = m().size2();
471  for (size_type iblock = 0; iblock < size1; iblock += blockSize){
472  for (size_type jblock = 0; jblock < size2; jblock += blockSize){
473  std::size_t blockSizei = std::min(blockSize,size1-iblock);
474  std::size_t blockSizej = std::min(blockSize,size2-jblock);
475 
476  //fill the block with the values of e
477  for (size_type j = 0; j < blockSizej; ++j){
478  for (size_type i = 0; i < blockSizei; ++i){
479  blockStorage[i][j] = e()(iblock+i,jblock+j);
480  }
481  }
482 
483  //compute block values and store in m
484  for (size_type i = 0; i < blockSizei; ++i){
485  for (size_type j = 0; j < blockSizej; ++j){
486  f(m()(iblock+i,jblock+j), blockStorage[i][j]);
487  }
488  }
489  }
490  }
491 }
492 
493 //dense-sparse case
494 template<template <class, class> class F,class M, class E>
495 void assign(
497  matrix_expression<E> const& e,
498  row_major, column_major,dense_random_access_iterator_tag, sparse_bidirectional_iterator_tag
499 ) {
500  for(std::size_t i = 0; i != m().size2(); ++i){
501  matrix_column<M> columnM(m(),i);
502  kernels::assign<F>(columnM,column(e,i));
503  }
504 }
505 
506 //sparse-dense case
507 template<template <class, class> class F,class M, class E>
508 void assign(
510  matrix_expression<E> const& e,
511  row_major, column_major, sparse_bidirectional_iterator_tag, dense_random_access_iterator_tag
512 ) {
513  for(std::size_t i = 0; i != m().size1(); ++i){
514  matrix_column<M> rowM(m(),i);
515  kernels::assign<F>(rowM,row(e,i));
516  }
517 }
518 
519 //sparse-sparse
520 template<template <class, class> class F,class M, class E>
521 void assign(
523  matrix_expression<E> const& e,
524  row_major, column_major,sparse_bidirectional_iterator_tag t,sparse_bidirectional_iterator_tag
525 ) {
526  typename matrix_temporary<M>::type eTrans = e;//explicit calculation of the transpose for now
527  assign<F>(m,eTrans,row_major(),row_major(),t,t);
528  //~ F<typename M::iterator::reference, typename E::value_type> f;
529  //~ //first evaluate e and fill the values togethe into a vector which
530  //~ //is then sorted by row_major order
531  //~ //this gives this algorithm a run time of O(eval(e)+k*log(k))
532  //~ //where eval(e) is the time to evaluate and k*log(k) the number of non-zero elements
533  //~ typedef typename M::value_type value_type;
534  //~ typedef typename M::size_type size_type;
535  //~ typedef row_major::sparse_element<value_type> Element;
536  //~ std::vector<Element> elements;
537 
538  //~ size_type size2 = m().size2();
539  //~ size_type size1 = m().size1();
540  //~ for(size_type j = 0; j != size2; ++j){
541  //~ typename E::const_column_iterator pos_e = e().column_begin(j);
542  //~ typename E::const_column_iterator end_e = e().column_end(j);
543  //~ for(; pos_e != end_e; ++pos_e){
544  //~ Element element;
545  //~ element.i = pos_e.index();
546  //~ element.j = j;
547  //~ element.value = *pos_e;
548  //~ elements.push_back(element);
549  //~ }
550  //~ }
551  //~ std::sort(elements.begin(),elements.end());
552 
553 
554  //~ //assign the contents to m, applying the functor every time
555  //~ //that is we assign it for every element for e and m
556  //~ m().reserve(elements.size());//reserve a bit of space, we might need more, though.
557  //~ std::vector<Element>::const_iterator elem = elements.begin();
558  //~ std::vector<Element>::const_iterator elem_end = elements.end();
559  //~ value_type zero = value_type();
560  //~ for(size_type row = 0; row != m().size2(); ++row){
561  //~ //todo pre-reserve enough space in the row of m()
562  //~ //merge both rows with f as functor
563  //~ typename M::row_iterator it = m().row_begin(row);
564  //~ while(it != m().row_end(row) && elem != elem_end && elem->i == row) {
565  //~ size_type it_index = it.index();
566  //~ size_type elem_index = elem->j;
567  //~ if (it_index == elem_index) {
568  //~ f(*it, *elem);
569  //~ ++ elem;
570  //~ } else if (it_index < elem_index) {
571  //~ f(*it, zero);
572  //~ } else{//it_index > elem_index so insert new element in v()
573  //~ it = m().set_element(it,elem_index,zero);
574  //~ f(*it, *elem);
575  //~ ++elem;
576  //~ }
577  //~ ++it;
578  //~ }
579  //~ //apply f to remaining elemms in the row
580  //~ for(;it != v().end();++it) {
581  //~ f(*it, zero);
582  //~ }
583  //~ //add missing elements
584  //~ for(;elem != elem_end;++it,++elem) {
585  //~ it = m().set_element(it,elem.index(),zero);
586  //~ f(*it, zero);
587  //~ }
588  //~ }
589 }
590 
591 
592 //kernels for packed
593 template<template <class, class> class F, class M, class E, class Triangular>
594 void assign(
596  matrix_expression<E> const& e,
597  packed<row_major,Triangular>, packed<row_major,Triangular>
598 ) {
599  typedef typename M::row_iterator MIter;
600  typedef typename E::const_row_iterator EIter;
601  typedef F<typename MIter::reference,typename EIter::value_type> Function;
602  //there is nothing we can do, if F does not leave the non-stored elements 0
603  //this is the case for all current aissgnment functors, but you never know :)
604  BOOST_STATIC_ASSERT( Function::left_zero_identity || Function::right_zero_identity);
605 
606  Function f;
607  for(std::size_t i = 0; i != m().size1(); ++i){
608  MIter mpos = m().row_begin(i);
609  EIter epos = e().row_begin(i);
610  MIter mend = m().row_end(i);
611  SIZE_CHECK(mpos.index() == epos.index());
612  for(; mpos!=mend; ++mpos,++epos){
613  f(*mpos,*epos);
614  }
615  }
616 }
617 
618 //todo: this is suboptimal as we do strided access!!!!
619 template<template <class, class> class F, class M, class E, class Triangular>
620 void assign(
622  matrix_expression<E> const& e,
623  packed<row_major,Triangular>, packed<column_major,Triangular>
624 ) {
625  typedef typename M::row_iterator MIter;
626  typedef typename E::const_row_iterator EIter;
627  typedef F<typename MIter::reference,typename EIter::value_type> Function;
628  //there is nothing we can do, if F does not leave the non-stored elements 0
629  BOOST_STATIC_ASSERT( Function::left_zero_identity);
630 
631  Function f;
632  for(std::size_t i = 0; i != m().size1(); ++i){
633  MIter mpos = m().row_begin(i);
634  MIter mend = m().row_end(i);
635  for(; mpos!=mend; ++mpos){
636  f(*mpos,e()(i,mpos.index()));
637  }
638  }
639 }
640 
641 //second level of dispatcher dispatches to actual computation kernels
642 
643 //first standard structured matrices which only differ in row_major/column_major
644 //for these we standardize input so that the first argument is row_major. Further
645 //we choose in the case if unknown_orientation the second argument to be the same
646 //as the first. more than one dispatcher can be called every time!
647 
648 
649 //everything fulfilled -> dispatch sparse/dense computing versions
650 template<template <class, class> class F, class M, class E>
651 void assign(
653  matrix_expression<E> const& e,
654  row_major, row_major
655 ) {
656  typedef typename M::const_row_iterator::iterator_category MCategory;
657  typedef typename E::const_row_iterator::iterator_category ECategory;
658  assign<F>(m,e,row_major(),row_major(),MCategory(),ECategory());
659 }
660 //everything fulfilled -> dispatch sparse/dense computing versions
661 template<template <class, class> class F, class M, class E>
662 void assign(
664  matrix_expression<E> const& e,
665  row_major, column_major
666 ) {
667  typedef typename M::const_row_iterator::iterator_category MCategory;
668  typedef typename E::const_column_iterator::iterator_category ECategory;
669  assign<F>(m,e,row_major(),column_major(),MCategory(),ECategory());
670 }
671 
672 //first argument is row_major, second is unknown->choose unknown to be row_major
673 template<template <class, class> class F,class M, class E>
674 void assign(
676  matrix_expression<E> const& e,
677  row_major, unknown_orientation
678 ) {
679  assign<F>(m,e,row_major(),row_major());
680 }
681 
682 
683 //first argument is column_major->transpose to row_major
684 template<template <class, class> class F, class M, class E,class EOrientation>
685 void assign(
687  matrix_expression<E> const& e,
688  column_major, EOrientation
689 ) {
690  typedef typename EOrientation::transposed_orientation TEOrientation;
691  detail::internal_transpose_proxy<M> transM(m());
692  detail::internal_transpose_proxy<E const> transE(e());
693  assign<F>(transM,transE,row_major(),TEOrientation());
694 }
695 
696 //now dispatch packed matrices. Again we dispatch so that the default is row_major
697 //we also ensure here that the triangular structure is compatible
698 
699 //first argument is column_major->transpose to row_major
700 template<template <class, class> class F, class M, class E,class EOrientation, class Triangular>
701 void assign(
703  matrix_expression<E> const& e,
704  packed<column_major,Triangular>, packed<EOrientation,Triangular>
705 ) {
706  typedef typename M::orientation::transposed_orientation TMPacked;
707  typedef typename E::orientation::transposed_orientation TEPacked;
708  detail::internal_transpose_proxy<M> transM(m());
709  detail::internal_transpose_proxy<E const> transE(e());
710  assign<F>(transM,transE,TMPacked(),TEPacked());
711 }
712 
713 
714 //First Level Dispatcher, dispatches by orientation
715 template<template <class,class> class F, class M, class E>
717  SIZE_CHECK(m().size1() == e().size1());
718  SIZE_CHECK(m().size2() == e().size2());
719  typedef typename M::orientation MOrientation;
720  typedef typename E::orientation EOrientation;
721 
722  assign<F>(m, e, MOrientation(),EOrientation());
723 }
724 
725 }}}
726 
727 #endif