Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
operatortrilinos.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@epfl.ch>
6  Date: 2006-03-06
7 
8  Copyright (C) 2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
30 #ifndef __operator_trilinos_matrix_H
31 #define __operator_trilinos_matrix_H 1
32 
33 
36 // #include <feel/feelalg/operatortrilinos.hpp>
37 
39 
40 #if defined( FEELPP_HAS_TRILINOS_EPETRA )
41 #include "Epetra_ConfigDefs.h"
42 #include "Epetra_Operator.h"
43 
44 
45 #if defined( FEELPP_HAS_MPI )
46 #include "mpi.h"
47 #include "Epetra_MpiComm.h"
48 #else
49 #include "Epetra_SerialComm.h"
50 #endif
51 #include "Epetra_Map.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_MultiVector.h"
54 #include "Epetra_Operator.h"
55 
56 
57 
58 
59 namespace Feel
60 {
61 class OperatorMatrix : public Epetra_Operator
62 {
63 public:
64 
65  typedef MatrixSparse<double> sparse_matrix_type;
66  typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
67 
68  typedef MatrixEpetra epetra_sparse_matrix_type;
69  typedef boost::shared_ptr<epetra_sparse_matrix_type> epetra_sparse_matrix_ptrtype;
70 
71  typedef Vector<double> vector_type;
72  typedef boost::shared_ptr<vector_type> vector_ptrtype;
73 
74  typedef VectorEpetra<double> epetra_vector_type;
75  typedef boost::shared_ptr<epetra_vector_type> epetra_vector_ptrtype;
76 
77  typedef Epetra_Operator prec_type;
78  typedef boost::shared_ptr<prec_type> prec_ptrtype;
79 
80  typedef SolverLinearTrilinos<double> solver_type;
81  typedef boost::shared_ptr<solver_type> solver_ptrtype;
82 
83  typedef Teuchos::ParameterList list_type;
84 
85  typedef solver_type::real_type real_type;
86 
87 
88  OperatorMatrix( sparse_matrix_ptrtype const& F, std::string _label, bool transpose = 0 )
89  :
90  M_F( F ),
91  M_Matrix( dynamic_cast<epetra_sparse_matrix_type const*>( F.get() ) ),
92  M_Prec( ),
93  M_Solver( new solver_type() ),
94  M_hasInverse( 1 ),
95  M_hasApply( 1 ),
96  M_useTranspose( transpose ),
97  M_domainMap( M_Matrix->getDomainMap() ),
98  M_rangeMap( M_Matrix->getRangeMap() ),
99  M_label( _label ),
100  M_maxiter( 1000 ),
101  M_tol( 1e-10 )
102  {
103  DVLOG(2) << "Create operator " << Label() << " ...\n";
104  }
105 
106  OperatorMatrix( sparse_matrix_ptrtype const& F,
107  list_type const& options,
108  std::string _label,
109  prec_ptrtype Prec, bool transpose = 0 )
110  :
111  M_F( F ),
112  M_Matrix( dynamic_cast<epetra_sparse_matrix_type const*>( F.get() ) ),
113  M_Prec( Prec ),
114  M_Solver( new solver_type() ),
115  M_hasInverse( 1 ),
116  M_hasApply( 1 ),
117  M_useTranspose( transpose ),
118  M_domainMap( M_Matrix->getDomainMap() ),
119  M_rangeMap( M_Matrix->getRangeMap() ),
120  M_label( _label ),
121  M_maxiter( 1000 ),
122  M_tol( 1e-10 )
123  {
124  DVLOG(2) << "Create operator " << Label() << " ...\n";
125 
126  M_Solver->setOptions( options );
127 
128  M_tol = M_Solver->getOptions().get( "tol", 1e-10 );
129  M_maxiter = M_Solver->getOptions().get( "max_iter", 100 );
130 
131 
132  }
133 
134  OperatorMatrix( const OperatorMatrix& tc )
135  :
136  M_F( tc.M_F ),
137  M_Matrix( tc.M_Matrix ),
138  M_Prec( tc.M_Prec ),
139  M_Solver( tc.M_Solver ),
140  M_hasInverse( tc.M_hasInverse ),
141  M_hasApply( tc.M_hasApply ),
142  M_useTranspose( tc.M_useTranspose ),
143  M_domainMap( tc.M_domainMap ),
144  M_rangeMap( tc.M_rangeMap ),
145  M_label( tc.M_label ),
146  M_maxiter( tc.M_maxiter ),
147  M_tol( tc.M_tol )
148  {
149  DVLOG(2) << "Copy operator " << Label() << " ...\n";
150  }
151 
152  bool hasInverse() const
153  {
154  return M_hasInverse;
155  }
156 
157  bool hasApply() const
158  {
159  return M_hasApply;
160  }
161 
162 
163  int Apply( const vector_ptrtype& X, vector_ptrtype& Y ) const
164  {
165  Apply( dynamic_cast<epetra_vector_type const&>( *X ).vec(),
166  dynamic_cast<epetra_vector_type&>( *Y ).vec() );
167 
168  return !hasApply();
169  }
170 
171  int Apply( const Epetra_MultiVector & X, Epetra_MultiVector & Y ) const
172  {
173  DVLOG(2) << "Apply Operator " << Label() << "\n";
174 
175  M_Matrix->multiply( false,X,Y );
176 
177  DVLOG(2) << "Apply Operator " << Label() << " successfully\n";
178  return !hasApply();
179  }
180 
181  int ApplyInverse ( const vector_ptrtype& X, vector_ptrtype& Y ) const
182  {
183  FEELPP_ASSERT( hasInverse() ).error( "This operator cannot be inverted." );
184 
185  int result = ApplyInverse( dynamic_cast<epetra_vector_type const&>( *X ).vec(),
186  dynamic_cast<epetra_vector_type&>( *Y ).vec() );
187 
188  return result;
189  }
190 
191  int ApplyInverse ( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
192  {
193  FEELPP_ASSERT( hasInverse() ).error( "This operator cannot be inverted." );
194 
195  std::pair<unsigned int, real_type> result;
196 
197  DVLOG(2) << "Apply Inverse Operator " << Label() << "\n";
198 
199  if ( M_Prec.get() != 0 )
200  result = M_Solver->solve( *M_Matrix, M_Prec, Y, X, M_tol, M_maxiter );
201 
202  else
203  result = M_Solver->solve( *M_Matrix, Y, X, M_tol, M_maxiter );
204 
205  DVLOG(2) << "Apply Inverse Operator " << Label() << " successfully\n";
206 
207  int result_integer = result.first;
208 
209  return result_integer;
210  }
211 
212  // other function
213  int SetUseTranspose( bool UseTranspose = 0 )
214  {
215  M_useTranspose = UseTranspose;
216 
217  return UseTranspose;
218  }
219 
220  double NormInf() const
221  {
222  return M_Matrix->linftyNorm();
223  }
224 
225  const char * Label () const
226  {
227  return( M_label.c_str() );
228  }
229 
230  bool UseTranspose() const
231  {
232  return M_useTranspose;
233  }
234 
235  bool HasNormInf () const
236  {
237  return( true );
238  }
239 
240  const Epetra_Comm& Comm() const
241  {
242  return( M_Matrix->getDomainMap().Comm() );
243  }
244 
245  const Epetra_Map& OperatorDomainMap() const
246  {
247  return( M_domainMap );
248  }
249 
250  const Epetra_Map& OperatorRangeMap() const
251  {
252  return( M_rangeMap );
253  }
254 
255  ~OperatorMatrix()
256  {
257  DVLOG(2) << "Destroyed matrix operator: " << this->Label() << " ...\n";
258  };
259 
260 private:
261 
262  sparse_matrix_ptrtype M_F;
263  epetra_sparse_matrix_type const* M_Matrix;
264 
265  prec_ptrtype M_Prec;
266 
267  solver_ptrtype M_Solver;
268 
269  bool M_hasInverse, M_hasApply, M_useTranspose;
270 
271  Epetra_Map M_domainMap, M_rangeMap;
272 
273  std::string M_label;
274 
275  int M_maxiter;
276 
277  double M_tol;
278 };
279 
280 
281 
282 
283 
284 
285 
286 
287 
288 template< typename operator_type >
289 class OperatorInverse : public Epetra_Operator
290 {
291 public:
292 
293  typedef boost::shared_ptr<operator_type> operator_ptrtype;
294 
295  // This constructor implements the F^-1 operator
296  OperatorInverse( operator_ptrtype& F )
297  :
298  M_F( F )
299  {
300  this->setName();
301  DVLOG(2) << "Create inverse operator " << this->Label() << "...\n";
302  }
303 
304  OperatorInverse( const OperatorInverse& tc )
305  :
306  M_F( tc.M_F ),
307  M_Label( tc.M_Label )
308  {
309  DVLOG(2) << "Copy inverse operator " << this->Label() << "...\n";
310  }
311 
312  bool hasInverse() const
313  {
314  return M_F->hasApply();
315  }
316 
317  bool hasApply() const
318  {
319  return M_F->hasInverse();
320  }
321 
322  int Apply( const Epetra_MultiVector & X, Epetra_MultiVector & Y ) const
323  {
324  DVLOG(2) << "Apply matrix " << Label() << "\n";
325 
326  FEELPP_ASSERT( hasApply() ).error( "This operator cannot be applied." );
327  M_F->ApplyInverse( X,Y );
328 
329  return !hasApply();
330  }
331 
332  int ApplyInverse ( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
333  {
334  DVLOG(2) << "ApplyInverse matrix " << Label() << "\n";
335 
336  FEELPP_ASSERT( hasInverse() ).error( "This operator cannot be inverted." );
337 
338  M_F->Apply( X,Y );
339 
340  return !hasInverse();
341  }
342 
343  // other function
344  int SetUseTranspose( bool UseTranspose )
345  {
346  return( false );
347  }
348 
349  double NormInf() const
350  {
351  return( false );
352  }
353 
354  const char * Label () const
355  {
356  return M_Label.c_str();
357  }
358 
359  bool UseTranspose() const
360  {
361  return( false );
362  }
363 
364  bool HasNormInf () const
365  {
366  return( false );
367  }
368 
369 
370  const Epetra_Comm & Comm() const
371  {
372  return( M_F->Comm() );
373  }
374 
375  const Epetra_Map & OperatorDomainMap() const
376  {
377  return( M_F->OperatorRangeMap() );
378  }
379 
380  const Epetra_Map & OperatorRangeMap() const
381  {
382  return( M_F->OperatorDomainMap() );
383  }
384 
385  ~OperatorInverse()
386  {
387  DVLOG(2) << "Destroyed inverse operator: " << this->Label() << " ...\n";
388  };
389 
390 private:
391 
392  operator_ptrtype M_F;
393 
394  std::string M_Label;
395 
396  void setName()
397  {
398  std::string L = M_F->Label();
399  L.append( ")" );
400  std::string temp( "inv(" );
401  temp.append( L );
402 
403  M_Label = temp;
404  }
405 };
406 
407 
408 
409 
410 
411 
412 template< typename op1_type, typename op2_type >
413 class OperatorCompose : public Epetra_Operator
414 {
415 public:
416 
417  typedef boost::shared_ptr<op1_type> op1_ptrtype;
418  typedef boost::shared_ptr<op2_type> op2_ptrtype;
419 
420  // This constructor implements the (F o G) operator
421 
422  OperatorCompose()
423  :
424  M_F(),
425  M_G()
426  {
427  std::string t( M_F->Label() );
428  std::string u( M_G->Label() );
429 
430  t.append( "*" );
431  t.append( u );
432  M_Label = t;
433  }
434 
435  OperatorCompose( op1_ptrtype& F, op2_ptrtype& G )
436  :
437  M_F( F ),
438  M_G( G )
439  {
440  std::string t( F->Label() );
441  std::string u( G->Label() );
442 
443  t.append( "*" );
444  t.append( u );
445  M_Label = t;
446  DVLOG(2) << "Create operator " << Label() << " ...\n";
447  }
448 
449  OperatorCompose( const OperatorCompose& tc )
450  :
451  M_F( tc.M_F ),
452  M_G( tc.M_G ),
453  M_Label( tc.M_Label )
454  {
455  DVLOG(2) << "Copy operator " << Label() << " ...\n";
456  }
457 
458  bool hasInverse() const
459  {
460  return M_F->hasInverse() * M_G->hasInverse();
461  }
462 
463  bool hasApply() const
464  {
465  return M_F->hasApply() * M_G->hasApply();
466  }
467 
468  int Apply( const Epetra_MultiVector & X, Epetra_MultiVector & Y ) const
469  {
470  FEELPP_ASSERT( hasApply() ).error( "This operator cannot be applied." );
471 
472  DVLOG(2) << "Apply operator " << Label() << " ...\n";
473 
474  Epetra_MultiVector Z( M_G->OperatorRangeMap(), 1 );
475 
476  M_G->Apply( X,Z );
477  M_F->Apply( Z,Y );
478 
479  return !hasApply();
480  }
481 
482  int ApplyInverse ( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
483  {
484  FEELPP_ASSERT( hasInverse() ).error( "This operator cannot be inverted." );
485 
486  DVLOG(2) << "Apply Inverse operator " << Label() << " ...\n";
487 
488  Epetra_MultiVector Z( X );
489 
490  M_F->ApplyInverse( X,Z );
491  M_G->ApplyInverse( Z,Y );
492 
493  return hasInverse();
494  }
495 
496  // other function
497  int SetUseTranspose( bool UseTranspose )
498  {
499  return( false );
500  }
501 
502  double NormInf() const
503  {
504  return 0;
505  }
506 
507  const char * Label () const
508  {
509 
510  return M_Label.c_str();
511  }
512 
513  bool UseTranspose() const
514  {
515  return( false );
516  }
517 
518  bool HasNormInf () const
519  {
520  return( false );
521  }
522 
523 
524  const Epetra_Comm & Comm() const
525  {
526  return( M_F->Comm() );
527  }
528 
529  const Epetra_Map & OperatorDomainMap() const
530  {
531  return( M_G->OperatorDomainMap() );
532  }
533 
534  const Epetra_Map & OperatorRangeMap() const
535  {
536  return( M_F->OperatorRangeMap() );
537  }
538 
539  ~OperatorCompose()
540  {
541  DVLOG(2) << "Destroyed compose operator: " << this->Label() << " ...\n";
542  };
543 
544 private:
545 
546  op1_ptrtype M_F;
547  op2_ptrtype M_G;
548 
549  std::string M_Label;
550 };
551 
552 
553 
554 
555 
556 
557 
558 
559 
560 
561 
562 
563 
564 template< typename op1_type>
565 class OperatorScale : public Epetra_Operator
566 {
567 public:
568 
569  typedef boost::shared_ptr<op1_type> op1_ptrtype;
570 
571  // This constructor implements the (\alpha F) operator
572  OperatorScale()
573  :
574  M_F(),
575  M_alpha( 0 )
576  {
577  }
578 
579  OperatorScale( op1_ptrtype& F )
580  :
581  M_F( F ),
582  M_alpha( 1 )
583  {
584  std::string temp( "alpha" );
585  std::string t = M_F->Label();
586 
587  temp.append( "." );
588  temp.append( t );
589 
590  M_Label = temp;
591 
592  DVLOG(2) << "Create scale operator " << Label() << " ...\n";
593  }
594 
595  OperatorScale( op1_ptrtype& F, double alpha )
596  :
597  M_F( F ),
598  M_alpha( alpha )
599  {
600  std::string temp( "alpha" );
601  std::string t = M_F->Label();
602 
603  temp.append( "." );
604  temp.append( t );
605 
606  M_Label = temp;
607 
608  DVLOG(2) << "Create scale operator " << Label() << " ...\n";
609  }
610 
611  OperatorScale( const OperatorScale& tc )
612  :
613  M_F( tc.M_F ),
614  M_alpha( tc.M_alpha ),
615  M_Label( tc.M_Label )
616  {
617  DVLOG(2) << "Copy scale operator " << Label() << " ...\n";
618  }
619 
620  bool hasInverse() const
621  {
622  return M_F->hasApply()*( M_alpha != 0 );
623  }
624 
625  bool hasApply() const
626  {
627  return M_F->hasApply();
628  }
629 
630 
631 
632  int Apply( const Epetra_MultiVector & X, Epetra_MultiVector & Y ) const
633  {
634  DVLOG(2) << "Apply scale operator " << Label() << "\n";
635 
636  M_F->Apply( X,Y );
637 
638  Y.Scale( M_alpha );
639 
640  return !hasApply();
641  }
642 
643  int ApplyInverse ( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
644  {
645  DVLOG(2) << "ApplyInverse scale operator " << Label() << "\n";
646 
647  FEELPP_ASSERT( hasInverse() && ( M_alpha != 0 ) ).error( "This operator cannot be inverted." );
648 
649  Epetra_MultiVector Z( X );
650  Z.Scale( 1./M_alpha );
651 
652  M_F->ApplyInverse( Z,Y );
653 
654  return !hasInverse();
655  }
656 
657  // other function
658  int SetUseTranspose( bool UseTranspose )
659  {
660  return( false );
661  }
662 
663  double NormInf() const
664  {
665  return 0;
666  }
667 
668  const char * Label () const
669  {
670  return M_Label.c_str();
671  }
672 
673  bool UseTranspose() const
674  {
675  return( false );
676  }
677 
678  bool HasNormInf () const
679  {
680  return( false );
681  }
682 
683 
684  const Epetra_Comm & Comm() const
685  {
686  return( M_F->Comm() );
687  }
688 
689  const Epetra_Map & OperatorDomainMap() const
690  {
691  return( M_F->OperatorDomainMap() );
692  }
693 
694  const Epetra_Map & OperatorRangeMap() const
695  {
696  return( M_F->OperatorRangeMap() );
697  }
698 
699  ~OperatorScale()
700  {
701  //DVLOG(2) << "Destroyed scale operator: " << Label() << " ...\n";
702  };
703 
704 private:
705 
706  op1_ptrtype M_F;
707 
708  double M_alpha;
709 
710  std::string M_Label;
711 };
712 
713 
714 
715 
716 
717 
718 template< typename operator_type >
719 class OperatorFree : public Epetra_Operator
720 {
721 public:
722 
723  typedef boost::shared_ptr<operator_type> operator_ptrtype;
724 
725  typedef Epetra_Operator prec_type;
726  typedef boost::shared_ptr<prec_type> prec_ptrtype;
727 
728  typedef SolverLinearTrilinos<double> solver_type;
729  typedef boost::shared_ptr<solver_type> solver_ptrtype;
730 
731  typedef solver_type::real_type real_type;
732 
733  typedef Teuchos::ParameterList list_type;
734 
735  OperatorFree( operator_ptrtype F )
736  :
737  M_op( F ),
738  M_Solver( new solver_type() ),
739  M_hasInverse( 0 ),
740  M_hasApply( F->hasApply() ),
741  M_useTranspose( F->UseTranspose() ),
742  M_label( F->Label() )
743  {
744  DVLOG(2) << "Create operator " << Label() << " ...\n";
745  }
746 
747  OperatorFree( operator_ptrtype F, prec_ptrtype Prec, list_type options )
748  :
749  M_op( F ),
750  M_Solver( new solver_type() ),
751  M_hasInverse( 1 ),
752  M_hasApply( F->hasApply() ),
753  M_useTranspose( F->UseTranspose() ),
754  M_label( F->Label() ),
755  M_Prec( Prec )
756  {
757  DVLOG(2) << "Create operator " << Label() << " ...\n";
758  M_tol = options.get( "tol", 1e-7 );
759  M_maxiter = options.get( "max_iter", 100 );
760 
761  M_Solver->setOptions( options );
762  }
763 
764 
765  OperatorFree( const OperatorFree& tc )
766  :
767  M_op( tc.M_op ),
768  M_hasInverse( tc.M_hasInverse ),
769  M_hasApply( tc.M_hasApply ),
770  M_useTranspose( tc.M_useTranspose ),
771  M_Solver( tc.M_Solver ),
772  M_label( tc.M_label ),
773  M_Prec( tc.M_Prec )
774  {
775  DVLOG(2) << "Copy operator " << Label() << " ...\n";
776  }
777 
778 
779  bool hasInverse() const
780  {
781  return M_hasInverse;
782  }
783 
784  bool hasApply() const
785  {
786  return M_hasApply;
787  }
788 
789 
790  int Apply( const Epetra_MultiVector & X, Epetra_MultiVector & Y ) const
791  {
792  DVLOG(2) << "Apply operator " << Label() << "\n";
793  M_op->Apply( X,Y );
794  DVLOG(2) << "Finished Apply operator " << Label() << "\n";
795 
796  return !hasApply();
797  }
798 
799  int ApplyInverse ( const Epetra_MultiVector& X, Epetra_MultiVector& Y ) const
800  {
801  DVLOG(2) << "ApplyInverse operator " << Label() << "\n";
802 
803  FEELPP_ASSERT( hasInverse() ).error( "This operator cannot be inverted." );
804 
805  std::pair<unsigned int, real_type> result = M_Solver->solve( M_op, M_Prec, Y, X, M_tol, M_maxiter );
806 
807  DVLOG(2) << "Finished ApplyInverse operator " << Label() << "\n";
808  return !hasInverse();
809  }
810 
811  // other function
812  int SetUseTranspose( bool UseTranspose = 0 )
813  {
814  M_useTranspose = UseTranspose;
815 
816  return UseTranspose;
817  }
818 
819  double NormInf() const
820  {
821  return 0;
822  }
823 
824  const char * Label () const
825  {
826  return( M_label.c_str() );
827  }
828 
829  bool UseTranspose() const
830  {
831  return M_useTranspose;
832  }
833 
834  bool HasNormInf () const
835  {
836  return( true );
837  }
838 
839  const Epetra_Comm & Comm() const
840  {
841  return( M_op->OperatorDomainMap().Comm() );
842  }
843 
844  const Epetra_Map & OperatorDomainMap() const
845  {
846  return( M_op->OperatorDomainMap() );
847  }
848 
849  const Epetra_Map & OperatorRangeMap() const
850  {
851  return( M_op->OperatorRangeMap() );
852  }
853 
854  ~OperatorFree()
855  {
856  DVLOG(2) << "Destroyed operator: " << Label() << " ...\n";
857  };
858 
859 private:
860 
861  operator_ptrtype M_op;
862 
863  solver_ptrtype M_Solver;
864 
865  bool M_hasInverse, M_hasApply, M_useTranspose;
866 
867  int M_maxiter;
868  double M_tol;
869 
870  std::string M_label;
871 
872  prec_ptrtype M_Prec;
873 };
874 
875 
876 
877 
878 
879 
880 
881 
882 } // Feel
883 
884 #endif // FEELPP_HAS_TRILINOS_EPETRA
885 #endif /* __operator_trilinos_matrix_H */

Generated on Sun Oct 20 2013 08:25:03 for Feel++ by doxygen 1.8.4