solve.hpp
Go to the documentation of this file.
1/*!
2 * \brief Defines types for matrix decompositions
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_SOLVE_HPP
29#define REMORA_SOLVE_HPP
30
31#include "decompositions.hpp"
32
33namespace remora{
34
35
36//////////////////////////////////////////////////////
37////////Expression Types for Solving and inverse
38//////////////////////////////////////////////////////
39
40template<class MatA, class VecV, class SystemType, class Side>
41class matrix_vector_solve: public vector_expression<
42 matrix_vector_solve<MatA, VecV,SystemType, Side>,
43 typename MatA::device_type
44>{
45public:
46 typedef typename MatA::const_closure_type matrix_closure_type;
47 typedef typename VecV::const_closure_type vector_closure_type;
48 typedef decltype(
49 typename MatA::value_type() * typename VecV::value_type()
50 ) value_type;
51 typedef typename MatA::size_type size_type;
52 typedef value_type const_reference;
53 typedef const_reference reference;
54
55 typedef matrix_vector_solve const_closure_type;
56 typedef const_closure_type closure_type;
57 typedef unknown_storage storage_type;
58 typedef unknown_storage const_storage_type;
59 typedef blockwise<dense_tag> evaluation_category;
60 typedef typename MatA::device_type device_type;
61
62 size_type size() const {
63 return m_rhs.size();
64 }
65 typename device_traits<device_type>::queue_type& queue()const{
66 return m_matrix.queue();
67 }
68 matrix_vector_solve(
69 matrix_closure_type const& matrix, vector_closure_type const&rhs,
70 SystemType system_type = SystemType()
71 ):m_matrix(matrix), m_rhs(rhs), m_system_type(system_type){}
72
73 matrix_closure_type const& lhs()const{
74 return m_matrix;
75 }
76
77 vector_closure_type const& rhs()const{
78 return m_rhs;
79 }
80
81 SystemType const& system_type()const{
82 return m_system_type;
83 }
84
85 typedef no_iterator iterator;
86 typedef iterator const_iterator;
87
88 //dispatcher to computation kernels
89 template<class VecX>
90 void assign_to(vector_expression<VecX, device_type>& x)const{
91 assign(x,m_rhs);
92 solver<MatA,SystemType> alg(m_matrix, m_system_type);
93 alg.solve(x,Side());
94 }
95 template<class VecX>
96 void plus_assign_to(vector_expression<VecX, device_type>& x)const{
97 typename vector_temporary<VecX>::type temp(m_rhs);
98 solver<MatA,SystemType> alg(m_matrix, m_system_type);
99 alg.solve(temp,Side());
100 plus_assign(x,temp);
101 }
102private:
103 matrix_closure_type m_matrix;
104 vector_closure_type m_rhs;
105 SystemType m_system_type;
106};
107
108
109template<class MatA, class MatB, class SystemType, class Side>
110class matrix_matrix_solve: public matrix_expression<
111 matrix_matrix_solve<MatA, MatB, SystemType, Side>,
112 typename MatA::device_type
113>{
114public:
115 typedef typename MatA::const_closure_type matrixA_closure_type;
116 typedef typename MatB::const_closure_type matrixB_closure_type;
117 typedef decltype(
118 typename MatA::value_type() * typename MatB::value_type()
119 ) value_type;
120 typedef typename MatA::size_type size_type;
121 typedef value_type const_reference;
122 typedef const_reference reference;
123
124 typedef matrix_matrix_solve const_closure_type;
125 typedef const_closure_type closure_type;
126 typedef unknown_storage storage_type;
127 typedef unknown_storage const_storage_type;
128 typedef blockwise<dense_tag> evaluation_category;
129 typedef typename MatA::device_type device_type;
130 typedef unknown_orientation orientation;
131
132 matrix_matrix_solve(
133 matrixA_closure_type const& matrix, matrixB_closure_type const& rhs,
134 SystemType const& system_type =SystemType()
135 ):m_matrix(matrix), m_rhs(rhs), m_system_type(system_type){}
136
137 size_type size1() const {
138 return m_rhs.size1();
139 }
140 size_type size2() const {
141 return m_rhs.size2();
142 }
143 typename device_traits<device_type>::queue_type& queue()const{
144 return m_matrix.queue();
145 }
146
147 matrixA_closure_type const& lhs()const{
148 return m_matrix;
149 }
150
151 matrixB_closure_type const& rhs()const{
152 return m_rhs;
153 }
154
155 SystemType const& system_type()const{
156 return m_system_type;
157 }
158
159 typedef no_iterator major_iterator;
160 typedef major_iterator const_major_iterator;
161
162 //dispatcher to computation kernels
163 template<class MatX>
164 void assign_to(matrix_expression<MatX, device_type>& X)const{
165 assign(X,m_rhs);
166 solver<MatA,SystemType> alg(m_matrix,m_system_type);
167 alg.solve(X,Side());
168 }
169 template<class MatX>
170 void plus_assign_to(matrix_expression<MatX, device_type>& X)const{
171 typename matrix_temporary<MatX>::type temp(m_rhs);
172 solver<MatA,SystemType> alg(m_matrix,m_system_type);
173 alg.solve(temp,Side());
174 plus_assign(X,temp);
175 }
176private:
177 matrixA_closure_type m_matrix;
178 matrixB_closure_type m_rhs;
179 SystemType m_system_type;
180};
181
182template<class MatA, class SystemType>
183class matrix_inverse: public matrix_expression<
184 matrix_inverse<MatA, SystemType>,
185 typename MatA::device_type
186>{
187public:
188 typedef typename MatA::const_closure_type matrix_closure_type;
189 typedef typename MatA::value_type value_type;
190 typedef typename MatA::size_type size_type;
191 typedef value_type const_reference;
192 typedef const_reference reference;
193
194 typedef matrix_inverse const_closure_type;
195 typedef const_closure_type closure_type;
196 typedef unknown_storage storage_type;
197 typedef unknown_storage const_storage_type;
198 typedef blockwise<dense_tag> evaluation_category;
199 typedef typename MatA::device_type device_type;
200 typedef typename MatA::orientation orientation;
201
202 matrix_inverse(matrix_closure_type const& matrix, SystemType system_type)
203 :m_matrix(matrix), m_system_type(system_type){}
204
205 size_type size1() const {
206 return m_matrix.size1();
207 }
208 size_type size2() const {
209 return m_matrix.size1();
210 }
211
212 typename device_traits<device_type>::queue_type& queue()const{
213 return m_matrix.queue();
214 }
215
216 matrix_closure_type const& matrix()const{
217 return m_matrix;
218 }
219
220 SystemType const& system_type()const{
221 return m_system_type;
222 }
223
224 typedef no_iterator major_iterator;
225 typedef major_iterator const_major_iterator;
226
227 //dispatcher to computation kernels
228 template<class MatX>
229 void assign_to(matrix_expression<MatX, device_type>& X)const{
230 typedef scalar_vector<value_type, device_type> diag_vec;
231 assign(X,diagonal_matrix<diag_vec>(diag_vec(size1(),value_type(1))));
232 solver<MatA,SystemType> alg(m_matrix,m_system_type);
233 alg.solve(X,left());
234 }
235 template<class MatX>
236 void plus_assign_to(matrix_expression<MatX, device_type>& X)const{
237 typedef scalar_vector<value_type, device_type> diag_vec;
238 typename matrix_temporary<MatX>::type temp = diagonal_matrix<diag_vec>(diag_vec(size1(),value_type(1)));
239 solver<MatA,SystemType> alg(m_matrix,m_system_type);
240 alg.solve(temp,left());
241 plus_assign(X,temp);
242 }
243private:
244 matrix_closure_type m_matrix;
245 SystemType m_system_type;
246};
247
248//////////////////////////////////////////////////////
249////////Expression Optimizations
250//////////////////////////////////////////////////////
251
252namespace detail{
253
254////////////////////////////////////
255//// Vector Solve
256////////////////////////////////////
257template<class M, class V, class Tag, class Side>
258struct matrix_vector_solve_optimizer{
259 typedef matrix_vector_solve<M,V, Tag, Side> type;
260
261 static type create(
262 typename M::const_closure_type const& m,
263 typename V::const_closure_type const& v,
264 Tag t
265 ){
266 return type(m,v,t);
267 }
268};
269
270////////////////////////////////////
271//// Matrix Solve
272////////////////////////////////////
273template<class M1, class M2, class Tag, class Side>
274struct matrix_matrix_solve_optimizer{
275 typedef matrix_matrix_solve<M1,M2, Tag, Side> type;
276
277 static type create(
278 typename M1::const_closure_type const& lhs,
279 typename M2::const_closure_type const& rhs,
280 Tag t
281 ){
282 return type(lhs,rhs,t);
283 }
284};
285
286////////////////////////////////////
287//// Matrix Inverse
288////////////////////////////////////
289template<class M, class Tag>
290struct matrix_inverse_optimizer{
291 typedef matrix_inverse<M, Tag> type;
292
293 static type create(typename M::const_closure_type const& m, Tag t){
294 return type(m,t);
295 }
296};
297
298//////////////////////////////////
299/////Interactions with other expressions
300//////////////////////////////////
301
302//small helper needed for transpose
303template<class T>
304struct solve_tag_transpose_helper{
305 static T transpose(T t){return t;}
306};
307template<bool Upper, bool Unit>
308struct solve_tag_transpose_helper<triangular_tag<Upper,Unit> >{
309 static triangular_tag<!Upper,Unit> transpose(triangular_tag<Upper,Unit>){return triangular_tag<!Upper,Unit>();}
310};
311
312//trans(solve(A,B,left)) = solve(trans(A),trans(B),right)
313//trans(solve(A,B,right)) = solve(trans(A),trans(B),left)
314template<class M1, class M2, bool Left, class Tag>
315struct matrix_transpose_optimizer<matrix_matrix_solve<M1,M2, Tag, system_tag<Left> > >{
316 typedef matrix_transpose_optimizer<typename M2::const_closure_type> lhs_opt;
317 typedef matrix_transpose_optimizer<typename M2::const_closure_type> rhs_opt;
318 typedef matrix_matrix_solve_optimizer<
319 typename lhs_opt::type,typename rhs_opt::type,
320 typename Tag::transposed_tag, system_tag<!Left>
321 > opt;
322 typedef typename opt::type type;
323
324 static type create(matrix_matrix_solve<M1,M2, Tag, system_tag<Left> > const& m){
325 return opt::create(
326 lhs_opt::create(m.rhs()),rhs_opt::create(m.lhs()),
327 solve_tag_transpose_helper<Tag>::transpose(m.system_type())
328 );
329 }
330};
331
332
333template<class M, class Tag>
334struct matrix_transpose_optimizer<matrix_inverse<M, Tag> >{
335 typedef matrix_transpose_optimizer<typename M::const_closure_type> mat_opt;
336 typedef matrix_inverse_optimizer<
337 typename mat_opt::type, typename Tag::transposed_orientation
338 > opt;
339 typedef typename opt::type type;
340
341 static type create(matrix_inverse<M, Tag> const& m){
342 return opt::create(
343 mat_opt::create(m.matrix()),
344 solve_tag_transpose_helper<Tag>::transpose(m.system_type())
345 );
346 }
347};
348
349
350
351//prod(inv(A),b) = solve(A,b,left)
352template<class M, class V, class Tag>
353struct matrix_vector_prod_optimizer<matrix_inverse<M,Tag>, V>{
354 typedef matrix_vector_solve_optimizer<M,V, Tag, left> opt;
355 typedef typename opt::type type;
356
357 static type create(matrix_inverse<M,Tag> const& inv, typename V::const_closure_type const& v){
358 return opt::create(inv.matrix(),v,inv.system_type());
359 }
360};
361
362//prod(solve(A,B,left),c) = solve(A, prod(B,c),right)
363template<class M1, class M2,class V, class Tag>
364struct matrix_vector_prod_optimizer<matrix_matrix_solve<M1,M2, Tag, left >, V >{
365 typedef matrix_vector_prod_optimizer<M2,V> prod_opt;
366 typedef matrix_vector_solve_optimizer<M1, typename prod_opt::type, Tag, left> opt;
367 typedef typename opt::type type;
368
369 static type create(matrix_matrix_solve<M1,M2, Tag, left> const& m, typename V::const_closure_type const& v){
370 return opt::create(
371 m.lhs(), prod_opt::create(m.rhs(),v), m.system_type()
372 );
373 }
374};
375
376//prod(solve(A,B,right),c) = prod(B,solve(A,c, left))
377template<class M1, class M2, class V, class Tag>
378struct matrix_vector_prod_optimizer<matrix_matrix_solve<M1,M2, Tag, right >, V >{
379 typedef matrix_vector_solve_optimizer<M1, V, Tag, left> solve_opt;
380 typedef matrix_vector_prod_optimizer<M2,typename solve_opt::type> opt;
381 typedef typename opt::type type;
382 static type create(matrix_matrix_solve<M1,M2, Tag, right> const& m, typename V::const_closure_type const& v){
383 return opt::create(
384 m.rhs(), solve_opt::create(m.lhs(),v,m.system_type())
385 );
386 }
387};
388
389//row(solve(A,B,left),i) = prod(solve(A,e_i,right),B) = prod(trans(B),solve(A,e_i,right))
390template<class M1, class M2,class Tag>
391struct matrix_row_optimizer<matrix_matrix_solve<M1,M2, Tag, left > >{
392 typedef matrix_transpose_optimizer<typename M2::const_closure_type> rhs_opt;
393 typedef unit_vector<typename M1::value_type, typename M1::device_type> unit;
394 typedef matrix_vector_solve_optimizer<M1, unit, Tag, right> solve_opt;
395 typedef matrix_vector_prod_optimizer<typename rhs_opt::type,typename solve_opt::type> opt;
396 typedef typename opt::type type;
397
398 static type create(matrix_matrix_solve<M1,M2, Tag, left> const& m, std::size_t i){
399 return opt::create(
400 rhs_opt::create(m.rhs()),
401 solve_opt::create(m.lhs(),unit(m.lhs().size2(),i), m.system_type())
402 );
403 }
404};
405
406//row(solve(A,B,right),i) = solve(A,row(B,i),right)
407template<class M1, class M2, class Tag>
408struct matrix_row_optimizer<matrix_matrix_solve<M1,M2, Tag, right > >{
409 typedef matrix_row_optimizer<typename M2::const_closure_type> rhs_opt;
410 typedef matrix_vector_solve_optimizer<M1, typename rhs_opt::type, Tag, right> opt;
411 typedef typename opt::type type;
412
413 static type create(matrix_matrix_solve<M1,M2, Tag, right> const& m, std::size_t i){
414 return opt::create(m.lhs(),rhs_opt::create(m.rhs(),i), m.system_type());
415 }
416};
417
418//row(inv(A),i) = solve(A,e_1,right)
419template<class M, class Tag>
420struct matrix_row_optimizer<matrix_inverse<M, Tag > >{
421 typedef unit_vector<typename M::value_type, typename M::device_type> unit;
422 typedef matrix_vector_solve_optimizer<M, unit, Tag, right> opt;
423 typedef typename opt::type type;
424
425 static type create(matrix_inverse<M, Tag > const& m, std::size_t i){
426 return opt::create(m.matrix(), unit(m.lhs().size2(),i), m.system_type());
427 }
428};
429
430//prod(inv(A),B) = solve(A,B,left)
431template<class M1, class M2, class Tag>
432struct matrix_matrix_prod_optimizer<matrix_inverse<M1,Tag>, M2>{
433 typedef matrix_matrix_solve_optimizer<M1,M2, Tag, left> opt;
434 typedef typename opt::type type;
435
436 static type create(matrix_inverse<M1,Tag> const& inv, typename M2::const_closure_type const& m){
437 return opt::create(inv.matrix(),m,inv.system_type());
438 }
439};
440
441//prod(B,inv(A)) = solve(A,B,right)
442template<class M1, class M2, class Tag>
443struct matrix_matrix_prod_optimizer<M1, matrix_inverse<M2,Tag> >{
444 typedef matrix_matrix_solve_optimizer<M2,M1, Tag, right> opt;
445 typedef typename opt::type type;
446
447 static type create(typename M1::const_closure_type const& m, matrix_inverse<M2,Tag> const& inv){
448 return opt::create(inv.matrix(),m,inv.system_type());
449 }
450};
451
452
453}
454
455
456//solvers for vector rhs
457template<class MatA, class VecB, bool Left, class Device, class SystemType>
458typename detail::matrix_vector_solve_optimizer<MatA,VecB, SystemType, system_tag<Left> >::type
459solve(
460 matrix_expression<MatA, Device> const& A,
461 vector_expression<VecB, Device> const& b,
462 SystemType t, system_tag<Left>
463){
464 REMORA_SIZE_CHECK(A().size1() == A().size2());
465 REMORA_SIZE_CHECK(A().size1() == b().size());
466 typedef detail::matrix_vector_solve_optimizer<MatA,VecB, SystemType, system_tag<Left> > opt;
467 return opt::create(A(),b(), t);
468}
469//solvers for matrix rhs
470template<class MatA, class MatB, bool Left, class Device, class SystemType>
471typename detail::matrix_matrix_solve_optimizer<MatA,MatB, SystemType, system_tag<Left> >::type
472solve(
473 matrix_expression<MatA, Device> const& A,
474 matrix_expression<MatB, Device> const& B,
475 SystemType t, system_tag<Left>
476){
477 REMORA_SIZE_CHECK(A().size1() == A().size2());
478 typedef detail::matrix_matrix_solve_optimizer<MatA,MatB, SystemType, system_tag<Left> > opt;
479 return opt::create(A(),B(), t);
480}
481
482//matrix inverses
483template<class MatA, class Device, class SystemType>
484typename detail::matrix_inverse_optimizer<MatA, SystemType>::type
485inv(
486 matrix_expression<MatA, Device> const& A, SystemType t
487){
488 REMORA_SIZE_CHECK(A().size1() == A().size2());
489 typedef detail::matrix_inverse_optimizer<MatA,SystemType> opt;
490 return opt::create(A(), t);
491}
492
493
494}
495#endif