28 #ifndef SHARK_LINALG_BLAS_KERNELS_MATRIX_ASSIGN_HPP 29 #define SHARK_LINALG_BLAS_KERNELS_MATRIX_ASSIGN_HPP 42 class internal_transpose_proxy:
public matrix_expression<internal_transpose_proxy<M> > {
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;
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;
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;
65 explicit internal_transpose_proxy(matrix_closure_type m):
69 size_type size1()
const {
70 return m_expression.size2();
72 size_type size2()
const {
73 return m_expression.size1();
77 matrix_closure_type
const &expression()
const{
80 matrix_closure_type &expression(){
85 reference
operator()(size_type i, size_type j)
const{
86 return m_expression(j, i);
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;
95 const_row_iterator row_begin(std::size_t i)
const {
96 return m_expression.column_begin(i);
98 const_row_iterator row_end(std::size_t i)
const {
99 return m_expression.column_end(i);
101 const_column_iterator column_begin(std::size_t j)
const {
102 return m_expression.row_begin(j);
104 const_column_iterator column_end(std::size_t j)
const {
105 return m_expression.row_end(j);
108 row_iterator row_begin(std::size_t i) {
109 return m_expression.column_begin(i);
111 row_iterator row_end(std::size_t i) {
112 return m_expression.column_end(i);
114 column_iterator column_begin(std::size_t j) {
115 return m_expression.row_begin(j);
117 column_iterator column_end(std::size_t j) {
118 return m_expression.row_end(j);
121 typedef typename major_iterator<internal_transpose_proxy<M> >::type major_iterator;
123 major_iterator set_element(major_iterator pos, size_type index, value_type value){
124 return m_expression.set_element(pos,index,value);
127 major_iterator clear_range(major_iterator start, major_iterator end){
128 return m_expression.clear_range(start,end);
131 major_iterator clear_element(major_iterator elem){
132 return m_expression.clear_element(elem);
136 expression().clear();
139 void reserve(size_type non_zeros) {
140 m_expression.reserve(non_zeros);
143 void reserve_row(size_type
row, size_type non_zeros) {
144 m_expression.reserve_row(row,non_zeros);
146 void reserve_column(size_type
column, size_type non_zeros) {
147 m_expression.reserve_column(column,non_zeros);
150 matrix_closure_type m_expression;
161 template<
template <
class T1,
class T2>
class F,
class M>
164 typename M::value_type t,
167 for(std::size_t i = 0; i != m().size1(); ++i){
173 template<
template <
class T1,
class T2>
class F,
class M>
176 typename M::value_type t,
179 for(std::size_t i = 0; i != m().size2(); ++i){
181 assign<F>(columnM,t);
187 template<
template <
class T1,
class T2>
class F,
class M,
class Orientation,
class Triangular>
190 typename M::value_type t,
191 packed<Orientation,Triangular>
193 assign<F>(m,t,Orientation());
197 template<
template <
class T1,
class T2>
class F,
class M>
200 typename M::value_type t
202 typedef typename M::orientation orientation;
203 assign<F> (m, t, orientation());
212 template<
class M,
class E,
class TagE,
class TagM>
216 row_major, row_major,TagE, TagM
218 for(std::size_t i = 0; i != m().size1(); ++i){
227 template<
class M,
class E>
231 row_major, column_major,dense_random_access_iterator_tag, dense_random_access_iterator_tag
234 std::size_t
const blockSize = 16;
235 typename M::value_type blockStorage[blockSize][blockSize];
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);
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);
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];
263 template<
class M,
class E>
267 row_major, column_major,dense_random_access_iterator_tag, sparse_bidirectional_iterator_tag
269 for(std::size_t i = 0; i != m().size2(); ++i){
277 template<
class M,
class E>
281 row_major, column_major, sparse_bidirectional_iterator_tag, dense_random_access_iterator_tag
283 for(std::size_t i = 0; i != m().size1(); ++i){
290 template<
class M,
class E>
294 row_major, column_major,sparse_bidirectional_iterator_tag,sparse_bidirectional_iterator_tag
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;
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){
312 element.i = pos_e.index();
314 element.value = *pos_e;
315 elements.push_back(element);
318 std::sort(elements.begin(),elements.end());
321 m().reserve(elements.size());
322 for(size_type current = 0; current != elements.size();){
324 size_type
row = elements[current].i;
325 size_type row_end = current;
326 while(row_end != elements.size() && elements[row_end].i ==
row)
328 m().reserve_row(row,row_end - current);
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);
340 template<
class M,
class E,
class Triangular>
344 packed<row_major,Triangular>, packed<row_major,Triangular>,
345 packed_random_access_iterator_tag, packed_random_access_iterator_tag
347 typedef typename M::row_iterator MIter;
348 typedef typename E::const_row_iterator EIter;
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);
355 for(; mpos!=mend; ++mpos,++epos){
363 template<
class M,
class E,
class Triangular>
367 packed<row_major,Triangular>, packed<column_major,Triangular>,
368 packed_random_access_iterator_tag, packed_random_access_iterator_tag
370 typedef typename M::row_iterator MIter;
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());
383 template<
class M,
class E,
class TagE,
class TagM>
387 row_major, unknown_orientation ,TagE tagE, TagM tagM
389 assign(m,e,row_major(),row_major(),tagE,tagM);
394 template<
class M,
class E,
class EOrientation,
class TagE,
class TagM>
398 column_major, EOrientation,TagE tagE, TagM tagM
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);
407 template<
class M,
class E,
class EOrientation,
class Triangular,
class TagM,
class TagE>
411 packed<column_major,Triangular>, packed<EOrientation,Triangular>,
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);
422 template<
class M,
class E>
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());
442 template<
template <
class,
class>
class F,
class M,
class E,
class TagE,
class TagM>
446 row_major, row_major,TagM, TagE
448 for(std::size_t i = 0; i != m().size1(); ++i){
450 kernels::assign<F>(rowM,
row(e,i));
457 template<
template <
class,
class>
class F,
class M,
class E>
461 row_major, column_major,dense_random_access_iterator_tag, dense_random_access_iterator_tag
463 F<typename major_iterator<M>::type::reference,
typename E::value_type> f;
465 std::size_t
const blockSize = 16;
466 typename M::value_type blockStorage[blockSize][blockSize];
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);
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);
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]);
494 template<
template <
class,
class>
class F,
class M,
class E>
498 row_major, column_major,dense_random_access_iterator_tag, sparse_bidirectional_iterator_tag
500 for(std::size_t i = 0; i != m().size2(); ++i){
502 kernels::assign<F>(columnM,
column(e,i));
507 template<
template <
class,
class>
class F,
class M,
class E>
511 row_major, column_major, sparse_bidirectional_iterator_tag, dense_random_access_iterator_tag
513 for(std::size_t i = 0; i != m().size1(); ++i){
515 kernels::assign<F>(rowM,
row(e,i));
520 template<
template <
class,
class>
class F,
class M,
class E>
524 row_major, column_major,sparse_bidirectional_iterator_tag t,sparse_bidirectional_iterator_tag
526 typename matrix_temporary<M>::type eTrans = e;
527 assign<F>(m,eTrans,row_major(),row_major(),t,t);
593 template<
template <
class,
class>
class F,
class M,
class E,
class Triangular>
597 packed<row_major,Triangular>, packed<row_major,Triangular>
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;
604 BOOST_STATIC_ASSERT( Function::left_zero_identity || Function::right_zero_identity);
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);
612 for(; mpos!=mend; ++mpos,++epos){
619 template<
template <
class,
class>
class F,
class M,
class E,
class Triangular>
623 packed<row_major,Triangular>, packed<column_major,Triangular>
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;
629 BOOST_STATIC_ASSERT( Function::left_zero_identity);
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()));
650 template<
template <
class,
class>
class F,
class M,
class E>
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());
661 template<
template <
class,
class>
class F,
class M,
class E>
665 row_major, column_major
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());
673 template<
template <
class,
class>
class F,
class M,
class E>
677 row_major, unknown_orientation
679 assign<F>(m,e,row_major(),row_major());
684 template<
template <
class,
class>
class F,
class M,
class E,
class EOrientation>
688 column_major, EOrientation
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());
700 template<
template <
class,
class>
class F,
class M,
class E,
class EOrientation,
class Triangular>
704 packed<column_major,Triangular>, packed<EOrientation,Triangular>
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());
715 template<
template <
class,
class>
class F,
class M,
class E>
719 typedef typename M::orientation MOrientation;
720 typedef typename E::orientation EOrientation;
722 assign<F>(m, e, MOrientation(),EOrientation());