matrix_expression.hpp
Go to the documentation of this file.
1/*!
2 * \brief Expression templates for expressions involving matrices
3 *
4 * \author O. Krause
5 * \date 2016
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 REMORA_MATRIX_EXPRESSION_HPP
29#define REMORA_MATRIX_EXPRESSION_HPP
30
31#include "detail/expression_optimizers.hpp"
32#include "kernels/matrix_fold.hpp"
33#include "proxy_expressions.hpp"
34#include "vector_expression.hpp"
35//~ #include "vector_set_expressions.hpp"
36
37namespace remora{
38
39/////////////////////////////////////////////
40//////////Vector->Matrix Operations
41/////////////////////////////////////////////
42
43
44///\brief Creates a matrix from a vector by repeating the vector in every row of the matrix.
45///
46///TODO: cpu only!
47///example: vector = (1,2,3)
48///repeat(vector,3) results in
49///(1,2,3)
50///(1,2,3)
51///(1,2,3)
52///@param vector the vector which is to be repeated as the rows of the resulting matrix
53///@param rows the number of rows of the matrix
54template<class VecV, class Device>
55vector_repeater<VecV> repeat(vector_expression<VecV, Device> const& vector, std::size_t rows){
56 return vector_repeater<VecV>(vector(),rows);
57}
58
59/// \brief Repeats a single element to form a matrix of size rows x columns.
60///
61///TODO: cpu only!
62///@param scalar the value which is repeated
63///@param rows the number of rows of the resulting vector
64///@param columns the number of columns of the resulting vector
65template<class T>
66typename std::enable_if<std::is_arithmetic<T>::value, scalar_matrix<T,cpu_tag, row_major> >::type
67repeat(T scalar, std::size_t rows, std::size_t columns){
68 return scalar_matrix<T,cpu_tag, row_major>(rows, columns, scalar);
69}
70
71/// \brief An identity matrix with values of type \c T
72///
73/// Elements or cordinates \f$(i,i)\f$ are equal to 1 (one) and all others to 0 (zero).
74template<class T>
75class identity_matrix: public diagonal_matrix<scalar_vector<T, cpu_tag> > {
76 typedef diagonal_matrix<scalar_vector<T, cpu_tag> > base_type;
77public:
78 identity_matrix(){}
79 identity_matrix(std::size_t size):base_type(scalar_vector<T, cpu_tag>(size,T(1))){}
80};
81
82
83
84/// \brief Creates a nxn diagonal matrix with with diagonal given by a vector.
85template<class MatA, class Device>
86diagonal_matrix<MatA> to_diagonal(vector_expression<MatA, Device> const& v){
87 return diagonal_matrix<MatA>(v());
88}
89
90/////////////////////////////////////////////
91//////////Unary Matrix Transformations
92/////////////////////////////////////////////
93
94
95/// \brief Negates the matrix-expression A.
96///
97/// \f$ (-A)_{ij} = - e_{ij} \f$
98template<class MatA, class Device>
99typename detail::matrix_scalar_multiply_optimizer<MatA>::type
100operator-(matrix_expression<MatA, Device> const& A){
101 return detail::matrix_scalar_multiply_optimizer<MatA>::create(A(), typename MatA::value_type(-1));
102}
103
104#define REMORA_UNARY_MATRIX_TRANSFORMATION(name, F)\
105template<class MatA, class Device>\
106typename detail::matrix_unary_optimizer<MatA,typename device_traits<Device>:: template F<typename MatA::value_type> >::type \
107name(matrix_expression<MatA, Device> const& m){\
108 typedef typename device_traits<Device>:: template F<typename MatA::value_type> functor_type;\
109 return detail::matrix_unary_optimizer<MatA, functor_type >::create(m(), functor_type());\
110}
111
112REMORA_UNARY_MATRIX_TRANSFORMATION(abs, abs)
113REMORA_UNARY_MATRIX_TRANSFORMATION(log, log)
114REMORA_UNARY_MATRIX_TRANSFORMATION(exp, exp)
115REMORA_UNARY_MATRIX_TRANSFORMATION(tanh,tanh)
116REMORA_UNARY_MATRIX_TRANSFORMATION(sin,sin)
117REMORA_UNARY_MATRIX_TRANSFORMATION(cos,cos)
118REMORA_UNARY_MATRIX_TRANSFORMATION(tan,tan)
119REMORA_UNARY_MATRIX_TRANSFORMATION(asin,asin)
120REMORA_UNARY_MATRIX_TRANSFORMATION(acos,acos)
121REMORA_UNARY_MATRIX_TRANSFORMATION(atan,atan)
122REMORA_UNARY_MATRIX_TRANSFORMATION(erf,erf)
123REMORA_UNARY_MATRIX_TRANSFORMATION(erfc,erfc)
124REMORA_UNARY_MATRIX_TRANSFORMATION(sqr, sqr)
125REMORA_UNARY_MATRIX_TRANSFORMATION(sqrt, sqrt)
126REMORA_UNARY_MATRIX_TRANSFORMATION(cbrt, cbrt)
127REMORA_UNARY_MATRIX_TRANSFORMATION(sigmoid, sigmoid)
128REMORA_UNARY_MATRIX_TRANSFORMATION(softPlus, soft_plus)
129REMORA_UNARY_MATRIX_TRANSFORMATION(elem_inv, inv)
130#undef REMORA_UNARY_MATRIX_TRANSFORMATION
131
132#define REMORA_UNARY_VECTOR_SET_TRANSFORMATION(name)\
133template<class S, class Device>\
134auto name(vector_set_expression<S, Device> const& set)\
135-> decltype(as_set(name( set().expression()), typename S::point_orientation())){\
136 return as_set(name( set().expression()), typename S::point_orientation());\
137}
138REMORA_UNARY_VECTOR_SET_TRANSFORMATION(abs)
139REMORA_UNARY_VECTOR_SET_TRANSFORMATION(log)
140REMORA_UNARY_VECTOR_SET_TRANSFORMATION(exp)
141REMORA_UNARY_VECTOR_SET_TRANSFORMATION(tanh)
142REMORA_UNARY_VECTOR_SET_TRANSFORMATION(sin)
143REMORA_UNARY_VECTOR_SET_TRANSFORMATION(cos)
144REMORA_UNARY_VECTOR_SET_TRANSFORMATION(tan)
145REMORA_UNARY_VECTOR_SET_TRANSFORMATION(asin)
146REMORA_UNARY_VECTOR_SET_TRANSFORMATION(acos)
147REMORA_UNARY_VECTOR_SET_TRANSFORMATION(atan)
148REMORA_UNARY_VECTOR_SET_TRANSFORMATION(erf)
149REMORA_UNARY_VECTOR_SET_TRANSFORMATION(erfc)
150REMORA_UNARY_VECTOR_SET_TRANSFORMATION(sqr)
151REMORA_UNARY_VECTOR_SET_TRANSFORMATION(sqrt)
152REMORA_UNARY_VECTOR_SET_TRANSFORMATION(cbrt)
153REMORA_UNARY_VECTOR_SET_TRANSFORMATION(sigmoid)
154REMORA_UNARY_VECTOR_SET_TRANSFORMATION(softPlus)
155REMORA_UNARY_VECTOR_SET_TRANSFORMATION(elem_inv)
156REMORA_UNARY_VECTOR_SET_TRANSFORMATION(operator-)
157#undef REMORA_UNARY_VECTOR_SET_TRANSFORMATION
158
159
160/////////////////////////////////////////////
161//////////Matrix-Scalar Operations
162/////////////////////////////////////////////
163
164/// \brief Computes the multiplication of a matrix-expression A with a scalar t.
165///
166/// \f$ (A*t)_{ij} = e_{ij}*t \f$
167template<class MatA, class T, class Device>
168typename std::enable_if<
169 std::is_convertible<T, typename MatA::value_type >::value,
170 typename detail::matrix_scalar_multiply_optimizer<MatA>::type
171>::type
172operator* (matrix_expression<MatA, Device> const& A, T scalar){
173 return detail::matrix_scalar_multiply_optimizer<MatA>::create(A(), typename MatA::value_type(scalar));
174}
175
176/// \brief Computes the multiplication of a matrix-expression A with a scalar t.
177///
178/// \f$ (t*A)_{ij} = t*e_{ij} \f$
179template<class T, class MatA, class Device>
180typename std::enable_if<
181 std::is_convertible<T, typename MatA::value_type >::value,
182 typename detail::matrix_scalar_multiply_optimizer<MatA>::type
183>::type
184operator* (T scalar, matrix_expression<MatA, Device> const& A){
185 return detail::matrix_scalar_multiply_optimizer<MatA>::create(A(), typename MatA::value_type(scalar));
186}
187
188
189///\brief Adds a matrix plus a scalar which is interpreted as a constant matrix
190template<class MatA, class T, class Device>
191typename std::enable_if<
192 std::is_convertible<T, typename MatA::value_type>::value,
193 matrix_addition<MatA, scalar_matrix<T,Device, typename MatA::orientation> >
194>::type operator+ (
195 matrix_expression<MatA, Device> const& A,
196 T t
197){
198 return A + scalar_matrix<T,Device, typename MatA::orientation>(A().size1(),A().size2(),t);
199}
200
201///\brief Adds a matrix plus a scalar which is interpreted as a constant matrix
202template<class T, class MatA, class Device>
203typename std::enable_if<
204 std::is_convertible<T, typename MatA::value_type>::value,
205 matrix_addition<MatA, scalar_matrix<T,Device, typename MatA::orientation> >
206>::type operator+ (
207 T t,
208 matrix_expression<MatA, Device> const& A
209){
210 return A + scalar_matrix<T,Device, typename MatA::orientation>(A().size1(),A().size2(),t);
211}
212
213///\brief Subtracts a scalar which is interpreted as a constant matrix from a matrix.
214template<class MatA, class T, class Device>
215typename std::enable_if<
216 std::is_convertible<T, typename MatA::value_type>::value ,
217 decltype(std::declval<MatA const&>() + T())
218>::type operator- (
219 matrix_expression<MatA, Device> const& A,
220 T t
221){
222 return A + (-t);
223}
224
225///\brief Subtracts a matrix from a scalar which is interpreted as a constant matrix
226template<class MatA, class T, class Device>
227typename std::enable_if<
228 std::is_convertible<T, typename MatA::value_type>::value,
229 decltype(T() + (-std::declval<MatA const&>()))
230>::type operator- (
231 T t,
232 matrix_expression<MatA, Device> const& A
233){
234 return t + (-A);
235}
236
237#define REMORA_MATRIX_SCALAR_TRANSFORMATION(name, F)\
238template<class T, class MatA, class Device> \
239typename std::enable_if< \
240 std::is_convertible<T, typename MatA::value_type >::value,\
241 matrix_binary<MatA, scalar_matrix<typename MatA::value_type,Device, typename MatA::orientation>,typename device_traits<Device>:: template F<typename MatA::value_type> > \
242>::type \
243name (matrix_expression<MatA, Device> const& m, T t){ \
244 typedef typename MatA::value_type type;\
245 typedef typename device_traits<Device>:: template F<type> functor_type;\
246 typedef scalar_matrix<type,Device, typename MatA::orientation> mat_type;\
247 return matrix_binary<MatA, mat_type, functor_type >(m(), mat_type(m().size1(), m().size2(), t) ,functor_type()); \
248}
249REMORA_MATRIX_SCALAR_TRANSFORMATION(operator/, divide)
250REMORA_MATRIX_SCALAR_TRANSFORMATION(operator<, less)
251REMORA_MATRIX_SCALAR_TRANSFORMATION(operator<=, less_equal)
252REMORA_MATRIX_SCALAR_TRANSFORMATION(operator>, greater)
253REMORA_MATRIX_SCALAR_TRANSFORMATION(operator>=, greater_equal)
254REMORA_MATRIX_SCALAR_TRANSFORMATION(operator==, equal)
255REMORA_MATRIX_SCALAR_TRANSFORMATION(operator!=, not_equal)
256REMORA_MATRIX_SCALAR_TRANSFORMATION(min, min)
257REMORA_MATRIX_SCALAR_TRANSFORMATION(max, max)
258REMORA_MATRIX_SCALAR_TRANSFORMATION(pow, pow)
259#undef REMORA_MATRIX_SCALAR_TRANSFORMATION
260
261// operations of the form op(t,v)[i,j] = op(t,v[i,j])
262#define REMORA_MATRIX_SCALAR_TRANSFORMATION_2(name, F)\
263template<class T, class MatA, class Device> \
264typename std::enable_if< \
265 std::is_convertible<T, typename MatA::value_type >::value,\
266 matrix_binary<scalar_matrix< typename MatA::value_type,Device, typename MatA::orientation>, MatA, typename device_traits<Device>:: template F< typename MatA::value_type> > \
267>::type \
268name (T t, matrix_expression<MatA, Device> const& m){ \
269 typedef typename MatA::value_type type;\
270 typedef typename device_traits<Device>:: template F<type> functor_type;\
271 typedef scalar_matrix<type,Device, typename MatA::orientation> mat_type;\
272 return matrix_binary<mat_type, MatA, functor_type >(mat_type(m().size1(), m().size2(), t), m(), functor_type()); \
273}
274REMORA_MATRIX_SCALAR_TRANSFORMATION_2(min, min)
275REMORA_MATRIX_SCALAR_TRANSFORMATION_2(max, max)
276#undef REMORA_MATRIX_SCALAR_TRANSFORMATION_2
277
278
279/////////////////////////////////////////////
280//////////Simple Matrix-Binary Operations
281/////////////////////////////////////////////
282
283
284///\brief Adds two Matrices
285template<class MatA, class MatB, class Device>
286matrix_addition<MatA, MatB > operator+ (
287 matrix_expression<MatA, Device> const& A,
288 matrix_expression<MatB, Device> const& B
289){
290 REMORA_SIZE_CHECK(A().size1() == B().size1());
291 REMORA_SIZE_CHECK(A().size2() == B().size2());
292 return matrix_addition<MatA, MatB>(A(),B());
293}
294
295///\brief Subtracts two Matrices
296template<class MatA, class MatB, class Device>
297auto operator- (
298 matrix_expression<MatA, Device> const& A,
299 matrix_expression<MatB, Device> const& B
300) -> decltype(A() + (-B)){
301 REMORA_SIZE_CHECK(A().size1() == B().size1());
302 REMORA_SIZE_CHECK(A().size2() == B().size2());
303 return A() + (-B);
304}
305
306template<class MatA, class MatB, class Device>
307matrix_binary<MatA, MatB,
308 typename device_traits<Device>:: template safe_divide<typename common_value_type<MatA,MatB>::type>
309>
310safe_div(
311 matrix_expression<MatA, Device> const& A,
312 matrix_expression<MatB, Device> const& B,
313 typename common_value_type<MatA,MatB>::type defaultValue
314){
315 REMORA_SIZE_CHECK(A().size1() == B().size1());
316 REMORA_SIZE_CHECK(A().size2() == B().size2());
317 typedef typename common_value_type<MatA,MatB>::type result_type;
318 typedef typename device_traits<Device>:: template safe_divide<result_type> functor_type;
319 return matrix_binary<MatA, MatB, functor_type>(A(),B(), functor_type(defaultValue));
320}
321
322
323#define REMORA_BINARY_MATRIX_EXPRESSION(name, F)\
324template<class MatA, class MatB, class Device>\
325matrix_binary<MatA, MatB, typename device_traits<Device>:: template F<typename common_value_type<MatA,MatB>::type> >\
326name(matrix_expression<MatA, Device> const& m1, matrix_expression<MatB, Device> const& m2){\
327 REMORA_SIZE_CHECK(m1().size1() == m2().size1());\
328 REMORA_SIZE_CHECK(m1().size2() == m2().size2());\
329 typedef typename common_value_type<MatA,MatB>::type type;\
330 typedef typename device_traits<Device>:: template F<type> functor_type;\
331 return matrix_binary<MatA, MatB, functor_type >(m1(),m2(), functor_type());\
332}
333REMORA_BINARY_MATRIX_EXPRESSION(operator*, multiply)
334REMORA_BINARY_MATRIX_EXPRESSION(element_prod, multiply)
335REMORA_BINARY_MATRIX_EXPRESSION(operator/, divide)
336REMORA_BINARY_MATRIX_EXPRESSION(element_div, divide)
337REMORA_BINARY_MATRIX_EXPRESSION(pow,pow)
338REMORA_BINARY_MATRIX_EXPRESSION(min,min)
339REMORA_BINARY_MATRIX_EXPRESSION(max,max)
340#undef REMORA_BINARY_MATRIX_EXPRESSION
341
342/////////////////////////////////////////
343/////////Matrix-Products
344/////////////////////////////////////////
345
346
347///\brief Computes the outer product of two vectors.
348///
349/// The outer product of two vectors v1 and v2 is defined as the matrix
350/// (outer_prod (v1, v2))_ij [i] [j] = v1[i] * v2 [j]
351template<class VecA, class VecB, class Device>
352outer_product<VecA, VecB >
353outer_prod(
354 vector_expression<VecA, Device> const& v1,
355 vector_expression<VecB, Device> const& v2
356) {
357 return outer_product<VecA, VecB>(v1(), v2());
358}
359
360
361/// \brief computes the matrix-vector product x+=Av
362///
363/// The call to prod does not compute the product itself but instead, as all other expressions,
364/// it returns an expression-object which can compute it. In contrast to other expression,
365/// this expression is optimized to make use of well known mathematical identities to reduce run time of the algorithm.
366template<class MatA, class VecV, class Device>
367typename detail::matrix_vector_prod_optimizer<MatA,VecV>::type prod(
368 matrix_expression<MatA, Device> const& A,vector_expression<VecV, Device> const& v
369) {
370 REMORA_SIZE_CHECK(A().size2() == v().size());
371 return detail::matrix_vector_prod_optimizer<MatA,VecV>::create(A(),v());
372}
373
374/// \brief computes the matrix-vector product x+=v^TA
375///
376/// it is computed via the identity (v^TA)^T= A^Tv
377///
378/// The call to prod does not compute the product itself but instead, as all other expressions,
379/// it returns an expression-object which can compute it. In contrast to other expression,
380/// this expression is optimized to make use of well known mathematical identities to reduce run time of the algorithm.
381template<class MatA, class VecV, class Device>
382auto prod(vector_expression<VecV, Device> const& v,matrix_expression<MatA, Device> const& A) -> decltype(prod(trans(A),v)){
383 REMORA_SIZE_CHECK(A().size1() == v().size());
384 return prod(trans(A),v);
385}
386
387/// \brief Operator syntax for computes the matrix-vector product
388///
389/// v%A= prod(v,A).
390template<class MatA, class VecV, class Device>
391auto operator%(vector_expression<VecV, Device> const& v,matrix_expression<MatA, Device> const& A) -> decltype(prod(trans(A),v)){
392 REMORA_SIZE_CHECK(A().size1() == v().size());
393 return prod(trans(A),v);
394}
395
396/// \brief Operator syntax for computes the matrix-vector product
397///
398/// A%v = prod(A,v).
399template<class MatA, class VecV, class Device>
400auto operator%(matrix_expression<MatA, Device> const& A,vector_expression<VecV, Device> const& v) -> decltype(prod(A,v)){
401 REMORA_SIZE_CHECK(A().size2() == v().size());
402 return prod(A,v);
403}
404
405/// \brief Computes the matrix-vector product x+= alpha * Av or x= alpha * Av
406///
407/// A is interpreted as triangular matrix.
408/// The first template argument governs the type
409/// of triangular matrix: lower, upper, unit_lower and unit_upper.
410///
411///Example: x += triangular_prod<lower>(A,v);
412template<class TriangularType, class MatA, class VecV, class Device>
413auto triangular_prod(
414 matrix_expression<MatA, Device> const& A,
415 vector_expression<VecV, Device>& v
416) -> decltype(prod(to_triangular(A, TriangularType()), v)){
417 REMORA_SIZE_CHECK(A().size2() == v().size());
418 return prod(to_triangular(A, TriangularType()), v);
419}
420
421/// \brief computes the matrix-matrix product X+=AB
422template<class MatA, class MatB, class Device>
423typename detail::matrix_matrix_prod_optimizer<MatA,MatB>::type prod(
424 matrix_expression<MatA, Device> const& A,
425 matrix_expression<MatB, Device> const& B
426) {
427 REMORA_SIZE_CHECK(A().size2() == B().size1());
428 return detail::matrix_matrix_prod_optimizer<MatA,MatB>::create(A(),B());
429}
430
431/// \brief Operator syntax for computes the matrix-matrix product
432///
433/// A%B= prod(A,B).
434template<class MatA, class MatB, class Device>
435auto operator%(
436 matrix_expression<MatA, Device> const& A,
437 matrix_expression<MatB, Device> const& B
438) -> decltype(prod(A,B)){
439 REMORA_SIZE_CHECK(A().size2() == B().size1());
440 return prod(A,B);
441}
442
443/// \brief Computes the matrix-vector product x+= alpha * AB or x= alpha * AB
444///
445/// A is interpreted as triangular matrix.
446/// The first template argument governs the type
447/// of triangular matrix: lower, upper, unit_lower and unit_upper.
448/// B is interpreted as dense matrix.
449///
450///Example: x += triangular_prod<lower>(A,v);
451template<class TriangularType, class MatA, class MatB, class Device>
452auto triangular_prod(
453 matrix_expression<MatA, Device> const& A,
454 matrix_expression<MatB, Device> const& B
455) -> decltype(prod(to_triangular(A, TriangularType()), B)){
456 REMORA_SIZE_CHECK(A().size2() == B().size1());
457 return prod(to_triangular(A, TriangularType()), B);
458}
459
460
461//~ // inner_prod (set, v)_i = inner_prod(set_i,v)
462//~ template<class S, class V, class Device>
463//~ typename detail::vector_set_inner_prod_optimizer<S,V>::type
464//~ inner_prod(
465 //~ vector_set_expression<S, Device> const& set,
466 //~ vector_expression<V, Device> const& v
467//~ ){
468 //~ REMORA_SIZE_CHECK(set().point_size() == v().size());
469 //~ return detail::vector_set_inner_prod_optimizer<S,V>::create(set,v);
470//~ }
471
472//~ template<class S, class V, class Device>
473//~ typename detail::vector_set_inner_prod_optimizer<S,V>::type
474//~ inner_prod(
475 //~ vector_expression<V, Device> const& v
476 //~ vector_set_expression<S, Device> const& set
477//~ ){
478 //~ REMORA_SIZE_CHECK(set1().point_size() == v().size());
479 //~ return detail::vector_set_inner_prod_optimizer<S,V>::create(set,v);
480//~ }
481
482//~ template<class S, class M, class Device>
483//~ typename detail::vector_set_matrix_prod_optimizer<S,M>::type
484//~ operator%(
485 //~ vector_set_expression<S, Device> const& set,
486 //~ matrix_expression<M, Device> const& m
487//~ ){
488 //~ REMORA_SIZE_CHECK(set().point_size() == m().size1());
489 //~ return detail::vector_set_matrix_prod_optimizer<S,M>::create(set,m);
490//~ }
491
492//~ template<class S, class M, class Device>
493//~ auto operator%(
494 //~ matrix_expression<M, Device> const& m
495 //~ vector_set_expression<S, Device> const& set
496//~ )->decltype(set % trans(m)){
497 //~ REMORA_SIZE_CHECK(set().point_size() == m().size2());
498 //~ return set % trans(m);
499//~ }
500
501
502/////////////////////////////////////////
503//////////VECTOR-SET REDUCTIONS
504/////////////////////////////////////////
505
506/// \brief Computes the sum of elements for each point in the set
507///
508/// Formula: (sum S)_i = sum_j S_ij
509template<class S, class Device>
510typename detail::fold_vector_set_optimizer<
511 S, typename device_traits<Device>:: template add<typename S::value_type>
512 , typename device_traits<Device>:: template identity<typename S::value_type>
513>::type
514sum(vector_set_expression<S, Device> const& set) {
515 typedef typename device_traits<Device>:: template add<typename S::value_type> Add;
516 typedef typename device_traits<Device>:: template identity<typename S::value_type> Identity;
517 return detail::fold_vector_set_optimizer<S, Add, Identity>::create(set(), Add(), Identity());
518}
519
520/// \brief Computes the maximum of elements for each point in the set
521///
522/// Formula: (max S)_i = max_j S_ij
523template<class S, class Device>
524typename detail::fold_vector_set_optimizer<
525 S, typename device_traits<Device>:: template max<typename S::value_type>
526 , typename device_traits<Device>:: template identity<typename S::value_type>
527>::type
528max(vector_set_expression<S, Device> const& set) {
529 typedef typename device_traits<Device>:: template max<typename S::value_type> Max;
530 typedef typename device_traits<Device>:: template identity<typename S::value_type> Identity;
531 return detail::fold_vector_set_optimizer<S, Max, Identity>::create(set(), Max(), Identity());
532}
533
534/// \brief Computes the minimum of elements for each point in the set
535///
536/// Formula: (min S)_i = min_j S_ij
537template<class S, class Device>
538typename detail::fold_vector_set_optimizer<
539 S, typename device_traits<Device>:: template min<typename S::value_type>
540 , typename device_traits<Device>:: template identity<typename S::value_type>
541>::type
542min(vector_set_expression<S, Device> const& set) {
543 typedef typename device_traits<Device>:: template min<typename S::value_type> Min;
544 typedef typename device_traits<Device>:: template identity<typename S::value_type> Identity;
545 return detail::fold_vector_set_optimizer<S, Min, Identity>::create(set(), Min(), Identity());
546}
547
548//~ /// \brief arg_max v = arg max_i v_i
549//~ template<class M, class Device>
550//~ std::size_t arg_max(vector_set_expression<M, Device> const& set) {
551 //~ return kernels::vector_max(elem_result);
552//~ }
553
554//~ /// \brief arg_min v = arg min_i v_i
555//~ template<class VecV, class Device>
556//~ auto arg_min(vector_expression<VecV, Device> const& v) -> decltype(arg_max(-v)){
557 //~ return arg_max(-v);
558//~ }
559
560/// \brief Computes norm_1 for each element in the set.
561///
562/// Formula: norm_1(S)_i = norm_1(point(S,i))
563template<class S, class Device>
564auto norm_1(vector_set_expression<S, Device> const& set) ->decltype(sum(abs(set))) {
565 return sum(abs(set));
566}
567
568/// \brief Computes norm_sqr for each element in the set.
569///
570/// Formula: norm_sqr(S)_i = norm_sqr(point(S,i))
571template<class S, class Device>
572auto norm_sqr(vector_set_expression<S, Device> const& set) ->decltype(sum(sqr(set))){
573 return sum(sqr(set));
574}
575
576/// \brief Computes norm_2 for each element in the set.
577///
578/// Formula: norm_2(S)_i = norm_2(point(S,i))
579template<class S, class Device>
580auto norm_2(vector_set_expression<S, Device> const& set) ->decltype(sqrt(norm_sqr(set))){
581 return sqrt(norm_sqr(set));
582}
583
584/// \brief Computes norm_inf for each element in the set.
585///
586/// Formula: norm_inf(S)_i = norm_inf(point(S,i))
587template<class S, class Device>
588auto norm_inf(vector_set_expression<S, Device> const& set) ->decltype(max(abs(set))){
589 return max(abs(set));
590}
591
592
593/////////////////////////////////////////
594//////////MATRIX REDUCTIONS
595/////////////////////////////////////////
596
597/// \brief Computes the elementwise sum over all elements of A
598///
599/// returns a scalar s = sum_ij A_ij
600template<class MatA, class Device>
601typename MatA::value_type sum(matrix_expression<MatA, Device> const& A){
602 typedef typename MatA::value_type value_type;
603 typedef typename device_traits<Device>:: template add<value_type> functor_type;
604 auto const& elem_result = eval_block(A);
605 value_type result = 0;
606 kernels::matrix_fold<functor_type>(elem_result,result);
607 return result;
608}
609
610
611
612/// \brief Computes the elementwise maximum over all elements of A
613///
614/// returns a scalar s = max_ij A_ij
615template<class MatA, class Device>
616typename MatA::value_type max(matrix_expression<MatA, Device> const& A){
617 typedef typename MatA::value_type value_type;
618 typedef typename device_traits<Device>:: template max<value_type> functor_type;
619 auto const& elem_result = eval_block(A);
620 value_type result = 0;
621 kernels::matrix_fold<functor_type>(elem_result,result);
622 return result;
623}
624
625/// \brief Computes the elementwise minimum over all elements of A
626///
627/// returns a scalar s = min_ij A_ij
628template<class MatA, class Device>
629typename MatA::value_type min(matrix_expression<MatA, Device> const& A){
630 typedef typename MatA::value_type value_type;
631 typedef typename device_traits<Device>:: template min<value_type> functor_type;
632 auto const& elem_result = eval_block(A);
633 value_type result = 0;
634 kernels::matrix_fold<functor_type>(elem_result,result);
635 return result;
636}
637
638/// \brief Returns the frobenius inner-product between matrices exprssions 1 and B.
639///
640///The frobenius inner product is defined as \f$ <A,B>_F=\sum_{ij} A_ij*B_{ij} \f$. It induces the
641/// Frobenius norm \f$ ||A||_F = \sqrt{<A,A>_F} \f$
642template<class MatA, class MatB, class Device>
643decltype(typename MatA::value_type() * typename MatB::value_type())
644frobenius_prod(
645 matrix_expression<MatA, Device> const& A,
646 matrix_expression<MatB, Device> const& B
647){
648 REMORA_SIZE_CHECK(A().size1() == B().size1());
649 REMORA_SIZE_CHECK(A().size2() == B().size2());
650 return sum(eval_block(A*B));
651}
652
653/// \brief Computes the matrix 1-norm |A|_1
654///
655/// It is defined as \f$ \max_i \sum_j |A_{ij}| \f$
656template<class MatA, class Device>
657typename real_traits<typename MatA::value_type>::type
658norm_1(matrix_expression<MatA, Device> const& A) {
659 return max(norm_1(as_columns(A)));
660}
661
662/// \brief computes the frobenius norm |A|_F
663///
664/// It is defined as \f$ \sqrt{Tr(A^TA)}=\sqrt{\sum_{ij} A_{ij}^2} \f$
665template<class MatA, class Device>
666typename real_traits<typename MatA::value_type>::type
667norm_frobenius(matrix_expression<MatA, Device> const& A) {
668 using std::sqrt;
669 return sqrt(sum(sqr(eval_block(A))));
670}
671
672/// \brief Computes the matrix inf-norm |A|_inf
673///
674/// It is defined as \f$ \max_i \sum_j |A_{ij}| \f$
675template<class MatA, class Device>
676typename real_traits<typename MatA::value_type>::type
677norm_inf(matrix_expression<MatA, Device> const& A) {
678 return max(norm_1(as_rows(A)));
679}
680
681/// \brief Evaluates the trace of matrix A
682///
683/// The trace is defined as the sum of the diagonal elements of A,
684/// \f$ \text{trace}(A) = \sum_i A_{ii}\f$
685///
686/// \param A square matrix
687/// \return the sum of the values at the diagonal of \em A
688template < class MatA, class Device>
689typename MatA::value_type trace(matrix_expression<MatA, Device> const& A){
690 REMORA_SIZE_CHECK(A().size1() == A().size2());
691 return sum(diag(A));
692}
693
694/////////////////////////////////////////
695//////// Block-Matrix Creation
696/////////////////////////////////////////
697
698/// \brief Forms the block matrix (A|B) where B is to the right of A
699template<class MatA, class MatB, class Device>
700matrix_concat<MatA, MatB, true> operator|(
701 matrix_expression<MatA, Device> const& A,
702 matrix_expression<MatB, Device> const& B
703){
704 return matrix_concat<MatA, MatB, true>(A(),B());
705}
706
707/// \brief Forms the block matrix (A|v) where v is a column vector to the right of A
708template<class MatA, class VecV, class Device>
709auto operator|(
710 matrix_expression<MatA, Device> const& A,
711 vector_expression<VecV, Device> const& v
712) -> decltype(A | trans(repeat(v,1))){
713 return A | trans(repeat(v,1));
714}
715
716/// \brief Forms the block matrix (v|A) where v is a column vector to the left of A
717template<class MatA, class VecV, class Device>
718auto operator|(
719 vector_expression<VecV, Device> const& v,
720 matrix_expression<MatA, Device> const& A
721) -> decltype(trans(repeat(v,1)) | A){
722 return trans(repeat(v,1)) | A;
723}
724
725/// \brief Forms the block matrix (A|t)
726///
727/// The scalar t is interpreted as column vector
728template<class MatA, class T, class Device>
729typename std::enable_if<
730 std::is_convertible<T, typename MatA::value_type>::value,
731 matrix_concat<MatA, scalar_matrix<T, Device, row_major>, true >
732>::type operator|(
733 matrix_expression<MatA, Device> const& A,
734 T const& t
735){
736 return A | repeat(t,A().size1(), 1);
737}
738
739/// \brief Forms the block matrix (t|A)
740///
741/// The scalar t is interpreted as column vector
742template<class MatA, class T, class Device>
743typename std::enable_if<
744 std::is_convertible<T, typename MatA::value_type>::value,
745 matrix_concat<scalar_matrix<T, Device, row_major>, MatA, true >
746>::type operator|(
747 T const& t,
748 matrix_expression<MatA, Device> const& A
749){
750 return repeat(t,A().size1(), 1) | A;
751}
752
753///\brief Forms the block matrix A&B where A is on top of B
754template<class MatA, class MatB, class Device>
755matrix_concat<MatA, MatB, false> operator&(
756 matrix_expression<MatA, Device> const& A,
757 matrix_expression<MatB, Device> const& B
758){
759 return matrix_concat<MatA, MatB, false>(A(),B());
760}
761
762/// \brief Forms the block matrix (A & v) where v is a row vector on the bottom of A
763template<class MatA, class VecV, class Device>
764auto operator&(
765 matrix_expression<MatA, Device> const& A,
766 vector_expression<VecV, Device> const& v
767) -> decltype(A & repeat(v,1)){
768 return A & repeat(v,1);
769}
770
771/// \brief Forms the block matrix (A & v) where v is a row vector on the top of A
772template<class MatA, class VecV, class Device>
773auto operator&(
774 vector_expression<VecV, Device> const& v,
775 matrix_expression<MatA, Device> const& A
776) -> decltype(repeat(v,1) & A){
777 return repeat(v,1) & A;
778}
779
780/// \brief Forms the block matrix (A & t)
781///
782/// The scalar t is interpreted as row vector
783template<class MatA, class T, class Device>
784typename std::enable_if<
785 std::is_convertible<T, typename MatA::value_type>::value,
786 matrix_concat<MatA, scalar_matrix<T, Device, row_major>, false >
787>::type operator&(
788 matrix_expression<MatA, Device> const& A,
789 T const& t
790){
791 return A & repeat(t, 1, A().size2());
792}
793
794/// \brief Forms the block matrix (t & A)
795///
796/// The scalar t is interpreted as row vector
797template<class MatA, class T, class Device>
798typename std::enable_if<
799 std::is_convertible<T, typename MatA::value_type>::value,
800 matrix_concat<scalar_matrix<T, Device, row_major>, MatA, false >
801>::type operator&(
802 T const& t,
803 matrix_expression<MatA, Device> const& A
804){
805 return repeat(t,1, A().size2()) & A;
806}
807
808}
809
810#endif