37 #ifndef GMM_REAL_PART_H
38 #define GMM_REAL_PART_H
46 struct linalg_real_part {};
47 struct linalg_imag_part {};
48 template <
typename R,
typename PART>
struct which_part {};
50 template <
typename C>
typename number_traits<C>::magnitude_type
51 real_or_imag_part(C x, linalg_real_part) {
return gmm::real(x); }
52 template <
typename C>
typename number_traits<C>::magnitude_type
53 real_or_imag_part(C x, linalg_imag_part) {
return gmm::imag(x); }
54 template <
typename T,
typename C,
typename OP> C
55 complex_from(T x, C y, OP op, linalg_real_part) {
return std::complex<T>(op(std::real(y), x), std::imag(y)); }
56 template <
typename T,
typename C,
typename OP> C
57 complex_from(T x, C y, OP op,linalg_imag_part) {
return std::complex<T>(std::real(y), op(std::imag(y), x)); }
59 template<
typename T>
struct project2nd {
60 T operator()(T , T b)
const {
return b; }
63 template<
typename T,
typename R,
typename PART>
class ref_elt_vector<T, which_part<R, PART> > {
69 operator T()
const {
return real_or_imag_part(std::complex<T>(r), PART()); }
70 ref_elt_vector(R r_) : r(r_) {}
71 inline ref_elt_vector &operator =(T v)
72 { r = complex_from(v, std::complex<T>(r), gmm::project2nd<T>(), PART());
return *
this; }
73 inline bool operator ==(T v)
const {
return (r == v); }
74 inline bool operator !=(T v)
const {
return (r != v); }
75 inline ref_elt_vector &operator +=(T v)
76 { r = complex_from(v, std::complex<T>(r), std::plus<T>(), PART());
return *
this; }
77 inline ref_elt_vector &operator -=(T v)
78 { r = complex_from(v, std::complex<T>(r), std::minus<T>(), PART());
return *
this; }
79 inline ref_elt_vector &operator /=(T v)
80 { r = complex_from(v, std::complex<T>(r), std::divides<T>(), PART());
return *
this; }
81 inline ref_elt_vector &operator *=(T v)
82 { r = complex_from(v, std::complex<T>(r), std::multiplies<T>(), PART());
return *
this; }
83 inline ref_elt_vector &operator =(
const ref_elt_vector &re)
84 { *
this = T(re);
return *
this; }
89 T operator *(T v) {
return T(*
this)* v; }
90 T operator /(T v) {
return T(*
this)/ v; }
93 template<
typename reference>
struct ref_or_value_type {
94 template <
typename T,
typename W>
95 static W r(
const T &x, linalg_real_part, W) {
98 template <
typename T,
typename W>
99 static W r(
const T &x, linalg_imag_part, W) {
104 template<
typename U,
typename R,
typename PART>
105 struct ref_or_value_type<ref_elt_vector<U, which_part<R, PART> > > {
106 template<
typename T ,
typename W>
107 static const T &r(
const T &x, linalg_real_part, W)
109 template<
typename T,
typename W>
110 static const T &r(
const T &x, linalg_imag_part, W) {
113 template<
typename T ,
typename W>
114 static T &r(T &x, linalg_real_part, W)
116 template<
typename T,
typename W>
117 static T &r(T &x, linalg_imag_part, W) {
127 template <
typename IT,
typename MIT,
typename PART>
128 struct part_vector_iterator {
129 typedef typename std::iterator_traits<IT>::value_type vtype;
130 typedef typename gmm::number_traits<vtype>::magnitude_type value_type;
131 typedef value_type *pointer;
132 typedef ref_elt_vector<value_type, which_part<typename std::iterator_traits<IT>::reference, PART> > reference;
133 typedef typename std::iterator_traits<IT>::difference_type difference_type;
134 typedef typename std::iterator_traits<IT>::iterator_category
139 part_vector_iterator(
void) {}
140 explicit part_vector_iterator(
const IT &i) : it(i) {}
141 part_vector_iterator(
const part_vector_iterator<MIT, MIT, PART> &i)
143 part_vector_iterator &
operator =
144 (
const part_vector_iterator<MIT, MIT, PART> &i) { it = i.it;
return *
this; }
146 size_type index(
void)
const {
return it.index(); }
147 part_vector_iterator operator ++(
int)
148 { part_vector_iterator tmp = *
this; ++it;
return tmp; }
149 part_vector_iterator operator --(
int)
150 { part_vector_iterator tmp = *
this; --it;
return tmp; }
151 part_vector_iterator &operator ++() { ++it;
return *
this; }
152 part_vector_iterator &operator --() { --it;
return *
this; }
153 part_vector_iterator &operator +=(difference_type i)
154 { it += i;
return *
this; }
155 part_vector_iterator &operator -=(difference_type i)
156 { it -= i;
return *
this; }
157 part_vector_iterator
operator +(difference_type i)
const
158 { part_vector_iterator itb = *
this;
return (itb += i); }
159 part_vector_iterator
operator -(difference_type i)
const
160 { part_vector_iterator itb = *
this;
return (itb -= i); }
161 difference_type
operator -(
const part_vector_iterator &i)
const
162 {
return difference_type(it - i.it); }
164 reference operator *()
const {
return reference(*it); }
165 reference operator [](
size_type ii)
const {
return reference(it[ii]); }
167 bool operator ==(
const part_vector_iterator &i)
const
168 {
return (i.it == it); }
169 bool operator !=(
const part_vector_iterator &i)
const
170 {
return (i.it != it); }
171 bool operator < (
const part_vector_iterator &i)
const
172 {
return (it < i.it); }
176 template <
typename PT,
typename PART>
struct part_vector {
177 typedef part_vector<PT, PART> this_type;
178 typedef typename std::iterator_traits<PT>::value_type V;
180 typedef typename select_ref<typename linalg_traits<V>::const_iterator,
181 typename linalg_traits<V>::iterator, PT>::ref_type iterator;
182 typedef typename linalg_traits<this_type>::reference reference;
183 typedef typename linalg_traits<this_type>::value_type value_type;
184 typedef typename linalg_traits<this_type>::porigin_type porigin_type;
186 iterator begin_, end_;
190 size_type size(
void)
const {
return size_; }
192 reference operator[](
size_type i)
const {
193 return reference(ref_or_value_type<reference>::r(
194 linalg_traits<V>::access(origin, begin_, end_, i),
195 PART(), value_type()));
199 : begin_(vect_begin(v)), end_(vect_end(v)),
200 origin(linalg_origin(v)), size_(gmm::vect_size(v)) {}
201 part_vector(
const V &v)
202 : begin_(vect_begin(const_cast<V &>(v))),
203 end_(vect_end(const_cast<V &>(v))),
204 origin(linalg_origin(const_cast<V &>(v))), size_(gmm::vect_size(v)) {}
206 part_vector(
const part_vector<CPT, PART> &cr)
207 : begin_(cr.begin_),end_(cr.end_),origin(cr.origin), size_(cr.size_) {}
210 template <
typename IT,
typename MIT,
typename ORG,
typename PT,
211 typename PART>
inline
212 void set_to_begin(part_vector_iterator<IT, MIT, PART> &it,
213 ORG o, part_vector<PT, PART> *, linalg_modifiable) {
214 typedef part_vector<PT, PART> VECT;
215 typedef typename linalg_traits<VECT>::V_reference ref_t;
216 set_to_begin(it.it, o,
typename linalg_traits<VECT>::pV(), ref_t());
218 template <
typename IT,
typename MIT,
typename ORG,
typename PT,
219 typename PART>
inline
220 void set_to_begin(part_vector_iterator<IT, MIT, PART> &it,
221 ORG o,
const part_vector<PT, PART> *, linalg_modifiable) {
222 typedef part_vector<PT, PART> VECT;
223 typedef typename linalg_traits<VECT>::V_reference ref_t;
224 set_to_begin(it.it, o,
typename linalg_traits<VECT>::pV(), ref_t());
226 template <
typename IT,
typename MIT,
typename ORG,
typename PT,
227 typename PART>
inline
228 void set_to_end(part_vector_iterator<IT, MIT, PART> &it,
229 ORG o, part_vector<PT, PART> *, linalg_modifiable) {
230 typedef part_vector<PT, PART> VECT;
231 typedef typename linalg_traits<VECT>::V_reference ref_t;
232 set_to_end(it.it, o,
typename linalg_traits<VECT>::pV(), ref_t());
234 template <
typename IT,
typename MIT,
typename ORG,
235 typename PT,
typename PART>
inline
236 void set_to_end(part_vector_iterator<IT, MIT, PART> &it,
237 ORG o,
const part_vector<PT, PART> *,
239 typedef part_vector<PT, PART> VECT;
240 typedef typename linalg_traits<VECT>::V_reference ref_t;
241 set_to_end(it.it, o,
typename linalg_traits<VECT>::pV(), ref_t());
244 template <
typename PT,
typename PART> std::ostream &
operator <<
245 (std::ostream &o,
const part_vector<PT, PART>& m)
246 { gmm::write(o,m);
return o; }
254 template <
typename PT,
typename PART>
struct part_row_ref {
256 typedef part_row_ref<PT, PART> this_type;
257 typedef typename std::iterator_traits<PT>::value_type M;
259 typedef typename std::iterator_traits<PT>::reference ref_M;
260 typedef typename select_ref<typename linalg_traits<this_type>
261 ::const_row_iterator,
typename linalg_traits<this_type>
262 ::row_iterator, PT>::ref_type iterator;
263 typedef typename linalg_traits<this_type>::value_type value_type;
264 typedef typename linalg_traits<this_type>::reference reference;
265 typedef typename linalg_traits<this_type>::porigin_type porigin_type;
267 iterator begin_, end_;
271 part_row_ref(ref_M m)
272 : begin_(mat_row_begin(m)), end_(mat_row_end(m)),
273 origin(linalg_origin(m)), nr(mat_nrows(m)), nc(mat_ncols(m)) {}
275 part_row_ref(
const part_row_ref<CPT, PART> &cr) :
276 begin_(cr.begin_),end_(cr.end_), origin(cr.origin),nr(cr.nr),nc(cr.nc) {}
279 return reference(ref_or_value_type<reference>::r(
280 linalg_traits<M>::access(begin_+i, j),
281 PART(), value_type()));
285 template<
typename PT,
typename PART> std::ostream &
operator <<
286 (std::ostream &o,
const part_row_ref<PT, PART>& m)
287 { gmm::write(o,m);
return o; }
289 template <
typename PT,
typename PART>
struct part_col_ref {
291 typedef part_col_ref<PT, PART> this_type;
292 typedef typename std::iterator_traits<PT>::value_type M;
294 typedef typename std::iterator_traits<PT>::reference ref_M;
295 typedef typename select_ref<typename linalg_traits<this_type>
296 ::const_col_iterator,
typename linalg_traits<this_type>
297 ::col_iterator, PT>::ref_type iterator;
298 typedef typename linalg_traits<this_type>::value_type value_type;
299 typedef typename linalg_traits<this_type>::reference reference;
300 typedef typename linalg_traits<this_type>::porigin_type porigin_type;
302 iterator begin_, end_;
306 part_col_ref(ref_M m)
307 : begin_(mat_col_begin(m)), end_(mat_col_end(m)),
308 origin(linalg_origin(m)), nr(mat_nrows(m)), nc(mat_ncols(m)) {}
310 part_col_ref(
const part_col_ref<CPT, PART> &cr) :
311 begin_(cr.begin_),end_(cr.end_), origin(cr.origin),nr(cr.nr),nc(cr.nc) {}
314 return reference(ref_or_value_type<reference>::r(
315 linalg_traits<M>::access(begin_+j, i),
316 PART(), value_type()));
322 template<
typename PT,
typename PART> std::ostream &
operator <<
323 (std::ostream &o,
const part_col_ref<PT, PART>& m)
324 { gmm::write(o,m);
return o; }
331 template <
typename TYPE,
typename PART,
typename PT>
332 struct part_return_ {
333 typedef abstract_null_type return_type;
335 template <
typename PT,
typename PART>
336 struct part_return_<row_major, PART, PT> {
337 typedef typename std::iterator_traits<PT>::value_type L;
338 typedef typename select_return<part_row_ref<const L *, PART>,
339 part_row_ref< L *, PART>, PT>::return_type return_type;
341 template <
typename PT,
typename PART>
342 struct part_return_<col_major, PART, PT> {
343 typedef typename std::iterator_traits<PT>::value_type L;
344 typedef typename select_return<part_col_ref<const L *, PART>,
345 part_col_ref<L *, PART>, PT>::return_type return_type;
348 template <
typename PT,
typename PART,
typename LT>
struct part_return__{
349 typedef abstract_null_type return_type;
352 template <
typename PT,
typename PART>
353 struct part_return__<PT, PART, abstract_matrix> {
354 typedef typename std::iterator_traits<PT>::value_type L;
355 typedef typename part_return_<
typename principal_orientation_type<
356 typename linalg_traits<L>::sub_orientation>::potype, PART,
357 PT>::return_type return_type;
360 template <
typename PT,
typename PART>
361 struct part_return__<PT, PART, abstract_vector> {
362 typedef typename std::iterator_traits<PT>::value_type L;
363 typedef typename select_return<part_vector<const L *, PART>,
364 part_vector<L *, PART>, PT>::return_type return_type;
367 template <
typename PT,
typename PART>
struct part_return {
368 typedef typename std::iterator_traits<PT>::value_type L;
369 typedef typename part_return__<PT, PART,
370 typename linalg_traits<L>::linalg_type>::return_type return_type;
373 template <
typename L>
inline
374 typename part_return<const L *, linalg_real_part>::return_type
375 real_part(
const L &l) {
376 return typename part_return<const L *, linalg_real_part>::return_type
377 (linalg_cast(
const_cast<L &
>(l)));
380 template <
typename L>
inline
381 typename part_return<L *, linalg_real_part>::return_type
383 return typename part_return<L *, linalg_real_part>::return_type(linalg_cast(l));
386 template <
typename L>
inline
387 typename part_return<const L *, linalg_imag_part>::return_type
388 imag_part(
const L &l) {
389 return typename part_return<const L *, linalg_imag_part>::return_type
390 (linalg_cast(
const_cast<L &
>(l)));
393 template <
typename L>
inline
394 typename part_return<L *, linalg_imag_part>::return_type
396 return typename part_return<L *, linalg_imag_part>::return_type(linalg_cast(l));
400 template <
typename PT,
typename PART>
401 struct linalg_traits<part_vector<PT, PART> > {
402 typedef part_vector<PT, PART> this_type;
403 typedef this_type * pthis_type;
405 typedef typename std::iterator_traits<PT>::value_type V;
406 typedef typename linalg_traits<V>::index_sorted index_sorted;
407 typedef typename linalg_traits<V>::is_reference V_reference;
408 typedef typename linalg_traits<V>::origin_type origin_type;
409 typedef typename select_ref<
const origin_type *, origin_type *,
410 PT>::ref_type porigin_type;
411 typedef typename which_reference<PT>::is_reference is_reference;
412 typedef abstract_vector linalg_type;
413 typedef typename linalg_traits<V>::value_type vtype;
414 typedef typename number_traits<vtype>::magnitude_type value_type;
415 typedef typename select_ref<value_type, ref_elt_vector<value_type,
416 which_part<typename linalg_traits<V>::reference,
417 PART> >, PT>::ref_type reference;
418 typedef typename select_ref<typename linalg_traits<V>::const_iterator,
419 typename linalg_traits<V>::iterator, PT>::ref_type pre_iterator;
420 typedef typename select_ref<abstract_null_type,
421 part_vector_iterator<pre_iterator, pre_iterator, PART>,
422 PT>::ref_type iterator;
423 typedef part_vector_iterator<typename linalg_traits<V>::const_iterator,
424 pre_iterator, PART> const_iterator;
425 typedef typename linalg_traits<V>::storage_type storage_type;
426 static size_type size(
const this_type &v) {
return v.size(); }
427 static iterator begin(this_type &v) {
428 iterator it; it.it = v.begin_;
429 if (!is_const_reference(is_reference()) && is_sparse(storage_type()))
430 set_to_begin(it, v.origin, pthis_type(), is_reference());
433 static const_iterator begin(
const this_type &v) {
434 const_iterator it(v.begin_);
435 if (!is_const_reference(is_reference()) && is_sparse(storage_type()))
436 { set_to_begin(it, v.origin, pthis_type(), is_reference()); }
439 static iterator end(this_type &v) {
441 if (!is_const_reference(is_reference()) && is_sparse(storage_type()))
442 set_to_end(it, v.origin, pthis_type(), is_reference());
445 static const_iterator end(
const this_type &v) {
446 const_iterator it(v.end_);
447 if (!is_const_reference(is_reference()) && is_sparse(storage_type()))
448 set_to_end(it, v.origin, pthis_type(), is_reference());
451 static origin_type* origin(this_type &v) {
return v.origin; }
452 static const origin_type* origin(
const this_type &v) {
return v.origin; }
454 static void clear(origin_type* o,
const iterator &begin_,
455 const iterator &end_, abstract_sparse) {
456 std::deque<size_type> ind;
457 iterator it = begin_;
458 for (; it != end_; ++it) ind.push_front(it.index());
459 for (; !(ind.empty()); ind.pop_back())
460 access(o, begin_, end_, ind.back()) = value_type(0);
462 static void clear(origin_type* o,
const iterator &begin_,
463 const iterator &end_, abstract_skyline) {
464 clear(o, begin_, end_, abstract_sparse());
466 static void clear(origin_type* o,
const iterator &begin_,
467 const iterator &end_, abstract_dense) {
468 for (iterator it = begin_; it != end_; ++it) *it = value_type(0);
471 static void clear(origin_type* o,
const iterator &begin_,
472 const iterator &end_)
473 {
clear(o, begin_, end_, storage_type()); }
474 static void do_clear(this_type &v) {
clear(v.origin, begin(v), end(v)); }
475 static value_type access(
const origin_type *o,
const const_iterator &it,
476 const const_iterator &ite,
size_type i) {
477 return real_or_imag_part(linalg_traits<V>::access(o, it.it, ite.it,i),
480 static reference access(origin_type *o,
const iterator &it,
482 {
return reference(linalg_traits<V>::access(o, it.it, ite.it,i)); }
485 template <
typename PT,
typename PART>
486 struct linalg_traits<part_row_ref<PT, PART> > {
487 typedef part_row_ref<PT, PART> this_type;
488 typedef typename std::iterator_traits<PT>::value_type M;
489 typedef typename linalg_traits<M>::origin_type origin_type;
490 typedef typename select_ref<
const origin_type *, origin_type *,
491 PT>::ref_type porigin_type;
492 typedef typename which_reference<PT>::is_reference is_reference;
493 typedef abstract_matrix linalg_type;
494 typedef typename linalg_traits<M>::value_type vtype;
495 typedef typename number_traits<vtype>::magnitude_type value_type;
496 typedef typename linalg_traits<M>::storage_type storage_type;
497 typedef abstract_null_type sub_col_type;
498 typedef abstract_null_type const_sub_col_type;
499 typedef abstract_null_type col_iterator;
500 typedef abstract_null_type const_col_iterator;
501 typedef typename org_type<typename linalg_traits<M>::const_sub_row_type>::t
502 pre_const_sub_row_type;
503 typedef typename org_type<typename linalg_traits<M>::sub_row_type>::t pre_sub_row_type;
504 typedef part_vector<const pre_const_sub_row_type *, PART>
506 typedef typename select_ref<abstract_null_type,
507 part_vector<pre_sub_row_type *, PART>, PT>::ref_type sub_row_type;
508 typedef typename linalg_traits<M>::const_row_iterator const_row_iterator;
509 typedef typename select_ref<abstract_null_type,
typename
510 linalg_traits<M>::row_iterator, PT>::ref_type row_iterator;
511 typedef typename select_ref<
512 typename linalg_traits<const_sub_row_type>::reference,
513 typename linalg_traits<sub_row_type>::reference,
514 PT>::ref_type reference;
515 typedef row_major sub_orientation;
516 typedef typename linalg_traits<M>::index_sorted index_sorted;
517 static size_type ncols(
const this_type &v) {
return v.nc; }
518 static size_type nrows(
const this_type &v) {
return v.nr; }
519 static const_sub_row_type row(
const const_row_iterator &it)
520 {
return const_sub_row_type(linalg_traits<M>::row(it)); }
521 static sub_row_type row(
const row_iterator &it)
522 {
return sub_row_type(linalg_traits<M>::row(it)); }
523 static row_iterator row_begin(this_type &m) {
return m.begin_; }
524 static row_iterator row_end(this_type &m) {
return m.end_; }
525 static const_row_iterator row_begin(
const this_type &m)
527 static const_row_iterator row_end(
const this_type &m) {
return m.end_; }
528 static origin_type* origin(this_type &v) {
return v.origin; }
529 static const origin_type* origin(
const this_type &v) {
return v.origin; }
530 static void do_clear(this_type &v);
531 static value_type access(
const const_row_iterator &itrow,
size_type i)
532 {
return real_or_imag_part(linalg_traits<M>::access(itrow, i), PART()); }
533 static reference access(
const row_iterator &itrow,
size_type i) {
534 return reference(ref_or_value_type<reference>::r(
535 linalg_traits<M>::access(itrow, i),
536 PART(), value_type()));
540 template <
typename PT,
typename PART>
541 struct linalg_traits<part_col_ref<PT, PART> > {
542 typedef part_col_ref<PT, PART> this_type;
543 typedef typename std::iterator_traits<PT>::value_type M;
544 typedef typename linalg_traits<M>::origin_type origin_type;
545 typedef typename select_ref<
const origin_type *, origin_type *,
546 PT>::ref_type porigin_type;
547 typedef typename which_reference<PT>::is_reference is_reference;
548 typedef abstract_matrix linalg_type;
549 typedef typename linalg_traits<M>::value_type vtype;
550 typedef typename number_traits<vtype>::magnitude_type value_type;
551 typedef typename linalg_traits<M>::storage_type storage_type;
552 typedef abstract_null_type sub_row_type;
553 typedef abstract_null_type const_sub_row_type;
554 typedef abstract_null_type row_iterator;
555 typedef abstract_null_type const_row_iterator;
556 typedef typename org_type<typename linalg_traits<M>::const_sub_col_type>::t
557 pre_const_sub_col_type;
558 typedef typename org_type<typename linalg_traits<M>::sub_col_type>::t pre_sub_col_type;
559 typedef part_vector<const pre_const_sub_col_type *, PART>
561 typedef typename select_ref<abstract_null_type,
562 part_vector<pre_sub_col_type *, PART>, PT>::ref_type sub_col_type;
563 typedef typename linalg_traits<M>::const_col_iterator const_col_iterator;
564 typedef typename select_ref<abstract_null_type,
typename
565 linalg_traits<M>::col_iterator, PT>::ref_type col_iterator;
566 typedef typename select_ref<
567 typename linalg_traits<const_sub_col_type>::reference,
568 typename linalg_traits<sub_col_type>::reference,
569 PT>::ref_type reference;
570 typedef col_major sub_orientation;
571 typedef typename linalg_traits<M>::index_sorted index_sorted;
572 static size_type nrows(
const this_type &v) {
return v.nr; }
573 static size_type ncols(
const this_type &v) {
return v.nc; }
574 static const_sub_col_type col(
const const_col_iterator &it)
575 {
return const_sub_col_type(linalg_traits<M>::col(it)); }
576 static sub_col_type col(
const col_iterator &it)
577 {
return sub_col_type(linalg_traits<M>::col(it)); }
578 static col_iterator col_begin(this_type &m) {
return m.begin_; }
579 static col_iterator col_end(this_type &m) {
return m.end_; }
580 static const_col_iterator col_begin(
const this_type &m)
582 static const_col_iterator col_end(
const this_type &m) {
return m.end_; }
583 static origin_type* origin(this_type &v) {
return v.origin; }
584 static const origin_type* origin(
const this_type &v) {
return v.origin; }
585 static void do_clear(this_type &v);
586 static value_type access(
const const_col_iterator &itcol,
size_type i)
587 {
return real_or_imag_part(linalg_traits<M>::access(itcol, i), PART()); }
588 static reference access(
const col_iterator &itcol,
size_type i) {
589 return reference(ref_or_value_type<reference>::r(
590 linalg_traits<M>::access(itcol, i),
591 PART(), value_type()));
595 template <
typename PT,
typename PART>
596 void linalg_traits<part_col_ref<PT, PART> >::do_clear(this_type &v) {
597 col_iterator it = mat_col_begin(v), ite = mat_col_end(v);
598 for (; it != ite; ++it)
clear(col(it));
601 template <
typename PT,
typename PART>
602 void linalg_traits<part_row_ref<PT, PART> >::do_clear(this_type &v) {
603 row_iterator it = mat_row_begin(v), ite = mat_row_end(v);
604 for (; it != ite; ++it)
clear(row(it));
608 #endif // GMM_REAL_PART_H