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
matrixpetsc.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@feelpp.org>
6  Date: 2005-10-18
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2007 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __MatrixPetsc_H
31 #define __MatrixPetsc_H 1
32 
33 #include <feel/feelconfig.h>
34 
35 #if defined(FEELPP_HAS_PETSC_H)
36 
37 
39 
42 
43 
44 extern "C"
45 {
46 #include <petscmat.h>
47 }
48 
49 #ifndef USE_COMPLEX_NUMBERS
50 extern "C" {
51 # include <petscversion.h>
52 # if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
53 # include <petscsles.h>
54 # else
55 # include <petscksp.h>
56 # endif
57 }
58 #else
59 # include <petscversion.h>
60 # if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
61 # include <petscsles.h>
62 # else
63 # include <petscksp.h>
64 # endif
65 #endif
66 
67 
68 
69 namespace Feel
70 {
71 template<typename T> class VectorPetsc;
72 
85 template<typename T>
86 class MatrixPetsc : public MatrixSparse<T>
87 {
88  typedef MatrixSparse<T> super;
89 public:
93 
94  static const bool is_row_major = true;
95 
97 
101 
102  typedef typename super::value_type value_type;
103  typedef typename super::real_type real_type;
104 
105  typedef std::vector<std::set<size_type> > pattern_type;
106 
107  typedef typename super::graph_type graph_type;
108  typedef typename super::graph_ptrtype graph_ptrtype;
109 
110  typedef typename super::datamap_type datamap_type;
111  typedef typename super::datamap_ptrtype datamap_ptrtype;
112 
114 
118 
134  MatrixPetsc();
135 
136  MatrixPetsc( datamap_ptrtype const& dmRow, datamap_ptrtype const& dmCol, WorldComm const& worldComm=Environment::worldComm() );
137 
138 
146  MatrixPetsc ( Mat m );
147  MatrixPetsc ( Mat m, datamap_ptrtype const& dmRow, datamap_ptrtype const& dmCol );
148  MatrixPetsc ( MatrixSparse<value_type> const& M, IS& isrow, IS& iscol );
149  MatrixPetsc ( MatrixSparse<value_type> const& M, std::vector<int> const& rowIndex, std::vector<int> const& colIndex );
155  ~MatrixPetsc();
156 
158 
162 
174  value_type operator () ( const size_type i,
175  const size_type j ) const;
176 
181  MatrixPetsc& operator=( MatrixSparse<value_type> const& M );
182 
184 
188 
193  size_type size1 () const;
194 
199  size_type size2 () const;
200 
205  size_type rowStart () const;
206 
211  size_type rowStop () const;
212 
214 
218 
219 
221 
225 
234  void init ( const size_type m,
235  const size_type n,
236  const size_type m_l,
237  const size_type n_l,
238  const size_type nnz=30,
239  const size_type noz=10 );
240 
244  void init ( const size_type m,
245  const size_type n,
246  const size_type m_l,
247  const size_type n_l,
248  graph_ptrtype const& graph );
249 
253  void setIndexSplit( std::vector< std::vector<size_type> > const &indexSplit );
254 
258  void resize( size_type m, size_type n, bool /*preserve*/ )
259  {
260  this->setInitialized( false );
261  init( m, n, m, n );
262  }
263 
264  void fill( pattern_type const& patt )
265  {}
272  void clear ();
273 
278  void zero ();
279 
285  void zero ( size_type start1, size_type stop1, size_type start2, size_type stop2 );
286 
292  void close () const;
293 
298  bool closed() const;
299 
307  real_type l1Norm () const;
308 
319  real_type linftyNorm () const;
320 
327  void set ( const size_type i,
328  const size_type j,
329  const value_type& value );
330 
339  void add ( const size_type i,
340  const size_type j,
341  const value_type& value );
342 
349  void addMatrix ( const ublas::matrix<value_type> &dm,
350  const std::vector<size_type> &rows,
351  const std::vector<size_type> &cols );
352 
359  void addMatrix ( int* rows, int nrows,
360  int* cols, int ncols,
361  value_type* data );
362 
367  void addMatrix ( const ublas::matrix<value_type> &dm,
368  const std::vector<size_type> &dof_indices )
369  {
370  this->addMatrix ( dm, dof_indices, dof_indices );
371  }
372 
384  void addMatrix ( const T a, MatrixSparse<T> &X );
385 
391  void matMatMult ( MatrixSparse<T> const& In, MatrixSparse<T> &Res );
392 
397  void scale( T const a );
398 
402  void diagonal ( Vector<value_type>& dest ) const;
403 
410  void transpose( MatrixSparse<value_type>& Mt ) const;
411 
415  virtual void symmetricPart( MatrixSparse<value_type>& Ms ) const;
416 
422  Mat mat () const
423  {
424  FEELPP_ASSERT ( M_mat != NULL ).error( "null petsc matrix" );
425  return M_mat;
426  }
427  Mat& mat ()
428  {
429  FEELPP_ASSERT ( M_mat != NULL ).warn( "null petsc matrix" );
430  return M_mat;
431  }
432 
439  void printMatlab( const std::string name="NULL" ) const;
440 
444  value_type
445  energy( Vector<value_type> const& __v,
446  Vector<value_type> const& __u,
447  bool transpose = false ) const;
448 
456  void zeroRows( std::vector<int> const& rows, Vector<value_type> const& values, Vector<value_type>& rhs, Context const& on_context );
457 
461  void updateBlockMat( boost::shared_ptr<MatrixSparse<T> > m, std::vector<size_type> start_i, std::vector<size_type> start_j );
462 
463 
464  void updatePCFieldSplit( PC & pc );
465 
466  std::vector<IS> const& petscSplitIS() const { return M_petscIS; }
467  std::map<PC*,bool > & mapSplitPC() { return M_mapPC; }
468 
469  std::vector<PetscInt> ia() { return M_ia; }
470  std::vector<PetscInt> ja() { return M_ja; }
471 
472 
474 
475  /*
476  * Set zero entries diagonal if missing : only for PETSC!
477  */
478  void zeroEntriesDiagonal();
479 
480 private:
481 
482  // disable
483  MatrixPetsc( MatrixPetsc const & );
484 
485 protected:
486 
490  Mat M_mat;
491 
492 private:
493 
494 
495  std::vector<IS> M_petscIS;
496 
497  std::map<PC*,bool > M_mapPC;
498 
503  const bool M_destroy_mat_on_exit;
504  std::vector<PetscInt> M_ia,M_ja;
505 };
506 
507 
508 
509 template<typename T>
510 class MatrixPetscMPI : public MatrixPetsc<T>
511 {
512  typedef MatrixPetsc<T> super;
513 
514 public :
515 
516  typedef typename super::graph_type graph_type;
517  typedef typename super::graph_ptrtype graph_ptrtype;
518  typedef typename super::value_type value_type;
519 
520  typedef typename super::datamap_type datamap_type;
521  typedef typename super::datamap_ptrtype datamap_ptrtype;
522 
523  MatrixPetscMPI();
524 
525  MatrixPetscMPI( datamap_ptrtype const& dmRow, datamap_ptrtype const& dmCol, WorldComm const& worldComm=Environment::worldComm() );
526 
527  MatrixPetscMPI( Mat m, datamap_ptrtype const& dmRow, datamap_ptrtype const& dmCol );
528 
529  ~MatrixPetscMPI()
530  {
531  this->clear();
532  }
533 
534  void init ( const size_type m,
535  const size_type n,
536  const size_type m_l,
537  const size_type n_l,
538  const size_type nnz=30,
539  const size_type noz=10 );
540 
544  void init ( const size_type m,
545  const size_type n,
546  const size_type m_l,
547  const size_type n_l,
548  graph_ptrtype const& graph );
549 
550 
551  size_type size1() const;
552  size_type size2() const;
553  size_type rowStart() const;
554  size_type rowStop() const;
555  size_type colStart() const;
556  size_type colStop() const;
557 
558  void set( const size_type i,
559  const size_type j,
560  const value_type& value );
561 
562  void add ( const size_type i,
563  const size_type j,
564  const value_type& value );
565 
566  void addMatrix( const ublas::matrix<value_type>& dm,
567  const std::vector<size_type>& rows,
568  const std::vector<size_type>& cols );
569 
570  void addMatrix( int* rows, int nrows,
571  int* cols, int ncols,
572  value_type* data );
573 
574  void addMatrix( const T a, MatrixSparse<T> &X );
575 
576 
577  void zero();
578  void zero( size_type start1, size_type stop1, size_type start2, size_type stop2 );
579  //void zeroEntriesDiagonal();
580  void zeroRows( std::vector<int> const& rows,
581  Vector<value_type> const& values,
582  Vector<value_type>& rhs,
583  Context const& on_context );
584 
585  value_type energy( Vector<value_type> const& __v,
586  Vector<value_type> const& __u,
587  bool transpose = false ) const;
588 
589 private :
590 
591  void addMatrixSameNonZeroPattern( const T a, MatrixSparse<T> &X );
592 
593 
594 };
595 
596 
597 
598 } // Feel
599 #endif /* FEELPP_HAS_PETSC */
600 #endif /* __MatrixPetsc_H */

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