ViennaCL - The Vienna Computing Library  1.5.2
matrix_operations.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_MATRIX_OPERATIONS_HPP_
2 #define VIENNACL_LINALG_MATRIX_OPERATIONS_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
25 #include "viennacl/forwards.h"
26 #include "viennacl/scalar.hpp"
27 #include "viennacl/vector.hpp"
29 #include "viennacl/tools/tools.hpp"
33 #include "viennacl/traits/size.hpp"
37 #include "viennacl/vector.hpp"
39 
40 #ifdef VIENNACL_WITH_OPENCL
42 #endif
43 
44 #ifdef VIENNACL_WITH_CUDA
46 #endif
47 
48 namespace viennacl
49 {
50  namespace linalg
51  {
52 
53  template <typename NumericT, typename F,
54  typename ScalarType1>
56  matrix_base<NumericT, F> const & mat2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
57  {
58  switch (viennacl::traits::handle(mat1).get_active_handle_id())
59  {
61  viennacl::linalg::host_based::am(mat1, mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha);
62  break;
63 #ifdef VIENNACL_WITH_OPENCL
65  viennacl::linalg::opencl::am(mat1, mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha);
66  break;
67 #endif
68 #ifdef VIENNACL_WITH_CUDA
70  viennacl::linalg::cuda::am(mat1, mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha);
71  break;
72 #endif
74  throw memory_exception("not initialised!");
75  default:
76  throw memory_exception("not implemented");
77  }
78  }
79 
80 
81  template <typename NumericT, typename F,
82  typename ScalarType1, typename ScalarType2>
84  matrix_base<NumericT, F> const & mat2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
85  matrix_base<NumericT, F> const & mat3, ScalarType2 const & beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
86  {
87  switch (viennacl::traits::handle(mat1).get_active_handle_id())
88  {
91  mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
92  mat3, beta, len_beta, reciprocal_beta, flip_sign_beta);
93  break;
94 #ifdef VIENNACL_WITH_OPENCL
97  mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
98  mat3, beta, len_beta, reciprocal_beta, flip_sign_beta);
99  break;
100 #endif
101 #ifdef VIENNACL_WITH_CUDA
104  mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
105  mat3, beta, len_beta, reciprocal_beta, flip_sign_beta);
106  break;
107 #endif
109  throw memory_exception("not initialised!");
110  default:
111  throw memory_exception("not implemented");
112  }
113  }
114 
115 
116  template <typename NumericT, typename F,
117  typename ScalarType1, typename ScalarType2>
119  matrix_base<NumericT, F> const & mat2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
120  matrix_base<NumericT, F> const & mat3, ScalarType2 const & beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
121  {
122  switch (viennacl::traits::handle(mat1).get_active_handle_id())
123  {
126  mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
127  mat3, beta, len_beta, reciprocal_beta, flip_sign_beta);
128  break;
129 #ifdef VIENNACL_WITH_OPENCL
132  mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
133  mat3, beta, len_beta, reciprocal_beta, flip_sign_beta);
134  break;
135 #endif
136 #ifdef VIENNACL_WITH_CUDA
139  mat2, alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
140  mat3, beta, len_beta, reciprocal_beta, flip_sign_beta);
141  break;
142 #endif
144  throw memory_exception("not initialised!");
145  default:
146  throw memory_exception("not implemented");
147  }
148  }
149 
150 
151  template <typename NumericT, typename F>
152  void matrix_assign(matrix_base<NumericT, F> & mat, NumericT s, bool clear = false)
153  {
154  switch (viennacl::traits::handle(mat).get_active_handle_id())
155  {
158  break;
159 #ifdef VIENNACL_WITH_OPENCL
162  break;
163 #endif
164 #ifdef VIENNACL_WITH_CUDA
167  break;
168 #endif
170  throw memory_exception("not initialised!");
171  default:
172  throw memory_exception("not implemented");
173  }
174  }
175 
176 
177  template <typename NumericT, typename F>
179  {
180  switch (viennacl::traits::handle(mat).get_active_handle_id())
181  {
184  break;
185 #ifdef VIENNACL_WITH_OPENCL
188  break;
189 #endif
190 #ifdef VIENNACL_WITH_CUDA
193  break;
194 #endif
196  throw memory_exception("not initialised!");
197  default:
198  throw memory_exception("not implemented");
199  }
200  }
201 
202 
204  template <typename NumericT, typename F>
206  {
207  switch (viennacl::traits::handle(v).get_active_handle_id())
208  {
211  break;
212 #ifdef VIENNACL_WITH_OPENCL
215  break;
216 #endif
217 #ifdef VIENNACL_WITH_CUDA
220  break;
221 #endif
223  throw memory_exception("not initialised!");
224  default:
225  throw memory_exception("not implemented");
226  }
227  }
228 
230  template <typename NumericT, typename F>
232  {
233  switch (viennacl::traits::handle(A).get_active_handle_id())
234  {
237  break;
238 #ifdef VIENNACL_WITH_OPENCL
241  break;
242 #endif
243 #ifdef VIENNACL_WITH_CUDA
246  break;
247 #endif
249  throw memory_exception("not initialised!");
250  default:
251  throw memory_exception("not implemented");
252  }
253  }
254 
255  template <typename NumericT, typename F>
256  void matrix_row(const matrix_base<NumericT, F> & A, unsigned int i, vector_base<NumericT> & v)
257  {
258  switch (viennacl::traits::handle(A).get_active_handle_id())
259  {
262  break;
263 #ifdef VIENNACL_WITH_OPENCL
266  break;
267 #endif
268 #ifdef VIENNACL_WITH_CUDA
271  break;
272 #endif
274  throw memory_exception("not initialised!");
275  default:
276  throw memory_exception("not implemented");
277  }
278  }
279 
280  template <typename NumericT, typename F>
281  void matrix_column(const matrix_base<NumericT, F> & A, unsigned int j, vector_base<NumericT> & v)
282  {
283  switch (viennacl::traits::handle(A).get_active_handle_id())
284  {
287  break;
288 #ifdef VIENNACL_WITH_OPENCL
291  break;
292 #endif
293 #ifdef VIENNACL_WITH_CUDA
296  break;
297 #endif
299  throw memory_exception("not initialised!");
300  default:
301  throw memory_exception("not implemented");
302  }
303  }
304 
310  template <typename T, typename F>
312  scalar<T> & result)
313  {
314  typedef typename matrix_base<T, F>::handle_type HandleType;
315  viennacl::vector_base<T> temp(const_cast<HandleType &>(A.handle()), A.internal_size(), 0, 1);
316  norm_2_impl(temp, result);
317  }
318 
324  template <typename T, typename F>
326  T & result)
327  {
328  typedef typename matrix_base<T, F>::handle_type HandleType;
329  viennacl::vector_base<T> temp(const_cast<HandleType &>(A.handle()), A.internal_size(), 0, 1);
330  norm_2_cpu(temp, result);
331  }
332 
333  //
335  //
336 
337 
338 
339  // A * x
340 
349  template <typename NumericT, typename F>
351  const vector_base<NumericT> & vec,
352  vector_base<NumericT> & result)
353  {
354  assert( (viennacl::traits::size1(mat) == viennacl::traits::size(result)) && bool("Size check failed at v1 = prod(A, v2): size1(A) != size(v1)"));
355  assert( (viennacl::traits::size2(mat) == viennacl::traits::size(vec)) && bool("Size check failed at v1 = prod(A, v2): size2(A) != size(v2)"));
356 
357  switch (viennacl::traits::handle(mat).get_active_handle_id())
358  {
360  viennacl::linalg::host_based::prod_impl(mat, vec, result);
361  break;
362 #ifdef VIENNACL_WITH_OPENCL
364  viennacl::linalg::opencl::prod_impl(mat, vec, result);
365  break;
366 #endif
367 #ifdef VIENNACL_WITH_CUDA
369  viennacl::linalg::cuda::prod_impl(mat, vec, result);
370  break;
371 #endif
373  throw memory_exception("not initialised!");
374  default:
375  throw memory_exception("not implemented");
376  }
377  }
378 
379 
380  // trans(A) * x
381 
390  template <typename NumericT, typename F>
392  const vector_base<NumericT> & vec,
393  vector_base<NumericT> & result)
394  {
395  assert( (viennacl::traits::size1(mat_trans.lhs()) == viennacl::traits::size(vec)) && bool("Size check failed at v1 = trans(A) * v2: size1(A) != size(v2)"));
396  assert( (viennacl::traits::size2(mat_trans.lhs()) == viennacl::traits::size(result)) && bool("Size check failed at v1 = trans(A) * v2: size2(A) != size(v1)"));
397 
398  switch (viennacl::traits::handle(mat_trans.lhs()).get_active_handle_id())
399  {
401  viennacl::linalg::host_based::prod_impl(mat_trans, vec, result);
402  break;
403 #ifdef VIENNACL_WITH_OPENCL
405  viennacl::linalg::opencl::prod_impl(mat_trans, vec, result);
406  break;
407 #endif
408 #ifdef VIENNACL_WITH_CUDA
410  viennacl::linalg::cuda::prod_impl(mat_trans, vec, result);
411  break;
412 #endif
414  throw memory_exception("not initialised!");
415  default:
416  throw memory_exception("not implemented");
417  }
418  }
419 
420 
421  //
423  //
424 
430  template <typename NumericT, typename F1, typename F2, typename F3, typename ScalarType >
432  const matrix_base<NumericT, F2> & B,
434  ScalarType alpha,
435  ScalarType beta)
436  {
437  assert( (viennacl::traits::size1(A) == viennacl::traits::size1(C)) && bool("Size check failed at C = prod(A, B): size1(A) != size1(C)"));
438  assert( (viennacl::traits::size2(A) == viennacl::traits::size1(B)) && bool("Size check failed at C = prod(A, B): size2(A) != size1(B)"));
439  assert( (viennacl::traits::size2(B) == viennacl::traits::size2(C)) && bool("Size check failed at C = prod(A, B): size2(B) != size2(C)"));
440 
441 
442  switch (viennacl::traits::handle(A).get_active_handle_id())
443  {
445  viennacl::linalg::host_based::prod_impl(A, B, C, alpha, beta);
446  break;
447 #ifdef VIENNACL_WITH_OPENCL
449  viennacl::linalg::opencl::prod_impl(A, B, C, alpha, beta);
450  break;
451 #endif
452 #ifdef VIENNACL_WITH_CUDA
454  viennacl::linalg::cuda::prod_impl(A, B, C, alpha, beta);
455  break;
456 #endif
458  throw memory_exception("not initialised!");
459  default:
460  throw memory_exception("not implemented");
461  }
462  }
463 
464 
465 
471  template <typename NumericT, typename F1, typename F2, typename F3, typename ScalarType >
474  op_trans> & A,
475  const matrix_base<NumericT, F2> & B,
477  ScalarType alpha,
478  ScalarType beta)
479  {
480  assert(viennacl::traits::size2(A.lhs()) == viennacl::traits::size1(C) && bool("Size check failed at C = prod(trans(A), B): size2(A) != size1(C)"));
481  assert(viennacl::traits::size1(A.lhs()) == viennacl::traits::size1(B) && bool("Size check failed at C = prod(trans(A), B): size1(A) != size1(B)"));
482  assert(viennacl::traits::size2(B) == viennacl::traits::size2(C) && bool("Size check failed at C = prod(trans(A), B): size2(B) != size2(C)"));
483 
484  switch (viennacl::traits::handle(A.lhs()).get_active_handle_id())
485  {
487  viennacl::linalg::host_based::prod_impl(A, B, C, alpha, beta);
488  break;
489 #ifdef VIENNACL_WITH_OPENCL
491  viennacl::linalg::opencl::prod_impl(A, B, C, alpha, beta);
492  break;
493 #endif
494 #ifdef VIENNACL_WITH_CUDA
496  viennacl::linalg::cuda::prod_impl(A, B, C, alpha, beta);
497  break;
498 #endif
500  throw memory_exception("not initialised!");
501  default:
502  throw memory_exception("not implemented");
503  }
504  }
505 
506 
507 
508 
514  template <typename NumericT, typename F1, typename F2, typename F3, typename ScalarType >
518  ScalarType alpha,
519  ScalarType beta)
520  {
521  assert(viennacl::traits::size1(A) == viennacl::traits::size1(C) && bool("Size check failed at C = prod(A, trans(B)): size1(A) != size1(C)"));
522  assert(viennacl::traits::size2(A) == viennacl::traits::size2(B.lhs()) && bool("Size check failed at C = prod(A, trans(B)): size2(A) != size2(B)"));
523  assert(viennacl::traits::size1(B.lhs()) == viennacl::traits::size2(C) && bool("Size check failed at C = prod(A, trans(B)): size1(B) != size2(C)"));
524 
526  {
528  viennacl::linalg::host_based::prod_impl(A, B, C, alpha, beta);
529  break;
530 #ifdef VIENNACL_WITH_OPENCL
532  viennacl::linalg::opencl::prod_impl(A, B, C, alpha, beta);
533  break;
534 #endif
535 #ifdef VIENNACL_WITH_CUDA
537  viennacl::linalg::cuda::prod_impl(A, B, C, alpha, beta);
538  break;
539 #endif
541  throw memory_exception("not initialised!");
542  default:
543  throw memory_exception("not implemented");
544  }
545  }
546 
547 
548 
554  template <typename NumericT, typename F1, typename F2, typename F3, typename ScalarType >
558  ScalarType alpha,
559  ScalarType beta)
560  {
561  assert(viennacl::traits::size2(A.lhs()) == viennacl::traits::size1(C) && bool("Size check failed at C = prod(trans(A), trans(B)): size2(A) != size1(C)"));
562  assert(viennacl::traits::size1(A.lhs()) == viennacl::traits::size2(B.lhs()) && bool("Size check failed at C = prod(trans(A), trans(B)): size1(A) != size2(B)"));
563  assert(viennacl::traits::size1(B.lhs()) == viennacl::traits::size2(C) && bool("Size check failed at C = prod(trans(A), trans(B)): size1(B) != size2(C)"));
564 
565  switch (viennacl::traits::handle(A.lhs()).get_active_handle_id())
566  {
568  viennacl::linalg::host_based::prod_impl(A, B, C, alpha, beta);
569  break;
570 #ifdef VIENNACL_WITH_OPENCL
572  viennacl::linalg::opencl::prod_impl(A, B, C, alpha, beta);
573  break;
574 #endif
575 #ifdef VIENNACL_WITH_CUDA
577  viennacl::linalg::cuda::prod_impl(A, B, C, alpha, beta);
578  break;
579 #endif
581  throw memory_exception("not initialised!");
582  default:
583  throw memory_exception("not implemented");
584  }
585  }
586 
587 
589 
590 
591 
597  template <typename T, typename F, typename OP>
599  matrix_expression<const matrix_base<T, F>, const matrix_base<T, F>, OP> const & proxy)
600  {
601  assert( (viennacl::traits::size1(A) == viennacl::traits::size1(proxy)) && bool("Size check failed at A = element_op(B): size1(A) != size1(B)"));
602  assert( (viennacl::traits::size2(A) == viennacl::traits::size2(proxy)) && bool("Size check failed at A = element_op(B): size2(A) != size2(B)"));
603 
604  switch (viennacl::traits::handle(A).get_active_handle_id())
605  {
608  break;
609 #ifdef VIENNACL_WITH_OPENCL
612  break;
613 #endif
614 #ifdef VIENNACL_WITH_CUDA
617  break;
618 #endif
620  throw memory_exception("not initialised!");
621  default:
622  throw memory_exception("not implemented");
623  }
624  }
625 
626 
627 #define VIENNACL_MAKE_BINARY_OP(OPNAME)\
628  template <typename T, typename F>\
629  viennacl::matrix_expression<const matrix_base<T, F>, const matrix_base<T, F>, op_element_binary<op_##OPNAME> >\
630  element_##OPNAME(matrix_base<T, F> const & A, matrix_base<T, F> const & B)\
631  {\
632  return viennacl::matrix_expression<const matrix_base<T, F>, const matrix_base<T, F>, op_element_binary<op_##OPNAME> >(A, B);\
633  }\
634 \
635  template <typename M1, typename M2, typename OP, typename T, typename F>\
636  viennacl::matrix_expression<const matrix_expression<const M1, const M2, OP>,\
637  const matrix_base<T, F>,\
638  op_element_binary<op_##OPNAME> >\
639  element_##OPNAME(matrix_expression<const M1, const M2, OP> const & proxy, matrix_base<T, F> const & B)\
640  {\
641  return viennacl::matrix_expression<const matrix_expression<const M1, const M2, OP>,\
642  const matrix_base<T, F>,\
643  op_element_binary<op_##OPNAME> >(proxy, B);\
644  }\
645 \
646  template <typename T, typename F, typename M2, typename M3, typename OP>\
647  viennacl::matrix_expression<const matrix_base<T, F>,\
648  const matrix_expression<const M2, const M3, OP>,\
649  op_element_binary<op_##OPNAME> >\
650  element_##OPNAME(matrix_base<T, F> const & A, matrix_expression<const M2, const M3, OP> const & proxy)\
651  {\
652  return viennacl::matrix_expression<const matrix_base<T, F>,\
653  const matrix_expression<const M2, const M3, OP>,\
654  op_element_binary<op_##OPNAME> >(A, proxy);\
655  }\
656 \
657  template <typename M1, typename M2, typename OP1,\
658  typename M3, typename M4, typename OP2>\
659  viennacl::matrix_expression<const matrix_expression<const M1, const M2, OP1>,\
660  const matrix_expression<const M3, const M4, OP2>,\
661  op_element_binary<op_##OPNAME> >\
662  element_##OPNAME(matrix_expression<const M1, const M2, OP1> const & proxy1,\
663  matrix_expression<const M3, const M4, OP2> const & proxy2)\
664  {\
665  return viennacl::matrix_expression<const matrix_expression<const M1, const M2, OP1>,\
666  const matrix_expression<const M3, const M4, OP2>,\
667  op_element_binary<op_##OPNAME> >(proxy1, proxy2);\
668  }
669 
673 
674 #undef VIENNACL_GENERATE_BINARY_OP_OVERLOADS
675 
676 
677 
678 #define VIENNACL_MAKE_UNARY_ELEMENT_OP(funcname) \
679  template <typename T, typename F> \
680  viennacl::matrix_expression<const matrix_base<T, F>, const matrix_base<T, F>, op_element_unary<op_##funcname> > \
681  element_##funcname(matrix_base<T, F> const & A) \
682  { \
683  return viennacl::matrix_expression<const matrix_base<T, F>, const matrix_base<T, F>, op_element_unary<op_##funcname> >(A, A); \
684  } \
685  template <typename LHS, typename RHS, typename OP> \
686  viennacl::matrix_expression<const matrix_expression<const LHS, const RHS, OP>, \
687  const matrix_expression<const LHS, const RHS, OP>, \
688  op_element_unary<op_##funcname> > \
689  element_##funcname(matrix_expression<const LHS, const RHS, OP> const & proxy) \
690  { \
691  return viennacl::matrix_expression<const matrix_expression<const LHS, const RHS, OP>, \
692  const matrix_expression<const LHS, const RHS, OP>, \
693  op_element_unary<op_##funcname> >(proxy, proxy); \
694  } \
695 
713 
714 #undef VIENNACL_MAKE_UNARY_ELEMENT_OP
715 
716 
717  //
719  //
720 
721 
727  template <typename NumericT>
728  viennacl::matrix_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_prod>
729  outer_prod(const vector_base<NumericT> & vec1, const vector_base<NumericT> & vec2)
730  {
731  return viennacl::matrix_expression< const vector_base<NumericT>, const vector_base<NumericT>, op_prod>(vec1, vec2);
732  }
733 
734 
747  template <typename NumericT, typename F, typename S1>
749  S1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
750  const vector_base<NumericT> & vec1,
751  const vector_base<NumericT> & vec2)
752  {
753  switch (viennacl::traits::handle(mat1).get_active_handle_id())
754  {
757  alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
758  vec1, vec2);
759  break;
760 #ifdef VIENNACL_WITH_OPENCL
763  alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
764  vec1, vec2);
765  break;
766 #endif
767 #ifdef VIENNACL_WITH_CUDA
770  alpha, len_alpha, reciprocal_alpha, flip_sign_alpha,
771  vec1, vec2);
772  break;
773 #endif
775  throw memory_exception("not initialised!");
776  default:
777  throw memory_exception("not implemented");
778  }
779  }
780 
781  } //namespace linalg
782 
783 
784 
785 
786  //
788  //
789 
790 
791  //v += A * x
797  template <typename NumericT, typename F>
798  vector<NumericT>
801  {
802  assert(viennacl::traits::size1(proxy.lhs()) == v1.size() && bool("Size check failed for v1 += A * v2: size1(A) != size(v1)"));
803 
804  vector<NumericT> result(viennacl::traits::size1(proxy.lhs()));
805  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
806  v1 += result;
807  return v1;
808  }
809 
815  template <typename NumericT, typename F>
816  vector<NumericT>
819  {
820  assert(viennacl::traits::size1(proxy.lhs()) == v1.size() && bool("Size check failed for v1 -= A * v2: size1(A) != size(v1)"));
821 
822  vector<NumericT> result(viennacl::traits::size1(proxy.lhs()));
823  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
824  v1 -= result;
825  return v1;
826  }
827 
828 
829 
830 
831 
832  //free functions:
838  template <typename NumericT, typename F>
842  {
843  assert(viennacl::traits::size1(proxy.lhs()) == viennacl::traits::size(v1) && bool("Size check failed for v1 + A * v2: size1(A) != size(v1)"));
844 
846  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
847  result += v1;
848  return result;
849  }
850 
856  template <typename NumericT, typename F>
860  {
861  assert(viennacl::traits::size1(proxy.lhs()) == viennacl::traits::size(v1) && bool("Size check failed for v1 - A * v2: size1(A) != size(v1)"));
862 
864  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
865  result = v1 - result;
866  return result;
867  }
868 
869 
871 
872 
873  //v += A^T * x
879  template <typename NumericT, typename F>
880  vector<NumericT>
883  const vector_base<NumericT>,
884  op_prod> & proxy)
885  {
886  assert(viennacl::traits::size2(proxy.lhs()) == v1.size() && bool("Size check failed in v1 += trans(A) * v2: size2(A) != size(v1)"));
887 
888  vector<NumericT> result(viennacl::traits::size2(proxy.lhs()));
889  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
890  v1 += result;
891  return v1;
892  }
893 
894  //v -= A^T * x
900  template <typename NumericT, typename F>
901  vector<NumericT>
904  const vector_base<NumericT>,
905  op_prod> & proxy)
906  {
907  assert(viennacl::traits::size2(proxy.lhs()) == v1.size() && bool("Size check failed in v1 += trans(A) * v2: size2(A) != size(v1)"));
908 
909  vector<NumericT> result(viennacl::traits::size2(proxy.lhs()));
910  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
911  v1 -= result;
912  return v1;
913  }
914 
915 
916  //free functions:
922  template <typename NumericT, typename F>
923  vector<NumericT>
926  const vector_base<NumericT>,
927  op_prod> & proxy)
928  {
929  assert(viennacl::traits::size2(proxy.lhs()) == viennacl::traits::size(v1) && bool("Size check failed in v1 + trans(A) * v2: size2(A) != size(v1)"));
930 
932  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
933  result += v1;
934  return result;
935  }
936 
942  template <typename NumericT, typename F>
943  vector<NumericT>
946  const vector_base<NumericT>,
947  op_prod> & proxy)
948  {
949  assert(viennacl::traits::size2(proxy.lhs()) == viennacl::traits::size(v1) && bool("Size check failed in v1 - trans(A) * v2: size2(A) != size(v1)"));
950 
952  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
953  result = v1 - result;
954  return result;
955  }
956 
957 
958 } //namespace viennacl
959 
960 
961 #endif
Implementations of dense matrix related operations, including matrix-vector products, using a plain single-threaded or OpenMP-enabled execution on CPU.
void ambm_m(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:118
std::size_t vcl_size_t
Definition: forwards.h:58
size_type internal_size() const
Returns the total amount of allocated memory in multiples of sizeof(SCALARTYPE)
Definition: matrix.hpp:651
void matrix_column(const matrix_base< NumericT, F > &mat, unsigned int j, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:566
void ambm_m(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:148
Implementations of dense matrix related operations, including matrix-vector products, using OpenCL.
Definition: forwards.h:478
void scaled_rank_1_update(matrix_base< NumericT, F > &mat1, S1 const &alpha, vcl_size_t, bool reciprocal_alpha, bool flip_sign_alpha, const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
The implementation of the operation mat += alpha * vec1 * vec2^T, i.e. a scaled rank 1 update...
Definition: matrix_operations.hpp:1120
Exception class in case of memory errors.
Definition: forwards.h:485
Generic size and resize functionality for different vector and matrix types.
Extracts the underlying OpenCL start index handle from a vector, a matrix, an expression etc...
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:172
void matrix_diag_to_vector(const matrix_base< NumericT, F > &A, int k, vector_base< NumericT > &v)
Dispatcher interface for v = diag(A, k)
Definition: matrix_operations.hpp:231
Various little tools used here and there in ViennaCL.
void matrix_row(const matrix_base< NumericT, F > &mat, unsigned int i, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:539
void element_op(matrix_base< NumericT, F > &A, matrix_expression< const matrix_base< NumericT, F >, const matrix_base< NumericT, F >, op_element_binary< OP > > const &proxy)
Implementation of the element-wise operations A = B .* C and A = B ./ C (using MATLAB syntax) ...
Definition: matrix_operations.hpp:604
void matrix_column(const matrix_base< NumericT, F > &mat, unsigned int j, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:406
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
void matrix_assign(matrix_base< NumericT, F > &mat, NumericT s, bool clear=false)
Definition: matrix_operations.hpp:395
void am(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
Definition: matrix_operations.hpp:55
void ambm(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:83
void matrix_diagonal_assign(matrix_base< NumericT, F > &mat, NumericT s)
Definition: matrix_operations.hpp:289
Common base class for dense vectors, vector ranges, and vector slices.
Definition: forwards.h:205
void clear(VectorType &vec)
Generic routine for setting all entries of a vector to zero. This is the version for non-ViennaCL obj...
Definition: clear.hpp:57
This file provides the forward declarations for the main types used within ViennaCL.
Determines row and column increments for matrices and matrix proxies.
void matrix_row(const matrix_base< NumericT, F > &A, unsigned int i, vector_base< NumericT > &v)
Definition: matrix_operations.hpp:256
void prod_impl(const matrix_base< NumericT, F > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:737
void matrix_diag_from_vector(const vector_base< NumericT > &vec, int k, matrix_base< NumericT, F > &mat)
Definition: matrix_operations.hpp:318
void matrix_assign(matrix_base< NumericT, F > &mat, NumericT s, bool clear=false)
Definition: matrix_operations.hpp:257
void matrix_assign(matrix_base< NumericT, F > &mat, NumericT s, bool clear=false)
Definition: matrix_operations.hpp:152
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:245
Definition: forwards.h:481
void matrix_diag_from_vector(const vector_base< NumericT > &vec, int k, matrix_base< NumericT, F > &mat)
Definition: matrix_operations.hpp:466
memory_types get_active_handle_id() const
Returns an ID for the currently active memory buffer. Other memory buffers might contain old or no da...
Definition: mem_handle.hpp:91
A dense matrix class.
Definition: forwards.h:290
void matrix_diag_to_vector(const matrix_base< NumericT, F > &mat, int k, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:504
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
void matrix_diag_from_vector(const vector_base< NumericT > &vec, int k, matrix_base< NumericT, F > &mat)
Definition: matrix_operations.hpp:245
void am(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
Definition: matrix_operations.hpp:66
void ambm_m(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:182
void scaled_rank_1_update(matrix_base< NumericT, F > &mat1, S1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
The implementation of the operation mat += alpha * vec1 * vec2^T, i.e. a scaled rank 1 update...
Definition: matrix_operations.hpp:2470
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
vector< NumericT > operator-=(vector_base< NumericT > &v1, const viennacl::vector_expression< const matrix_base< NumericT, F >, const vector_base< NumericT >, viennacl::op_prod > &proxy)
Implementation of the operation v1 -= A * v2, where A is a matrix.
Definition: matrix_operations.hpp:817
void element_op(matrix_base< T, F > &A, matrix_expression< const matrix_base< T, F >, const matrix_base< T, F >, op_element_binary< OP > > const &proxy)
Definition: matrix_operations.hpp:489
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
Definition: forwards.h:480
void scaled_rank_1_update(matrix_base< NumericT, F > &mat1, S1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
The implementation of the operation mat += alpha * vec1 * vec2^T, i.e. a scaled rank 1 update...
Definition: matrix_operations.hpp:748
void norm_frobenius_cpu(matrix_base< T, F > const &vec, T &result)
Computes the Frobenius norm of a vector with final reduction on the CPU.
Definition: matrix_operations.hpp:325
#define VIENNACL_MAKE_UNARY_ELEMENT_OP(funcname)
Definition: matrix_operations.hpp:678
void element_op(matrix_base< T, F > &A, matrix_expression< const matrix_base< T, F >, const matrix_base< T, F >, op_element_binary< OP > > const &proxy)
Implementation of binary element-wise operations A = OP(B,C)
Definition: matrix_operations.hpp:460
vector< NumericT > operator+=(vector_base< NumericT > &v1, const viennacl::vector_expression< const matrix_base< NumericT, F >, const vector_base< NumericT >, viennacl::op_prod > &proxy)
Implementation of the operation v1 += A * v2, where A is a matrix.
Definition: matrix_operations.hpp:799
void matrix_assign(matrix_base< NumericT, F > &mat, NumericT s, bool clear=false)
Definition: matrix_operations.hpp:198
void matrix_row(const matrix_base< NumericT, F > &mat, unsigned int i, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:363
An expression template class that represents a binary operation that yields a vector.
Definition: forwards.h:181
#define VIENNACL_MAKE_BINARY_OP(OPNAME)
Definition: matrix_operations.hpp:627
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector.hpp:837
void matrix_diagonal_assign(matrix_base< NumericT, F > &mat, NumericT s)
Definition: matrix_operations.hpp:223
void element_op(matrix_base< T, F > &A, matrix_expression< const matrix_base< T, F >, const matrix_base< T, F >, OP > const &proxy)
Implementation of the element-wise operation A = B .* C and A = B ./ C for matrices (using MATLAB syn...
Definition: matrix_operations.hpp:598
Expression template class for representing a tree of expressions which ultimately result in a matrix...
Definition: forwards.h:283
void norm_2_cpu(vector_base< T > const &vec, T &result)
Computes the l^2-norm of a vector with final reduction on the CPU - dispatcher interface.
Definition: vector_operations.hpp:669
handle_type & handle()
Returns the OpenCL handle, non-const-version.
Definition: matrix.hpp:654
viennacl::vector< NumericT > operator-(const vector_base< NumericT > &v1, const vector_expression< const matrix_base< NumericT, F >, const vector_base< NumericT >, op_prod > &proxy)
Implementation of the operation 'result = v1 - A * v2', where A is a matrix.
Definition: matrix_operations.hpp:858
Proxy classes for vectors.
All the predicates used within ViennaCL. Checks for expressions to be vectors, etc.
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
viennacl::matrix_expression< const vector_base< NumericT >, const vector_base< NumericT >, op_prod > outer_prod(const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
Returns a proxy class for the operation mat += vec1 * vec2^T, i.e. a rank 1 update.
Definition: matrix_operations.hpp:729
A tag class representing matrix-vector products and element-wise multiplications. ...
Definition: forwards.h:76
Implementations of dense matrix related operations, including matrix-vector products, using CUDA.
void matrix_column(const matrix_base< NumericT, F > &A, unsigned int j, vector_base< NumericT > &v)
Definition: matrix_operations.hpp:281
void matrix_diag_to_vector(const matrix_base< NumericT, F > &mat, int k, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:370
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
void scaled_rank_1_update(matrix_base< NumericT, F > &mat1, S1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
The implementation of the operation mat += alpha * vec1 * vec2^T, i.e. a scaled rank 1 update...
Definition: matrix_operations.hpp:955
void am(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t, bool reciprocal_alpha, bool flip_sign_alpha)
Definition: matrix_operations.hpp:52
Definition: forwards.h:479
void matrix_diag_from_vector(const vector_base< NumericT > &v, int k, matrix_base< NumericT, F > &A)
Dispatcher interface for A = diag(v, k)
Definition: matrix_operations.hpp:205
void norm_2_impl(vector_base< T > const &vec, scalar< T > &result)
Computes the l^2-norm of a vector - dispatcher interface.
Definition: vector_operations.hpp:624
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
A tag class representing transposed matrices.
Definition: forwards.h:165
void prod_impl(const matrix_base< NumericT, F > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:547
void prod_impl(const matrix_base< NumericT, F > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:350
void norm_frobenius_impl(matrix_base< T, F > const &vec, scalar< T > &result)
Computes the Frobenius norm of a matrix - dispatcher interface.
Definition: matrix_operations.hpp:311
void matrix_diagonal_assign(matrix_base< NumericT, F > &mat, NumericT s)
Definition: matrix_operations.hpp:440
Extracts the underlying OpenCL handle from a vector, a matrix, an expression etc. ...
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
void am(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
Definition: matrix_operations.hpp:58
void matrix_row(const matrix_base< NumericT, F > &mat, unsigned int i, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:417
Implementation of the ViennaCL scalar class.
void ambm(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:108
viennacl::vector< NumericT > operator+(const vector_base< NumericT > &v1, const vector_expression< const matrix_base< NumericT, F >, const vector_base< NumericT >, op_prod > &proxy)
Implementation of the operation 'result = v1 + A * v2', where A is a matrix.
Definition: matrix_operations.hpp:840
void ambm_m(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:261
A collection of compile time type deductions.
void matrix_diagonal_assign(matrix_base< NumericT, F > &mat, NumericT s)
Definition: matrix_operations.hpp:178
void matrix_column(const matrix_base< NumericT, F > &mat, unsigned int j, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:450
void ambm(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:132
void ambm(matrix_base< NumericT, F > &mat1, matrix_base< NumericT, F > const &mat2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, matrix_base< NumericT, F > const &mat3, ScalarType2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: matrix_operations.hpp:98
void prod_impl(const matrix_base< NumericT, F > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:1364
Simple enable-if variant that uses the SFINAE pattern.
void matrix_diag_to_vector(const matrix_base< NumericT, F > &mat, int k, vector_base< NumericT > &vec)
Definition: matrix_operations.hpp:306