1#ifndef DUNE_FEM_PETSCLINEAROPERATOR_HH
2#define DUNE_FEM_PETSCLINEAROPERATOR_HH
15#include <dune/fem/function/petscdiscretefunction/petscdiscretefunction.hh>
16#include <dune/fem/misc/functor.hh>
17#include <dune/fem/misc/fmatrixconverter.hh>
18#include <dune/fem/misc/petsc/petsccommon.hh>
19#include <dune/fem/misc/threads/threadsafevalue.hh>
20#include <dune/fem/operator/common/localcontribution.hh>
21#include <dune/fem/operator/common/localmatrix.hh>
22#include <dune/fem/operator/common/localmatrixwrapper.hh>
23#include <dune/fem/operator/common/temporarylocalmatrix.hh>
24#include <dune/fem/operator/common/operator.hh>
25#include <dune/fem/operator/common/stencil.hh>
26#include <dune/fem/operator/matrix/columnobject.hh>
27#include <dune/fem/operator/matrix/functor.hh>
28#include <dune/fem/space/mapper/nonblockmapper.hh>
29#include <dune/fem/space/mapper/petsc.hh>
30#include <dune/fem/storage/objectstack.hh>
32#include <dune/fem/solver/parameter.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 setUnitRowsImpl( unitRows_, auxRows_ );
221 bool finalizedAlready()
const
223 PetscBool assembled = PETSC_FALSE ;
224 ::Dune::Petsc::MatAssembled( petscMatrix_, &assembled );
225 return assembled == PETSC_TRUE;
228 void finalizeAssembly ()
const
230 const_cast< ThisType&
> (*this).finalize();
234 const DomainSpaceType &domainSpace ()
const {
return domainMappers_.space(); }
235 const RangeSpaceType &rangeSpace ()
const {
return rangeMappers_.space(); }
240 template <
class DF,
class RF>
241 void apply (
const DF &arg, RF &dest )
const
244 petscArg_.reset(
new PetscDomainFunctionType(
"PetscOp-arg", domainSpace() ) );
246 petscDest_.reset(
new PetscRangeFunctionType(
"PetscOp-arg", rangeSpace() ) );
248 petscArg_->assign( arg );
249 apply( *petscArg_, *petscDest_ );
250 dest.assign( *petscDest_ );
254 void apply (
const PetscDomainFunctionType &arg, PetscRangeFunctionType &dest )
const
258 ::Dune::Petsc::MatMult( petscMatrix_, *arg.petscVec() , *dest.petscVec() );
261 void operator() (
const DomainFunctionType &arg, RangeFunctionType &dest )
const
268 reserve( SimpleStencil<DomainSpaceType,RangeSpaceType>(0),
true );
272 void reserve (
const std::vector< Set >& sparsityPattern )
274 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ),
true );
278 template <
class StencilType>
279 void reserve (
const StencilType &stencil,
const bool isSimpleStencil =
false )
281 domainMappers_.update();
282 rangeMappers_.update();
284 if(sequence_ != domainSpace().sequence())
294 ::Dune::Petsc::MatCreate( domainSpace().gridPart().comm(), &petscMatrix_ );
297 const PetscInt localCols = domainMappers_.ghostMapper().interiorSize() * domainLocalBlockSize;
298 const PetscInt localRows = rangeMappers_.ghostMapper().interiorSize() * rangeLocalBlockSize;
300 const PetscInt globalCols = domainMappers_.parallelMapper().size() * domainLocalBlockSize;
301 const PetscInt globalRows = rangeMappers_.parallelMapper().size() * rangeLocalBlockSize;
303 ::Dune::Petsc::MatSetSizes( petscMatrix_, localRows, localCols, globalRows, globalCols );
308 ::Dune::Petsc::MatSetType( petscMatrix_, MATAIJVIENNACL );
310 else if( blockedMode_ )
312 bs = domainLocalBlockSize ;
313 ::Dune::Petsc::MatSetType( petscMatrix_, MATBAIJ );
315 ::Dune::Petsc::MatSetBlockSize( petscMatrix_, bs );
319 ::Dune::Petsc::MatSetType( petscMatrix_, MATAIJ );
337 if( isSimpleStencil || std::is_same< StencilType,SimpleStencil<DomainSpaceType,RangeSpaceType> >::value )
339 ::Dune::Petsc::MatSetUp( petscMatrix_, bs, domainLocalBlockSize * stencil.maxNonZerosEstimate() );
343 DomainAuxiliaryDofsType domainAuxiliaryDofs( domainMappers_.ghostMapper() );
344 RangeAuxiliaryDofsType rangeAuxiliaryDofs( rangeMappers_.ghostMapper() );
346 std::vector< PetscInt > d_nnz( localRows / bs, 0 );
347 std::vector< PetscInt > o_nnz( localRows / bs, 0 );
348 for(
const auto& entry : stencil.globalStencil() )
350 const int petscIndex = rangeMappers_.ghostIndex( entry.first );
351 if( rangeAuxiliaryDofs.contains( petscIndex ) )
354 for (
unsigned int rb = 0; rb<rangeLocalBlockSize/bs; ++rb)
356 const size_t row = petscIndex*rangeLocalBlockSize/bs + rb;
359 assert( row <
size_t(localRows/bs) );
360 d_nnz[ row ] = o_nnz[ row ] = 0;
361 for(
const auto local : entry.second )
363 if( !domainAuxiliaryDofs.contains( domainMappers_.ghostIndex( local ) ) )
364 d_nnz[ row ] += domainLocalBlockSize/bs;
366 o_nnz[ row ] += domainLocalBlockSize/bs;
370 ::Dune::Petsc::MatSetUp( petscMatrix_, bs, &d_nnz[0], &o_nnz[0] );
372 sequence_ = domainSpace().sequence();
376 status_ = statAssembled ;
379 virtual void clear ()
382 ::Dune::Petsc::MatZeroEntries( petscMatrix_ );
389 void setUnitRowsImpl(
const std::vector< PetscInt >& r,
390 const std::vector< PetscInt >& a )
392 ::Dune::Petsc::MatZeroRows( petscMatrix_, r.size(), r.data(), 1.0 );
393 ::Dune::Petsc::MatZeroRows( petscMatrix_, a.size(), a.data(), 0.0 );
396 template <
class Container>
397 void setUnitRows(
const Container& unitRows,
const Container& auxRows )
399 std::vector< PetscInt > r, a;
400 r.reserve( unitRows.size() );
401 a.reserve( auxRows.size() );
403 auto setupRows = [
this] (
const Container& rows, std::vector< PetscInt >& r )
405 for(
const auto& row : rows )
407 const PetscInt block = this->rangeMappers_.parallelIndex( row / rangeLocalBlockSize );
408 r.push_back( block * rangeLocalBlockSize + (row % rangeLocalBlockSize) );
412 setupRows( unitRows, r );
413 setupRows( auxRows, a );
415 if( finalizedAlready() )
417 setUnitRowsImpl( r, a );
421 unitRows_.reserve( unitRows_.size() + r.size() );
422 for(
const auto& row : r )
423 unitRows_.push_back( row );
425 auxRows_.reserve( auxRows_.size() + a.size() );
426 for(
const auto& row : a )
427 auxRows_.push_back( row );
432 ObjectType* newObject()
const
434 return new ObjectType( *
this, domainSpace(), rangeSpace() );
441 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
442 LocalMatrixType localMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
const
444 return LocalMatrixType(localMatrixStack_, domainEntity, rangeEntity);
447 LocalColumnObjectType localColumn(
const DomainEntityType &colEntity )
const
449 return LocalColumnObjectType ( *
this, colEntity );
453 void unitRow(
const PetscInt localRow,
const PetscScalar diag = 1.0 )
455 std::array< PetscInt, domainLocalBlockSize > rows;
456 const PetscInt row = rangeMappers_.parallelIndex( localRow );
457 for(
unsigned int i=0, r = row * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r )
462 if( finalizedAlready() )
465 ::Dune::Petsc::MatZeroRows( petscMatrix_, domainLocalBlockSize, rows.data(), diag );
470 assert( std::abs( diag - 1. ) < 1e-12 );
471 unitRows_.reserve( unitRows_.size() + domainLocalBlockSize );
472 for(
const auto& r : rows )
474 unitRows_.push_back( r );
479 bool blockedMode()
const {
return blockedMode_; }
482 template<
class PetscOp >
483 void applyToBlock (
const PetscInt localRow,
const PetscInt localCol,
const MatrixBlockType& block, PetscOp op )
486 const PetscInt localCols = domainMappers_.ghostMapper().interiorSize() * domainLocalBlockSize;
487 const PetscInt localRows = rangeMappers_.ghostMapper().interiorSize() * rangeLocalBlockSize;
488 assert( localRow < localRows );
489 assert( localCol < localCols );
495 const PetscInt row = rangeMappers_.parallelIndex( localRow );
496 const PetscInt col = rangeMappers_.parallelIndex( localCol );
497 ::Dune::Petsc::MatSetValuesBlocked( petscMatrix_, 1, &row, 1, &col, block.data(), op );
502 const PetscInt row = rangeMappers_.parallelIndex( localRow );
503 const PetscInt col = rangeMappers_.parallelIndex( localCol );
504 std::array< PetscInt, domainLocalBlockSize > rows;
505 std::array< PetscInt, domainLocalBlockSize > cols;
506 for(
unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
513 ::Dune::Petsc::MatSetValues( petscMatrix_, domainLocalBlockSize, rows.data(), domainLocalBlockSize, cols.data(), block.data(), op );
515 setStatus( statAssembled );
518 template<
class LocalBlock,
class PetscOp >
519 void applyToBlock (
const size_t row,
const size_t col,
const LocalBlock& block, PetscOp op )
521 assert( block.rows() == rangeLocalBlockSize );
522 assert( block.cols() == domainLocalBlockSize );
525 MatrixBlockType matBlock( block );
526 applyToBlock( row, col, matBlock, op );
529 template<
class LocalBlock >
530 void setBlock (
const size_t row,
const size_t col,
const LocalBlock& block )
532#ifndef USE_SMP_PARALLEL
533 assert( status_==statAssembled || status_==statInsert );
538 setStatus( statInsert );
539 applyToBlock(
static_cast< PetscInt
> (row),
static_cast< PetscInt
> (col), block, INSERT_VALUES );
542 template<
class LocalBlock >
543 void addBlock (
const size_t row,
const size_t col,
const LocalBlock& block )
545#ifndef USE_SMP_PARALLEL
546 assert( status_==statAssembled || status_==statInsert );
551 setStatus( statAdd );
552 applyToBlock(
static_cast< PetscInt
> (row),
static_cast< PetscInt
> (col), block, ADD_VALUES );
556 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
559 template <
class Operation>
560 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
561 const TemporaryLocalMatrixType &localMat,
const Operation&,
562 const std::integral_constant< bool, false> nonscaled )
564 return localMat.data();
568 template <
class LM,
class Operation>
569 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
570 const Assembly::Impl::LocalMatrixGetter< LM >& localMat,
const Operation&,
571 const std::integral_constant< bool, false> nonscaled )
573 return localMat.localMatrix().data();
577 template <
class LocalMatrix,
class Operation,
bool T>
578 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
579 const LocalMatrix &localMat,
const Operation& operation,
580 const std::integral_constant< bool, T> scaled )
582 std::vector< PetscScalar >& v = *(v_);
583 v.resize( rSize * cSize );
584 for(
unsigned int i = 0, ic = 0 ; i< rSize; ++i )
586 for(
unsigned int j =0; j< cSize; ++j, ++ic )
588 v[ ic ] = operation( localMat.get( i, j ) );
594 template<
class LocalMatrix,
class Operation,
class PetscOp,
bool T >
595 void applyLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
596 const LocalMatrix &localMat,
const Operation& operation,
598 const std::integral_constant<bool, T> scaled )
600 auto& rcTemp = *(rcTemp_);
601 std::vector< PetscInt >& r = rcTemp.first;
602 std::vector< PetscInt >& c = rcTemp.second;
606 setupIndicesBlocked( rangeMappers_, rangeEntity, r );
607 setupIndicesBlocked( domainMappers_, domainEntity, c );
610 const unsigned int rSize = r.size() * domainLocalBlockSize ;
611 const unsigned int cSize = c.size() * domainLocalBlockSize ;
613 const PetscScalar* values = getValues( rSize, cSize, localMat, operation, scaled );
614 ::Dune::Petsc::MatSetValuesBlocked( petscMatrix_, r.size(), r.data(), c.size(), c.data(), values, petscOp );
618 setupIndices( rangeMappers_, rangeEntity, r );
619 setupIndices( domainMappers_, domainEntity, c );
621 const unsigned int rSize = r.size();
622 const unsigned int cSize = c.size();
624 const PetscScalar* values = getValues( rSize, cSize, localMat, operation, scaled );
625 ::Dune::Petsc::MatSetValues( petscMatrix_, rSize, r.data(), cSize, c.data(), values, petscOp );
627 setStatus( statAssembled );
631 template<
class LocalMatrix >
632 void addLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
634#ifndef USE_SMP_PARALLEL
635 assert( status_==statAssembled || status_==statAdd );
637 setStatus( statAdd );
639 auto operation = [] (
const PetscScalar& value ) -> PetscScalar {
return value; };
641 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, ADD_VALUES, std::integral_constant< bool, false>() );
644 template<
class LocalMatrix,
class Scalar >
645 void addScaledLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat,
const Scalar &s )
647#ifndef USE_SMP_PARALLEL
648 assert( status_==statAssembled || status_==statAdd );
650 setStatus( statAdd );
652 auto operation = [ &s ] (
const PetscScalar& value ) -> PetscScalar {
return s * value; };
654 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, ADD_VALUES, std::integral_constant< bool, true>() );
657 template<
class LocalMatrix >
658 void setLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
660#ifndef USE_SMP_PARALLEL
661 assert( status_==statAssembled || status_==statInsert );
663 setStatus( statInsert );
665 auto operation = [] (
const PetscScalar& value ) -> PetscScalar {
return value; };
667 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, INSERT_VALUES, std::integral_constant< bool, false>() );
670 template<
class LocalMatrix >
671 void getLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity, LocalMatrix &localMat )
const
676#ifndef USE_SMP_PARALLEL
677 assert( status_==statAssembled || status_==statGet );
679 setStatus( statGet );
681 auto& rcTemp = *(rcTemp_);
682 std::vector< PetscInt >& r = rcTemp.first;
683 std::vector< PetscInt >& c = rcTemp.second;
684 std::vector< PetscScalar >& v = *(v_);
686 setupIndices( rangeMappers_, rangeEntity, r );
687 setupIndices( domainMappers_, domainEntity, c );
689 v.resize( r.size() * c.size() );
690 ::Dune::Petsc::MatGetValues( petscMatrix_, r.size(), r.data(), c.size(), c.data(), v.data() );
692 for( std::size_t i =0 ; i< r.size(); ++i )
693 for( std::size_t j =0; j< c.size(); ++j )
694 localMat.set( i, j, v[ i * c.size() +j ] );
696 setStatus( statAssembled );
699 void scaleLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const RangeFieldType &s )
704#ifndef USE_SMP_PARALLEL
705 assert( status_==statAssembled || status_==statGet );
707 setStatus( statGet );
709 auto& rcTemp = *(rcTemp_);
710 std::vector< PetscInt >& r = rcTemp.first;
711 std::vector< PetscInt >& c = rcTemp.second;
712 std::vector< PetscScalar >& v = *(v_);
714 setupIndices( rangeMappers_, rangeEntity, r );
715 setupIndices( domainMappers_, domainEntity, c );
717 const unsigned int rSize = r.size();
718 const unsigned int cSize = c.size();
719 const std::size_t
size = rSize * cSize;
723 ::Dune::Petsc::MatGetValues( petscMatrix_, r.size(), r.data(), c.size(), c.data(), v.data() );
726 for( std::size_t i=0; i<
size; ++i )
730 ::Dune::Petsc::MatSetValues( petscMatrix_, rSize, r.data(), cSize, c.data(), v.data(), INSERT_VALUES );
731 setStatus( statAssembled );
737 ::Dune::Petsc::MatView( petscMatrix_, PETSC_VIEWER_STDOUT_WORLD );
741 void print( std::ostream& s )
const
743 if( &s == &std::cout || &s == &std::cerr )
750 Mat& exportMatrix ()
const
758 PetscLinearOperator ();
763 if( status_ != statNothing )
765 ::Dune::Petsc::MatDestroy( &petscMatrix_ );
766 status_ = statNothing ;
773 void setStatus (
const Status &newstatus)
const
776#ifndef USE_SMP_PARALLEL
781 template<
class DFS,
class Entity >
782 static void setupIndices (
const PetscMappers< DFS > &mappers,
const Entity &entity, std::vector< PetscInt > &indices )
784 NonBlockMapper< const typename PetscMappers< DFS >::ParallelMapperType, DFS::localBlockSize > nonBlockMapper( mappers.parallelMapper() );
785 nonBlockMapper.map( entity, indices );
788 template<
class DFS,
class Entity >
789 static void setupIndicesBlocked (
const PetscMappers< DFS > &mappers,
const Entity &entity, std::vector< PetscInt > &indices )
791 mappers.parallelMapper().map( entity, indices );
797 DomainMappersType domainMappers_;
798 RangeMappersType rangeMappers_;
801 mutable Mat petscMatrix_;
803 mutable LocalMatrixStackType localMatrixStack_;
804 mutable Status status_;
806 const bool viennaCL_;
807 const bool blockedMode_;
809 mutable std::unique_ptr< PetscDomainFunctionType > petscArg_;
810 mutable std::unique_ptr< PetscRangeFunctionType > petscDest_;
812 mutable ThreadSafeValue< std::vector< PetscScalar > > v_;
813 mutable ThreadSafeValue< std::pair< std::vector< PetscInt >, std::vector< PetscInt > > > rcTemp_;
815 mutable std::vector< PetscInt > unitRows_;
816 mutable std::vector< PetscInt > auxRows_;
824 template<
typename DomainFunction,
typename RangeFunction >
825 class PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrix
826 :
public LocalMatrixDefault< typename PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrixTraits >
828 typedef LocalMatrix ThisType;
829 typedef LocalMatrixDefault< typename PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrixTraits > BaseType;
831 typedef PetscLinearOperator< DomainFunction, RangeFunction > PetscLinearOperatorType;
835 typedef PetscInt DofIndexType;
836 typedef std::vector< DofIndexType > IndexVectorType;
837 typedef typename DomainFunction::DiscreteFunctionSpaceType DomainSpaceType;
838 typedef typename RangeFunction::DiscreteFunctionSpaceType RangeSpaceType;
839 typedef typename DomainSpaceType::BasisFunctionSetType DomainBasisFunctionSetType;
840 typedef typename RangeSpaceType::BasisFunctionSetType RangeBasisFunctionSetType;
842 enum { rangeBlockSize = RangeSpaceType::localBlockSize };
843 enum { domainBlockSize = DomainSpaceType ::localBlockSize };
848 template<
typename PetscMapping >
849 struct PetscAssignFunctor
851 explicit PetscAssignFunctor (
const PetscMapping &petscMapping, IndexVectorType &indices )
852 : petscMapping_( petscMapping ),
856 template<
typename T >
857 void operator() (
const std::size_t localIndex, T globalIndex ) { indices_[ localIndex ] = petscMapping_.globalMapping( globalIndex ); }
860 const PetscMapping &petscMapping_;
861 IndexVectorType &indices_;
865 [[deprecated(
"Use TemporaryLocal Matrix and {add,set,get}LocalMatrix")]]
866 LocalMatrix (
const PetscLinearOperatorType &petscLinOp,
867 const DomainSpaceType &domainSpace,
868 const RangeSpaceType &rangeSpace )
869 : BaseType( domainSpace, rangeSpace ),
870 petscLinearOperator_( petscLinOp )
877 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
880 BaseType :: init( domainEntity, rangeEntity );
890 setupIndices( petscLinearOperator_.rangeMappers_, rangeEntity, rowIndices_ );
891 setupIndices( petscLinearOperator_.domainMappers_, domainEntity, colIndices_ );
893 status_ = statAssembled;
894 petscLinearOperator_.setStatus(status_);
897 inline void add (
const int localRow,
const int localCol,
const RangeFieldType &value )
899 assert( status_==statAssembled || status_==statAdd );
901 petscLinearOperator_.setStatus(status_);
902 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIndex( localRow ), globalColIndex( localCol ) , value, ADD_VALUES );
905 inline void set(
const int localRow,
const int localCol,
const RangeFieldType &value )
907 assert( status_==statAssembled || status_==statInsert );
908 status_ = statInsert;
909 petscLinearOperator_.setStatus(status_);
910 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIndex( localRow ), globalColIndex( localCol ) , value, INSERT_VALUES );
914 void clearRow (
const int localRow )
916 assert( status_==statAssembled || status_==statInsert );
917 status_ = statInsert;
918 petscLinearOperator_.setStatus(status_);
919 const int col = this->columns();
920 const int globalRowIdx = globalRowIndex( localRow );
921 for(
int localCol=0; localCol<col; ++localCol)
923 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIdx, globalColIndex( localCol ), 0.0, INSERT_VALUES );
927 inline const RangeFieldType
get (
const int localRow,
const int localCol )
const
929 assert( status_==statAssembled || status_==statGet );
931 petscLinearOperator_.setStatus(status_);
933 const int r[] = {globalRowIndex( localRow )};
934 const int c[] = {globalColIndex( localCol )};
935 ::Dune::Petsc::MatGetValues( petscMatrix(), 1, r, 1, c, v );
939 inline void scale (
const RangeFieldType &factor )
941 const_cast< PetscLinearOperatorType&
> (petscLinearOperator_).scaleLocalMatrix( this->domainEntity(), this->rangeEntity(), factor );
947 Mat& petscMatrix () {
return petscLinearOperator_.petscMatrix_; }
948 const Mat& petscMatrix ()
const {
return petscLinearOperator_.petscMatrix_; }
951 int rows()
const {
return rowIndices_.size(); }
952 int columns()
const {
return colIndices_.size(); }
955 DofIndexType globalRowIndex(
const int localRow )
const
957 assert( localRow <
static_cast< int >( rowIndices_.size() ) );
958 return rowIndices_[ localRow ];
961 DofIndexType globalColIndex(
const int localCol )
const
963 assert( localCol <
static_cast< int >( colIndices_.size() ) );
964 return colIndices_[ localCol ];
970 const PetscLinearOperatorType &petscLinearOperator_;
971 IndexVectorType rowIndices_;
972 IndexVectorType colIndices_;
973 mutable Status status_;
This file implements a dense matrix with dynamic numbers of rows and columns.
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22