1#ifndef DUNE_FEM_PETSCLINEAROPERATOR_HH
2#define DUNE_FEM_PETSCLINEAROPERATOR_HH
13#include <dune/common/dynmatrix.hh>
44 struct PetscSolverParameter :
public LocalParameter< SolverParameter, PetscSolverParameter >
46 typedef LocalParameter< SolverParameter, PetscSolverParameter > BaseType;
49 using BaseType :: parameter ;
50 using BaseType :: keyPrefix ;
52 PetscSolverParameter(
const ParameterReader ¶meter = Parameter::container() )
53 : BaseType( parameter )
56 PetscSolverParameter(
const SolverParameter& sp )
57 : PetscSolverParameter( sp.keyPrefix(), sp.parameter() )
60 PetscSolverParameter(
const std::string &keyPrefix,
const ParameterReader ¶meter = Parameter::container() )
61 : BaseType( keyPrefix, parameter )
64 bool isPetscSolverParameter()
const {
return true; }
66 static const int boomeramg = 0;
67 static const int parasails = 1;
68 static const int pilut = 2;
70 int hypreMethod()
const
72 const std::string hyprePCNames[] = {
"boomer-amg",
"parasails",
"pilu-t" };
74 if (parameter().exists(
"petsc.preconditioning.method"))
76 hypreType = parameter().getEnum(
"petsc.preconditioning.hypre.method", hyprePCNames, 0 );
77 std::cout <<
"WARNING: using deprecated parameter 'petsc.preconditioning.hypre.method' use "
78 << keyPrefix() <<
"preconditioning.hypre.method instead\n";
81 hypreType = parameter().getEnum( keyPrefix()+
"hypre.method", hyprePCNames, 0 );
85 int superluMethod()
const
87 const std::string factorizationNames[] = {
"petsc",
"superlu",
"mumps" };
89 if (parameter().exists(
"petsc.preconditioning.lu.method"))
91 factorType = parameter().getEnum(
"petsc.preconditioning.lu.method", factorizationNames, 0 );
92 std::cout <<
"WARNING: using deprecated parameter 'petsc.preconditioning.lu.method' use "
93 << keyPrefix() <<
"preconditioning.lu.method instead\n";
96 factorType = parameter().getEnum( keyPrefix()+
"preconditioning.lu.method", factorizationNames, 0 );
100 bool viennaCL ()
const {
101 return parameter().getValue<
bool > ( keyPrefix() +
"petsc.viennacl", false );
103 bool blockedMode ()
const {
104 return parameter().getValue<
bool > ( keyPrefix() +
"petsc.blockedmode", true );
113 template<
typename DomainFunction,
typename RangeFunction >
114 class PetscLinearOperator
115 :
public Fem::AssembledOperator< DomainFunction, RangeFunction >
117 typedef PetscLinearOperator< DomainFunction, RangeFunction > ThisType;
119 typedef Mat MatrixType;
120 typedef DomainFunction DomainFunctionType;
121 typedef RangeFunction RangeFunctionType;
122 typedef typename DomainFunctionType::RangeFieldType DomainFieldType;
123 typedef typename RangeFunctionType::RangeFieldType RangeFieldType;
124 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
125 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
127 typedef PetscDiscreteFunction< DomainSpaceType > PetscDomainFunctionType;
128 typedef PetscDiscreteFunction< RangeSpaceType > PetscRangeFunctionType;
130 typedef typename DomainSpaceType::GridPartType::template Codim< 0 >::EntityType DomainEntityType;
131 typedef typename RangeSpaceType::GridPartType::template Codim< 0 >::EntityType RangeEntityType;
133 static const unsigned int domainLocalBlockSize = DomainSpaceType::localBlockSize;
134 static const unsigned int rangeLocalBlockSize = RangeSpaceType::localBlockSize;
136 static constexpr bool squareBlocked = domainLocalBlockSize == rangeLocalBlockSize ;
137 static constexpr bool blockedMatrix = domainLocalBlockSize > 1 && squareBlocked;
141 typedef FlatFieldMatrix< RangeFieldType, domainLocalBlockSize, rangeLocalBlockSize > MatrixBlockType;
142 typedef MatrixBlockType block_type;
145 enum Status {statAssembled=0,statAdd=1,statInsert=2,statGet=3,statNothing=4};
147 typedef PetscMappers< DomainSpaceType > DomainMappersType;
148 typedef PetscMappers< RangeSpaceType > RangeMappersType;
150 typedef AuxiliaryDofs< typename DomainSpaceType::GridPartType, typename DomainMappersType::GhostMapperType > DomainAuxiliaryDofsType;
151 typedef AuxiliaryDofs< typename RangeSpaceType::GridPartType, typename RangeMappersType::GhostMapperType > RangeAuxiliaryDofsType;
157 struct LocalMatrixTraits
159 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
160 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
161 typedef LocalMatrix LocalMatrixType;
162 typedef typename RangeSpaceType::RangeFieldType RangeFieldType;
165 typedef RangeFieldType LittleBlockType;
169 typedef LocalMatrix ObjectType;
170 typedef ThisType LocalMatrixFactoryType;
174 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
175 typedef ColumnObject< ThisType > LocalColumnObjectType;
177 PetscLinearOperator (
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
178 const PetscSolverParameter& param = PetscSolverParameter() )
179 : domainMappers_( domainSpace ),
180 rangeMappers_( rangeSpace ),
181 localMatrixStack_( *this ),
182 status_(statNothing),
183 viennaCL_( param.viennaCL() ),
184 blockedMode_( blockedMatrix && (!viennaCL_) && param.blockedMode() )
187 PetscLinearOperator (
const std::string &,
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
188 const PetscSolverParameter& param = PetscSolverParameter() )
189 : PetscLinearOperator( domainSpace, rangeSpace, param )
193 ~PetscLinearOperator ()
200 ::Dune::Petsc::MatAssemblyBegin( petscMatrix_, MAT_FLUSH_ASSEMBLY );
201 ::Dune::Petsc::MatAssemblyEnd ( petscMatrix_, MAT_FLUSH_ASSEMBLY );
203 status_ = statAssembled;
208 if( ! finalizedAlready() )
210 ::Dune::Petsc::MatAssemblyBegin( petscMatrix_, MAT_FINAL_ASSEMBLY );
211 ::Dune::Petsc::MatAssemblyEnd ( petscMatrix_, MAT_FINAL_ASSEMBLY );
213 if( ! unitRows_.empty() )
215 ::Dune::Petsc::MatZeroRows( petscMatrix_, unitRows_.size(), unitRows_.data(), 1. );
223 bool finalizedAlready()
const
225 PetscBool assembled = PETSC_FALSE ;
226 ::Dune::Petsc::MatAssembled( petscMatrix_, &assembled );
227 return assembled == PETSC_TRUE;
230 void finalizeAssembly ()
const
232 const_cast< ThisType&
> (*this).finalize();
236 const DomainSpaceType &domainSpace ()
const {
return domainMappers_.space(); }
237 const RangeSpaceType &rangeSpace ()
const {
return rangeMappers_.space(); }
242 template <
class DF,
class RF>
243 void apply (
const DF &arg, RF &dest )
const
246 petscArg_.reset(
new PetscDomainFunctionType(
"PetscOp-arg", domainSpace() ) );
248 petscDest_.reset(
new PetscRangeFunctionType(
"PetscOp-arg", rangeSpace() ) );
250 petscArg_->assign( arg );
251 apply( *petscArg_, *petscDest_ );
252 dest.assign( *petscDest_ );
256 void apply (
const PetscDomainFunctionType &arg, PetscRangeFunctionType &dest )
const
260 ::Dune::Petsc::MatMult( petscMatrix_, *arg.petscVec() , *dest.petscVec() );
263 void operator() (
const DomainFunctionType &arg, RangeFunctionType &dest )
const
270 reserve( SimpleStencil<DomainSpaceType,RangeSpaceType>(0),
true );
274 void reserve (
const std::vector< Set >& sparsityPattern )
276 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ),
true );
280 template <
class StencilType>
281 void reserve (
const StencilType &stencil,
const bool isSimpleStencil =
false )
283 domainMappers_.update();
284 rangeMappers_.update();
286 if(sequence_ != domainSpace().sequence())
298 ::Dune::Petsc::MatCreate( domainSpace().gridPart().comm(), &petscMatrix_ );
301 const PetscInt localCols = domainMappers_.ghostMapper().interiorSize() * domainLocalBlockSize;
302 const PetscInt localRows = rangeMappers_.ghostMapper().interiorSize() * rangeLocalBlockSize;
304 const PetscInt globalCols = domainMappers_.parallelMapper().size() * domainLocalBlockSize;
305 const PetscInt globalRows = rangeMappers_.parallelMapper().size() * rangeLocalBlockSize;
307 ::Dune::Petsc::MatSetSizes( petscMatrix_, localRows, localCols, globalRows, globalCols );
312 ::Dune::Petsc::MatSetType( petscMatrix_, MATAIJVIENNACL );
314 else if( blockedMode_ )
316 bs = domainLocalBlockSize ;
317 ::Dune::Petsc::MatSetType( petscMatrix_, MATBAIJ );
319 ::Dune::Petsc::MatSetBlockSize( petscMatrix_, bs );
323 ::Dune::Petsc::MatSetType( petscMatrix_, MATAIJ );
341 if( isSimpleStencil || std::is_same< StencilType,SimpleStencil<DomainSpaceType,RangeSpaceType> >::value )
343 ::Dune::Petsc::MatSetUp( petscMatrix_, bs, domainLocalBlockSize * stencil.maxNonZerosEstimate() );
347 DomainAuxiliaryDofsType domainAuxiliaryDofs( domainMappers_.ghostMapper() );
348 RangeAuxiliaryDofsType rangeAuxiliaryDofs( rangeMappers_.ghostMapper() );
350 std::vector< PetscInt > d_nnz( localRows / bs, 0 );
351 std::vector< PetscInt > o_nnz( localRows / bs, 0 );
352 for(
const auto& entry : stencil.globalStencil() )
354 const int petscIndex = rangeMappers_.ghostIndex( entry.first );
355 if( rangeAuxiliaryDofs.contains( petscIndex ) )
358 for (
unsigned int rb = 0; rb<rangeLocalBlockSize/bs; ++rb)
360 const size_t row = petscIndex*rangeLocalBlockSize/bs + rb;
363 assert( row <
size_t(localRows/bs) );
364 d_nnz[ row ] = o_nnz[ row ] = 0;
365 for(
const auto local : entry.second )
367 if( !domainAuxiliaryDofs.contains( domainMappers_.ghostIndex( local ) ) )
368 d_nnz[ row ] += domainLocalBlockSize/bs;
370 o_nnz[ row ] += domainLocalBlockSize/bs;
374 ::Dune::Petsc::MatSetUp( petscMatrix_, bs, &d_nnz[0], &o_nnz[0] );
376 sequence_ = domainSpace().sequence();
380 status_ = statAssembled ;
386 ::Dune::Petsc::MatZeroEntries( petscMatrix_ );
390 template <
class Vector>
391 void setUnitRows(
const Vector &rows )
393 std::vector< PetscInt > r( rows.size() );
394 for( std::size_t i =0 ; i< rows.size(); ++i )
396 const PetscInt block = rangeMappers_.parallelIndex( rows[ i ] / rangeLocalBlockSize );
397 r[ i ] = block * rangeLocalBlockSize + (rows[ i ] % rangeLocalBlockSize);
400 if( finalizedAlready() )
402 ::Dune::Petsc::MatZeroRows( petscMatrix_, r.size(), r.data(), 1. );
406 unitRows_.reserve( unitRows_.size() + r.size() );
407 for(
const auto& row : r )
408 unitRows_.push_back( row );
413 ObjectType* newObject()
const
415 return new ObjectType( *
this, domainSpace(), rangeSpace() );
422 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
423 LocalMatrixType localMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
const
425 return LocalMatrixType(localMatrixStack_, domainEntity, rangeEntity);
428 LocalColumnObjectType localColumn(
const DomainEntityType &colEntity )
const
430 return LocalColumnObjectType ( *
this, colEntity );
434 void unitRow(
const PetscInt localRow,
const PetscScalar diag = 1.0 )
436 std::array< PetscInt, domainLocalBlockSize > rows;
437 const PetscInt row = rangeMappers_.parallelIndex( localRow );
438 for(
unsigned int i=0, r = row * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r )
443 if( finalizedAlready() )
446 ::Dune::Petsc::MatZeroRows( petscMatrix_, domainLocalBlockSize, rows.data(), diag );
451 assert(
std::abs( diag - 1. ) < 1e-12 );
452 unitRows_.reserve( unitRows_.size() + domainLocalBlockSize );
453 for(
const auto& r : rows )
455 unitRows_.push_back( r );
460 bool blockedMode()
const {
return blockedMode_; }
463 template<
class PetscOp >
464 void applyToBlock (
const PetscInt localRow,
const PetscInt localCol,
const MatrixBlockType& block, PetscOp op )
467 const PetscInt localCols = domainMappers_.ghostMapper().interiorSize() * domainLocalBlockSize;
468 const PetscInt localRows = rangeMappers_.ghostMapper().interiorSize() * rangeLocalBlockSize;
469 assert( localRow < localRows );
470 assert( localCol < localCols );
476 const PetscInt row = rangeMappers_.parallelIndex( localRow );
477 const PetscInt col = rangeMappers_.parallelIndex( localCol );
478 ::Dune::Petsc::MatSetValuesBlocked( petscMatrix_, 1, &row, 1, &col, block.data(), op );
483 const PetscInt row = rangeMappers_.parallelIndex( localRow );
484 const PetscInt col = rangeMappers_.parallelIndex( localCol );
485 std::array< PetscInt, domainLocalBlockSize > rows;
486 std::array< PetscInt, domainLocalBlockSize > cols;
487 for(
unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
494 ::Dune::Petsc::MatSetValues( petscMatrix_, domainLocalBlockSize, rows.data(), domainLocalBlockSize, cols.data(), block.data(), op );
496 setStatus( statAssembled );
499 template<
class LocalBlock,
class PetscOp >
500 void applyToBlock (
const size_t row,
const size_t col,
const LocalBlock& block, PetscOp op )
502 assert( block.rows() == rangeLocalBlockSize );
503 assert( block.cols() == domainLocalBlockSize );
506 MatrixBlockType matBlock( block );
507 applyToBlock( row, col, matBlock, op );
510 template<
class LocalBlock >
511 void setBlock (
const size_t row,
const size_t col,
const LocalBlock& block )
513#ifndef USE_SMP_PARALLEL
514 assert( status_==statAssembled || status_==statInsert );
516 assert( row < std::numeric_limits< int > :: max() );
517 assert( col < std::numeric_limits< int > :: max() );
519 setStatus( statInsert );
520 applyToBlock(
static_cast< PetscInt
> (row),
static_cast< PetscInt
> (col), block, INSERT_VALUES );
523 template<
class LocalBlock >
524 void addBlock (
const size_t row,
const size_t col,
const LocalBlock& block )
526#ifndef USE_SMP_PARALLEL
527 assert( status_==statAssembled || status_==statInsert );
529 assert( row < std::numeric_limits< int > :: max() );
530 assert( col < std::numeric_limits< int > :: max() );
532 setStatus( statAdd );
533 applyToBlock(
static_cast< PetscInt
> (row),
static_cast< PetscInt
> (col), block, ADD_VALUES );
537 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
540 template <
class Operation>
541 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
542 const TemporaryLocalMatrixType &localMat,
const Operation&,
543 const std::integral_constant< bool, false> nonscaled )
545 return localMat.data();
549 template <
class LM,
class Operation>
550 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
551 const Assembly::Impl::LocalMatrixGetter< LM >& localMat,
const Operation&,
552 const std::integral_constant< bool, false> nonscaled )
554 return localMat.localMatrix().data();
558 template <
class LocalMatrix,
class Operation,
bool T>
559 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
560 const LocalMatrix &localMat,
const Operation& operation,
561 const std::integral_constant< bool, T> scaled )
563 std::vector< PetscScalar >& v = *(v_);
564 v.resize( rSize * cSize );
565 for(
unsigned int i = 0, ic = 0 ; i< rSize; ++i )
567 for(
unsigned int j =0; j< cSize; ++j, ++ic )
569 v[ ic ] = operation( localMat.get( i, j ) );
575 template<
class LocalMatrix,
class Operation,
class PetscOp,
bool T >
576 void applyLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
577 const LocalMatrix &localMat,
const Operation& operation,
579 const std::integral_constant<bool, T> scaled )
581 auto& rcTemp = *(rcTemp_);
582 std::vector< PetscInt >& r = rcTemp.first;
583 std::vector< PetscInt >& c = rcTemp.second;
587 setupIndicesBlocked( rangeMappers_, rangeEntity, r );
588 setupIndicesBlocked( domainMappers_, domainEntity, c );
591 const unsigned int rSize = r.size() * domainLocalBlockSize ;
592 const unsigned int cSize = c.size() * domainLocalBlockSize ;
594 const PetscScalar* values = getValues( rSize, cSize, localMat, operation, scaled );
595 ::Dune::Petsc::MatSetValuesBlocked( petscMatrix_, r.size(), r.data(), c.size(), c.data(), values, petscOp );
599 setupIndices( rangeMappers_, rangeEntity, r );
600 setupIndices( domainMappers_, domainEntity, c );
602 const unsigned int rSize = r.size();
603 const unsigned int cSize = c.size();
605 const PetscScalar* values = getValues( rSize, cSize, localMat, operation, scaled );
606 ::Dune::Petsc::MatSetValues( petscMatrix_, rSize, r.data(), cSize, c.data(), values, petscOp );
608 setStatus( statAssembled );
612 template<
class LocalMatrix >
613 void addLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
615#ifndef USE_SMP_PARALLEL
616 assert( status_==statAssembled || status_==statAdd );
618 setStatus( statAdd );
620 auto operation = [] (
const PetscScalar& value ) -> PetscScalar {
return value; };
622 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, ADD_VALUES, std::integral_constant< bool, false>() );
625 template<
class LocalMatrix,
class Scalar >
626 void addScaledLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat,
const Scalar &s )
628#ifndef USE_SMP_PARALLEL
629 assert( status_==statAssembled || status_==statAdd );
631 setStatus( statAdd );
633 auto operation = [ &s ] (
const PetscScalar& value ) -> PetscScalar {
return s * value; };
635 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, ADD_VALUES, std::integral_constant< bool, true>() );
638 template<
class LocalMatrix >
639 void setLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
641#ifndef USE_SMP_PARALLEL
642 assert( status_==statAssembled || status_==statInsert );
644 setStatus( statInsert );
646 auto operation = [] (
const PetscScalar& value ) -> PetscScalar {
return value; };
648 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, INSERT_VALUES, std::integral_constant< bool, false>() );
651 template<
class LocalMatrix >
652 void getLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity, LocalMatrix &localMat )
const
657#ifndef USE_SMP_PARALLEL
658 assert( status_==statAssembled || status_==statGet );
660 setStatus( statGet );
662 auto& rcTemp = *(rcTemp_);
663 std::vector< PetscInt >& r = rcTemp.first;
664 std::vector< PetscInt >& c = rcTemp.second;
665 std::vector< PetscScalar >& v = *(v_);
667 setupIndices( rangeMappers_, rangeEntity, r );
668 setupIndices( domainMappers_, domainEntity, c );
670 v.resize( r.size() * c.size() );
671 ::Dune::Petsc::MatGetValues( petscMatrix_, r.size(), r.data(), c.size(), c.data(), v.data() );
673 for( std::size_t i =0 ; i< r.size(); ++i )
674 for( std::size_t j =0; j< c.size(); ++j )
675 localMat.set( i, j, v[ i * c.size() +j ] );
677 setStatus( statAssembled );
680 void scaleLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const RangeFieldType &s )
685#ifndef USE_SMP_PARALLEL
686 assert( status_==statAssembled || status_==statGet );
688 setStatus( statGet );
690 auto& rcTemp = *(rcTemp_);
691 std::vector< PetscInt >& r = rcTemp.first;
692 std::vector< PetscInt >& c = rcTemp.second;
693 std::vector< PetscScalar >& v = *(v_);
695 setupIndices( rangeMappers_, rangeEntity, r );
696 setupIndices( domainMappers_, domainEntity, c );
698 const unsigned int rSize = r.size();
699 const unsigned int cSize = c.size();
700 const std::size_t size = rSize * cSize;
704 ::Dune::Petsc::MatGetValues( petscMatrix_, r.size(), r.data(), c.size(), c.data(), v.data() );
707 for( std::size_t i=0; i<size; ++i )
711 ::Dune::Petsc::MatSetValues( petscMatrix_, rSize, r.data(), cSize, c.data(), v.data(), INSERT_VALUES );
712 setStatus( statAssembled );
718 ::Dune::Petsc::MatView( petscMatrix_, PETSC_VIEWER_STDOUT_WORLD );
722 void print( std::ostream& s )
const
724 if( &s == &std::cout || &s == &std::cerr )
731 Mat& exportMatrix ()
const
739 PetscLinearOperator ();
744 if( status_ != statNothing )
746 ::Dune::Petsc::MatDestroy( &petscMatrix_ );
747 status_ = statNothing ;
752 void setStatus (
const Status &newstatus)
const
755#ifndef USE_SMP_PARALLEL
760 template<
class DFS,
class Entity >
761 static void setupIndices (
const PetscMappers< DFS > &mappers,
const Entity &entity, std::vector< PetscInt > &indices )
763 NonBlockMapper< const typename PetscMappers< DFS >::ParallelMapperType, DFS::localBlockSize > nonBlockMapper( mappers.parallelMapper() );
764 nonBlockMapper.map( entity, indices );
767 template<
class DFS,
class Entity >
768 static void setupIndicesBlocked (
const PetscMappers< DFS > &mappers,
const Entity &entity, std::vector< PetscInt > &indices )
770 mappers.parallelMapper().map( entity, indices );
776 DomainMappersType domainMappers_;
777 RangeMappersType rangeMappers_;
780 mutable Mat petscMatrix_;
782 mutable LocalMatrixStackType localMatrixStack_;
783 mutable Status status_;
785 const bool viennaCL_;
786 const bool blockedMode_;
788 mutable std::unique_ptr< PetscDomainFunctionType > petscArg_;
789 mutable std::unique_ptr< PetscRangeFunctionType > petscDest_;
791 mutable ThreadSafeValue< std::vector< PetscScalar > > v_;
792 mutable ThreadSafeValue< std::pair< std::vector< PetscInt >, std::vector< PetscInt > > > rcTemp_;
794 mutable std::vector< PetscInt > unitRows_;
802 template<
typename DomainFunction,
typename RangeFunction >
803 class PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrix
804 :
public LocalMatrixDefault< typename PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrixTraits >
806 typedef LocalMatrix ThisType;
807 typedef LocalMatrixDefault< typename PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrixTraits > BaseType;
809 typedef PetscLinearOperator< DomainFunction, RangeFunction > PetscLinearOperatorType;
813 typedef PetscInt DofIndexType;
814 typedef std::vector< DofIndexType > IndexVectorType;
815 typedef typename DomainFunction::DiscreteFunctionSpaceType DomainSpaceType;
816 typedef typename RangeFunction::DiscreteFunctionSpaceType RangeSpaceType;
817 typedef typename DomainSpaceType::BasisFunctionSetType DomainBasisFunctionSetType;
818 typedef typename RangeSpaceType::BasisFunctionSetType RangeBasisFunctionSetType;
820 enum { rangeBlockSize = RangeSpaceType::localBlockSize };
821 enum { domainBlockSize = DomainSpaceType ::localBlockSize };
826 template<
typename PetscMapping >
827 struct PetscAssignFunctor
829 explicit PetscAssignFunctor (
const PetscMapping &petscMapping, IndexVectorType &indices )
830 : petscMapping_( petscMapping ),
834 template<
typename T >
835 void operator() (
const std::size_t localIndex, T globalIndex ) { indices_[ localIndex ] = petscMapping_.globalMapping( globalIndex ); }
838 const PetscMapping &petscMapping_;
839 IndexVectorType &indices_;
843 [[deprecated(
"Use TemporaryLocal Matrix and {add,set,get}LocalMatrix")]]
844 LocalMatrix (
const PetscLinearOperatorType &petscLinOp,
845 const DomainSpaceType &domainSpace,
846 const RangeSpaceType &rangeSpace )
847 : BaseType( domainSpace, rangeSpace ),
848 petscLinearOperator_( petscLinOp )
855 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
858 BaseType :: init( domainEntity, rangeEntity );
868 setupIndices( petscLinearOperator_.rangeMappers_, rangeEntity, rowIndices_ );
869 setupIndices( petscLinearOperator_.domainMappers_, domainEntity, colIndices_ );
871 status_ = statAssembled;
872 petscLinearOperator_.setStatus(status_);
875 inline void add (
const int localRow,
const int localCol,
const RangeFieldType &value )
877 assert( status_==statAssembled || status_==statAdd );
879 petscLinearOperator_.setStatus(status_);
880 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIndex( localRow ), globalColIndex( localCol ) , value, ADD_VALUES );
883 inline void set(
const int localRow,
const int localCol,
const RangeFieldType &value )
885 assert( status_==statAssembled || status_==statInsert );
886 status_ = statInsert;
887 petscLinearOperator_.setStatus(status_);
888 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIndex( localRow ), globalColIndex( localCol ) , value, INSERT_VALUES );
892 void clearRow (
const int localRow )
894 assert( status_==statAssembled || status_==statInsert );
895 status_ = statInsert;
896 petscLinearOperator_.setStatus(status_);
897 const int col = this->columns();
898 const int globalRowIdx = globalRowIndex( localRow );
899 for(
int localCol=0; localCol<col; ++localCol)
901 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIdx, globalColIndex( localCol ), 0.0, INSERT_VALUES );
905 inline const RangeFieldType
get (
const int localRow,
const int localCol )
const
907 assert( status_==statAssembled || status_==statGet );
909 petscLinearOperator_.setStatus(status_);
911 const int r[] = {globalRowIndex( localRow )};
912 const int c[] = {globalColIndex( localCol )};
913 ::Dune::Petsc::MatGetValues( petscMatrix(), 1, r, 1, c, v );
917 inline void scale (
const RangeFieldType &factor )
919 const_cast< PetscLinearOperatorType&
> (petscLinearOperator_).scaleLocalMatrix( this->domainEntity(), this->rangeEntity(), factor );
925 Mat& petscMatrix () {
return petscLinearOperator_.petscMatrix_; }
926 const Mat& petscMatrix ()
const {
return petscLinearOperator_.petscMatrix_; }
929 int rows()
const {
return rowIndices_.size(); }
930 int columns()
const {
return colIndices_.size(); }
933 DofIndexType globalRowIndex(
const int localRow )
const
935 assert( localRow <
static_cast< int >( rowIndices_.size() ) );
936 return rowIndices_[ localRow ];
939 DofIndexType globalColIndex(
const int localCol )
const
941 assert( localCol <
static_cast< int >( colIndices_.size() ) );
942 return colIndices_[ localCol ];
948 const PetscLinearOperatorType &petscLinearOperator_;
949 IndexVectorType rowIndices_;
950 IndexVectorType colIndices_;
951 mutable Status status_;
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122