30 #ifndef _OPERATORLINEAR_HPP_
31 #define _OPERATORLINEAR_HPP_
42 template<
class DomainSpace,
class DualImageSpace>
52 typedef typename super::domain_space_type domain_space_type;
53 typedef typename super::dual_image_space_type dual_image_space_type;
54 typedef typename super::domain_space_ptrtype domain_space_ptrtype;
55 typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
56 typedef typename domain_space_type::element_type domain_element_type;
57 typedef typename dual_image_space_type::element_type dual_image_element_type;
60 typedef typename super::backend_ptrtype backend_ptrtype;
63 typedef typename backend_type::vector_ptrtype vector_ptrtype;
64 typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
66 template<
typename T,
typename Storage>
67 struct domain_element:
public super::domain_space_type::template Element<T,Storage> {};
69 typedef typename domain_space_type::template Element<
typename domain_space_type::value_type,
72 typedef typename domain_space_type::template Element<
typename domain_space_type::value_type,
75 typedef typename dual_image_space_type::template Element<
typename dual_image_space_type::value_type,
78 typedef typename dual_image_space_type::template Element<
typename dual_image_space_type::value_type,
84 typedef FsFunctionalLinear<DualImageSpace> image_element_type;
91 M_pattern( Pattern::COUPLED ),
92 M_name(
"operatorlinear")
97 M_backend( ol.M_backend ),
99 M_pattern( ol.M_pattern ),
109 M_matrix = ol.M_matrix;
113 dual_image_space_ptrtype dualImageSpace,
114 backend_ptrtype backend,
115 bool buildMatrix =
true ,
119 M_backend( backend ),
120 M_pattern( pattern ),
121 M_name(
"operatorlinear" )
123 if ( buildMatrix ) M_matrix = M_backend->newMatrix( _trial=domainSpace, _test=dualImageSpace , _pattern=M_pattern );
129 init( domain_space_ptrtype domainSpace,
130 dual_image_space_ptrtype dualImageSpace,
131 backend_ptrtype backend,
132 bool buildMatrix =
true ,
135 this->setDomainSpace( domainSpace );
136 this->setDualImageSpace( dualImageSpace );
139 if ( buildMatrix ) M_matrix = M_backend->newMatrix( _trial=domainSpace, _test=dualImageSpace , _pattern=M_pattern );
142 void setName( std::string name ) { M_name = name; }
143 std::string name()
const {
return M_name ; }
147 template<
typename Storage>
149 apply(
const domain_element<typename domain_element_type::value_type, Storage>& de,
150 image_element_type& ie )
const
152 if ( ! M_matrix->closed() )
157 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
158 *_v1 = de;_v1->close();
159 vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
160 M_backend->prod( M_matrix, _v1, _v2 );
161 ie.container() = *_v2;
165 apply(
const domain_element_type& de,
166 image_element_type& ie )
const
168 if ( ! M_matrix->closed() )
173 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
174 *_v1 = de;_v1->close();
175 vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
176 M_backend->prod( M_matrix, _v1, _v2 );
177 ie.container() = *_v2;
181 energy(
const typename domain_space_type::element_type & de,
182 const typename dual_image_space_type::element_type & ie )
const
184 if ( ! M_matrix->closed() )
189 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
190 *_v1 = de;_v1->close();
191 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
193 vector_ptrtype _v3( M_backend->newVector( _test=ie.functionSpace() ) );
194 M_backend->prod( M_matrix, _v1, _v3 );
200 apply(
const typename domain_space_type::element_type & de,
201 typename dual_image_space_type::element_type & ie )
203 if ( ! M_matrix->closed() )
208 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
209 *_v1 = de;_v1->close();
210 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
211 M_backend->prod( M_matrix, _v1, _v2 );
212 ie.container() = *_v2;
219 apply(
const domain_element_range_type & de,
220 typename dual_image_space_type::element_type & ie )
222 if ( ! M_matrix->closed() )
227 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
228 *_v1 = de;_v1->close();
229 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
230 M_backend->prod( M_matrix, _v1, _v2 );
231 ie.container() = *_v2;
235 apply(
const typename domain_space_type::element_type & de,
236 dual_image_element_range_type & ie )
238 if ( ! M_matrix->closed() )
243 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
244 *_v1 = de;_v1->close();
245 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
246 M_backend->prod( M_matrix, _v1, _v2 );
247 ie.container() = *_v2;
252 apply(
const domain_element_range_type & de,
253 dual_image_element_range_type & ie )
255 if ( ! M_matrix->closed() )
260 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
261 *_v1 = de;_v1->close();
262 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
263 M_backend->prod( M_matrix, _v1, _v2 );
264 ie.container() = *_v2;
269 apply(
const domain_element_slice_type & de,
270 typename dual_image_space_type::element_type & ie )
272 if ( ! M_matrix->closed() )
277 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
278 *_v1 = de;_v1->close();
279 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
280 M_backend->prod( M_matrix, _v1, _v2 );
281 ie.container() = *_v2;
285 apply(
const typename domain_space_type::element_type & de,
286 dual_image_element_slice_type & ie )
288 if ( ! M_matrix->closed() )
293 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
294 *_v1 = de;_v1->close();
295 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
296 M_backend->prod( M_matrix, _v1, _v2 );
297 ie.container() = *_v2;
302 apply( domain_element_slice_type de,
303 dual_image_element_slice_type ie )
305 if ( ! M_matrix->closed() )
310 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
311 *_v1 = de;_v1->close();
312 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
313 M_backend->prod( M_matrix, _v1, _v2 );
314 ie.container() = *_v2;
318 apply(
const typename domain_space_type::element_type::component_type & de,
319 typename dual_image_space_type::element_type & ie )
321 if ( ! M_matrix->closed() )
326 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
327 *_v1 = de;_v1->close();
328 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
329 M_backend->prod( M_matrix, _v1, _v2 );
330 ie.container() = *_v2;
338 apply(
const domain_element_range_type & de,
339 dual_image_element_slice_type & ie )
341 if ( ! M_matrix->closed() )
346 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
347 *_v1 = de;_v1->close();
348 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
349 M_backend->prod( M_matrix, _v1, _v2 );
350 ie.container() = *_v2;
354 apply(
const domain_element_slice_type & de,
355 dual_image_element_range_type & ie )
357 if ( ! M_matrix->closed() )
362 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
363 *_v1 = de;_v1->close();
364 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
365 M_backend->prod( M_matrix, _v1, _v2 );
366 ie.container() = *_v2;
373 template<
typename Storage1,
typename Storage2>
375 apply(
const domain_element<typename domain_element_type::value_type, Storage1>& de,
376 dual_image_element<typename dual_image_element_type::value_type, Storage2>& ie )
378 if ( ! M_matrix->closed() )
383 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
384 *_v1 = de;_v1->close();
386 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
387 M_backend->prod( M_matrix, _v1, _v2 );
388 ie.container() = *_v2;
391 template <
typename T1 =
typename domain_space_type::element_type,
392 typename T2 =
typename dual_image_space_type::element_type >
394 operator()( T1 & de )
396 T2 elt_image( this->dualImageSpace(),
"oio" );
397 this->apply( de,elt_image );
402 template <
typename T1 =
typename domain_space_type::element_type,
403 typename T2 =
typename dual_image_space_type::element_type >
405 operator()( boost::shared_ptr<T1> de )
407 T2 elt_image( this->dualImageSpace(),
"oio" );
408 this->apply( *de,elt_image );
417 const image_element_type& ie )
419 if ( ! M_matrix->closed() )
424 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
425 vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
426 *_v2 = ie.container();
427 M_backend->solve( M_matrix, M_matrix, _v1, _v2 );
432 if ( !M_backend->converged() )
434 std::cerr <<
"[OperatorLinear::applyInverse] "
435 <<
"solver failed to converge" << std::endl;
445 template<
typename RhsExpr>
447 inv( RhsExpr
const& rhs_expr )
449 if ( ! M_matrix->closed() )
454 domain_element_type de = this->domainSpace()->element();
456 auto ie = M_backend->newVector( _test=this->dualImageSpace() );
457 form1( _test=this->dualImageSpace(), _vector=ie ) =
459 rhs_expr *
id( this->dualImageSpace()->element() ) );
461 M_backend->solve( M_matrix, de, ie );
468 template<
class ExprT>
469 this_type& operator=( ExprT
const& e )
472 form2( _trial=this->domainSpace(),
473 _test=this->dualImageSpace(),
474 _matrix=M_matrix ) = e;
478 this_type& operator=( this_type
const& m )
480 M_backend = m.M_backend;
482 M_matrix->addMatrix( 1.0, m.M_matrix );
483 M_pattern = m.M_pattern;
489 template<
class ExprT>
490 this_type& operator+=( ExprT
const& e )
492 form2( _trial=this->domainSpace(),
493 _test=this->dualImageSpace(),
494 _matrix=M_matrix ) += e;
500 if ( ! M_matrix->closed() )
513 matrix_type
const& mat()
const
519 matrix_ptrtype
const& matPtr()
const
525 matrix_ptrtype& matPtr()
530 virtual void matPtr( matrix_ptrtype & matrix )
536 OperatorLinear& add( T
const& scalar, OperatorLinear
const& ol )
539 M_matrix->addMatrix( scalar, *ol.M_matrix );
543 backend_ptrtype& backend()
554 OperatorLinear& add( T
const& scalar, boost::shared_ptr<OperatorLinear> ol )
557 this->add( scalar, *ol );
561 OperatorLinear& operator+( boost::shared_ptr<OperatorLinear> ol )
564 this->add( 1.0, *ol );
567 OperatorLinear& operator+( OperatorLinear
const& ol )
570 this->add( 1.0, ol );
576 backend_ptrtype M_backend;
577 matrix_ptrtype M_matrix;
584 template<
typename Args>
585 struct compute_opLinear_return
587 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::domainSpace>::type>::type::element_type domain_space_type;
588 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::imageSpace>::type>::type::element_type image_space_type;
590 typedef boost::shared_ptr<OperatorLinear<domain_space_type, image_space_type> > type;
593 BOOST_PARAMETER_FUNCTION(
594 (
typename compute_opLinear_return<Args>::type ),
598 ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
599 ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
602 ( backend, *, Backend<
typename compute_opLinear_return<Args>::domain_space_type::value_type>::build() )
603 ( pattern, *, (
size_type)Pattern::COUPLED )
607 Feel::detail::ignore_unused_variable_warning( args );
608 typedef OperatorLinear<typename compute_opLinear_return<Args>::domain_space_type,
609 typename compute_opLinear_return<Args>::image_space_type> operatorlinear_type;
611 boost::shared_ptr<operatorlinear_type> opI(
new operatorlinear_type( domainSpace,imageSpace,backend,pattern ) );
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
Linear Operator between function spaces, represented by a matrix.
Definition: operatorlinear.hpp:43
elements_type const & elements() const
Definition: elements.hpp:355
Definition: solverlinear.hpp:33
base class for all linear algebra backends
Definition: backend.hpp:133
static backend_ptrtype build(BackendType=BACKEND_GMM, WorldComm const &worldComm=Environment::worldComm())
Definition: backend.cpp:167
type_traits< T >::real_type inner_product(Vector< T > const &v1, Vector< T > const &v2)
Definition: vector.hpp:579
interface to vector
Definition: matrixeigendense.hpp:50
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Definition: matrixsparse.hpp:76
Operator between function spaces.
Definition: operator.hpp:44
virtual void applyInverse(domain_element_type &de, const image_element_type &ie)
apply the inverse of the operator:
Definition: operatorlinear.hpp:416