MultiFieldMatrix Class Reference

#include <MultiFieldMatrix.h>

Public Types
typedef MultiField::ArrayIndex	Index
typedef pair< Index, Index >	EntryIndex
typedef pair< const StorageSite *, const StorageSite * >	SSPair
typedef map< EntryIndex, shared_ptr< Matrix > >	MatrixMap
typedef map< Index, int >	MatrixSizeMap
typedef map< SSPair, shared_ptr< Array< int > > >	MatrixMappersMap
typedef map< Index, shared_ptr < StorageSite > >	StorageSiteMap
typedef map< const StorageSite *, shared_ptr< StorageSite > >	GhostStorageSiteMap
typedef map< Index, shared_ptr < CRConnectivity > >	CoarseToFineMappingMap
typedef map< EntryIndex, shared_ptr< CRConnectivity > >	CoarseConnectivitiesMap

Public Member Functions
	MultiFieldMatrix ()
virtual	~MultiFieldMatrix ()
	DEFINE_TYPENAME ("MultiFieldMatrix")
void	initAssembly ()
virtual void	multiply (IContainer &yB, const IContainer &xB) const
virtual void	multiplyAndAdd (IContainer &yB, const IContainer &xB) const
virtual void	forwardGS (IContainer &xB, const IContainer &bB, IContainer &temp) const
virtual void	reverseGS (IContainer &xB, const IContainer &bB, IContainer &temp) const
virtual void	Jacobi (IContainer &xB, const IContainer &bB, IContainer &tempB) const
virtual void	iluSolve (IContainer &xB, const IContainer &bB, IContainer &tempB) const
virtual void	spikeSolve (IContainer &xB, const IContainer &bB, IContainer &tempB, const SpikeStorage &spike_storage) const
virtual void	solveBoundary (IContainer &xB, const IContainer &bB, IContainer &temp) const
virtual void	computeResidual (const IContainer &xB, const IContainer &bB, IContainer &rB) const
bool	hasMatrix (const Index &rowIndex, const Index &colIndex) const
Matrix &	getMatrix (const Index &rowIndex, const Index &colIndex)
const Matrix &	getMatrix (const Index &rowIndex, const Index &colIndex) const
void	removeMatrix (const Index &rowIndex, const Index &colIndex)
void	createCoarsening (MultiField &coarseIndex, const int groupSize, const double weightRatioThreshold)
void	syncGhostCoarsening (MultiField &coarseIndexField)
void	createCoarseToFineMapping (const MultiField &coarseIndexField)
void	createCoarseConnectivity (MultiField &coarseIndex)
void	createCoarseMatrices (MultiField &coarseIndex)
void	transpose ()
void	injectResidual (const MultiField &coarseIndex, const MultiField &fineResidualField, MultiField &coarseBField)
void	correctSolution (const MultiField &coarseIndex, MultiField &fineSolutionField, MFRPtr scaleField, const MultiField &coarseSolutionField)
int	getSize () const
int	getLocalSize () const
void	addMatrix (const Index &rowI, const Index &colI, shared_ptr< Matrix > m)
MFRPtr	quadProduct (const MultiField &x) const

Private Attributes
MatrixMap	_matrices
MatrixSizeMap	_coarseSizes
MatrixSizeMap	_coarseGhostSizes
MatrixMappersMap	_coarseScatterMaps
MatrixMappersMap	_coarseGatherMaps
StorageSiteMap	_coarseSites
CoarseToFineMappingMap	_coarseToFineMappings
CoarseConnectivitiesMap	_coarseConnectivities
MatrixMap	_coarseMatrices

Friends
class	LinearSystem
class	LinearSystemMerger

Detailed Description

Definition at line 49 of file MultiFieldMatrix.h.

Member Typedef Documentation

typedef map<EntryIndex,shared_ptr<CRConnectivity> > MultiFieldMatrix::CoarseConnectivitiesMap

Definition at line 70 of file MultiFieldMatrix.h.

typedef map<Index,shared_ptr<CRConnectivity> > MultiFieldMatrix::CoarseToFineMappingMap

Definition at line 69 of file MultiFieldMatrix.h.

typedef pair<Index,Index> MultiFieldMatrix::EntryIndex

Definition at line 58 of file MultiFieldMatrix.h.

typedef map<const StorageSite*,shared_ptr<StorageSite> > MultiFieldMatrix::GhostStorageSiteMap

Definition at line 67 of file MultiFieldMatrix.h.

typedef MultiField::ArrayIndex MultiFieldMatrix::Index

Definition at line 56 of file MultiFieldMatrix.h.

typedef map<EntryIndex,shared_ptr<Matrix> > MultiFieldMatrix::MatrixMap

Definition at line 61 of file MultiFieldMatrix.h.

typedef map<SSPair,shared_ptr<Array<int> > > MultiFieldMatrix::MatrixMappersMap

Definition at line 64 of file MultiFieldMatrix.h.

typedef map<Index,int> MultiFieldMatrix::MatrixSizeMap

Definition at line 62 of file MultiFieldMatrix.h.

typedef pair<const StorageSite*, const StorageSite*> MultiFieldMatrix::SSPair

Definition at line 59 of file MultiFieldMatrix.h.

typedef map<Index,shared_ptr<StorageSite> > MultiFieldMatrix::StorageSiteMap

Definition at line 66 of file MultiFieldMatrix.h.

Constructor & Destructor Documentation

MultiFieldMatrix::MultiFieldMatrix

(

)

Definition at line 16 of file MultiFieldMatrix.cpp.

References logCtor.

                                   :
  _matrices(),
  _coarseSizes(),
  _coarseGhostSizes(),
  _coarseSites(),
  _coarseToFineMappings(),
  _coarseConnectivities(),
  _coarseMatrices()
{
  logCtor();
}

MultiFieldMatrix::~MultiFieldMatrix ( )

virtual

Definition at line 28 of file MultiFieldMatrix.cpp.

References logDtor.

{
  logDtor();
}

Member Function Documentation

void MultiFieldMatrix::addMatrix ( const Index &  rowI, const Index &  colI, shared_ptr< Matrix >  m  )

inline

Definition at line 137 of file MultiFieldMatrix.h.

References _matrices.

Referenced by ElectricModel< T >::Impl::initChargeTransportLinearization(), FlowModel< T >::Impl::initContinuityLinearization(), PlateModel< T >::Impl::initDeformationLinearization(), StructureModel< T >::Impl::initDeformationLinearization(), ElectricModel< T >::Impl::initElectroStaticsLinearization(), KineticModel< T >::initKineticModelLinearization(), ThermalModel< T >::Impl::initLinearization(), FractureModel< T >::Impl::initLinearization(), VacancyModel< T >::Impl::initLinearization(), SpeciesModel< T >::Impl::initLinearization(), KeModel< T >::Impl::initLinearization(), KeModel< T >::Impl::initLinearizationk(), FlowModel< T >::Impl::initMomentumLinearization(), BatteryModel< T >::Impl::initPCLinearization(), PhononModel< T >::initPhononModelLinearization(), BatteryModel< T >::Impl::initPotentialLinearization(), BatteryModel< T >::Impl::initSpeciesLinearization(), and BatteryModel< T >::Impl::initThermalLinearization().

  {
    EntryIndex e(rowI,colI);
    _matrices[e] = m;
  }

void MultiFieldMatrix::computeResidual ( const IContainer &  xB, const IContainer &  bB, IContainer &  rB  ) const

virtual

Definition at line 402 of file MultiFieldMatrix.cpp.

References ArrayBase::copyPartial(), MultiField::getArrayIndex(), StorageSite::getCount(), MultiField::getLength(), getMatrix(), StorageSite::getSelfCount(), hasMatrix(), Matrix::multiplyAndAdd(), and ArrayBase::zeroPartial().

Referenced by AMG::cycle(), CG::solve(), JacobiSolver::solve(), ILU0Solver::solve(), BCGStab::solve(), SpikeSolver::solve(), DirectSolver::solve(), and AMG::solve().

{
  const MultiField& x = dynamic_cast<const MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
  MultiField& r = dynamic_cast<MultiField&>(rB);

  const int xLen = x.getLength();

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const ArrayBase& bI = b[rowIndex];
          const StorageSite& rowSite = *rowIndex.second;
          ArrayBase& rI = r[rowIndex];
          rI.copyPartial(bI,0,rowSite.getSelfCount());
          rI.zeroPartial(rowSite.getSelfCount(),rowSite.getCount());
          
          for(int j=0; j<xLen; j++)
          {
              const Index colIndex = x.getArrayIndex(j);
              if (hasMatrix(rowIndex,colIndex))
              {
                  const Matrix& mIJ = getMatrix(rowIndex,colIndex);
                  mIJ.multiplyAndAdd(rI,x[colIndex]);
              }
          }
      }
  }
#ifdef FVM_PARALLEL
  //r.sync();
#endif
}

void MultiFieldMatrix::correctSolution	(	const MultiField &	coarseIndex,
		MultiField &	fineSolutionField,
		MFRPtr	scaleField,
		const MultiField &	coarseSolutionField
	)

Definition at line 765 of file MultiFieldMatrix.cpp.

References _coarseSites, ArrayBase::correct(), MultiField::getArrayIndex(), and MultiField::getLength().

Referenced by AMG::cycle().

{
  const int xLen = fineSolutionField.getLength();

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = fineSolutionField.getArrayIndex(i);

      const StorageSite& coarseRowSite = *_coarseSites[rowIndex];

      const Index coarseRowIndex(rowIndex.first,&coarseRowSite);
      ArrayBase* scale(0);
      if (scaleField)
        scale = &((*scaleField)[*(coarseRowIndex.first)]);
      
      ArrayBase& fineSolution = dynamic_cast<ArrayBase&>(fineSolutionField[rowIndex]);
      fineSolution.correct(coarseSolutionField[coarseRowIndex],
                           coarseIndex[rowIndex],
                           scale,
                           rowIndex.second->getSelfCount());
  }
}

void MultiFieldMatrix::createCoarseConnectivity

(

MultiField &

coarseIndex

)

Definition at line 670 of file MultiFieldMatrix.cpp.

References _coarseConnectivities, _coarseSites, _coarseToFineMappings, Matrix::createCoarseConnectivity(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), and hasMatrix().

Referenced by LinearSystem::createCoarse().

{

  const int xLen = coarseIndex.getLength();

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = coarseIndex.getArrayIndex(i);

      const StorageSite& coarseRowSite = *_coarseSites[rowIndex];

      const CRConnectivity& coarseToFine = *_coarseToFineMappings[rowIndex];

      for(int j=0; j<xLen; j++)
      {
          const Index colIndex = coarseIndex.getArrayIndex(j);
          if (hasMatrix(rowIndex,colIndex))
          {
              Matrix& mIJ = getMatrix(rowIndex,colIndex);

              const StorageSite& coarseColSite = *_coarseSites[colIndex];

              shared_ptr<CRConnectivity> coarseConnectivity
                (mIJ.createCoarseConnectivity(coarseIndex[rowIndex],
                                              coarseToFine,
                                              coarseRowSite,coarseColSite));
              EntryIndex e(rowIndex,colIndex);
              _coarseConnectivities[e]=coarseConnectivity;
          }
      }
  }
}

void MultiFieldMatrix::createCoarseMatrices

(

MultiField &

coarseIndex

)

Definition at line 705 of file MultiFieldMatrix.cpp.

References _coarseConnectivities, _coarseMatrices, _coarseToFineMappings, Matrix::createCoarseMatrix(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), and hasMatrix().

Referenced by LinearSystem::createCoarse().

{
  const int xLen = coarseIndex.getLength();
  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = coarseIndex.getArrayIndex(i);

      const CRConnectivity& coarseToFine = *_coarseToFineMappings[rowIndex];

      for(int j=0; j<xLen; j++)
      {
          const Index colIndex = coarseIndex.getArrayIndex(j);
          if (hasMatrix(rowIndex,colIndex))
          {
              Matrix& mIJ = getMatrix(rowIndex,colIndex);
              EntryIndex e(rowIndex,colIndex);

              const CRConnectivity& coarseConnectivity = *_coarseConnectivities[e];

              shared_ptr<Matrix> coarseMatrix 
                (mIJ.createCoarseMatrix(coarseIndex[rowIndex],
                                        coarseToFine,
                                        coarseConnectivity));
              _coarseMatrices[e]=coarseMatrix;
          }
      }
  }
}

void MultiFieldMatrix::createCoarsening	(	MultiField &	coarseIndex,
		const int	groupSize,
		const double	weightRatioThreshold
	)

Definition at line 449 of file MultiFieldMatrix.cpp.

References _coarseSizes, Matrix::createCoarsening(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), and hasMatrix().

Referenced by LinearSystem::createCoarse().

{

  const int xLen = coarseIndex.getLength();


  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = coarseIndex.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          Matrix& mII = getMatrix(rowIndex,rowIndex);
         // cout << " proc_id = " << MPI::COMM_WORLD.Get_rank() << " coarse size = " << 
         //            mII.createCoarsening(coarseIndex[rowIndex], groupSize,weightRatioThreshold) <<endl;
          _coarseSizes[rowIndex] =
            mII.createCoarsening(coarseIndex[rowIndex],
                                 groupSize,weightRatioThreshold);
      }
  }
}

void MultiFieldMatrix::createCoarseToFineMapping

(

const MultiField &

coarseIndexField

)

Definition at line 627 of file MultiFieldMatrix.cpp.

References _coarseSites, _coarseToFineMappings, MultiField::getArrayIndex(), StorageSite::getCountLevel1(), and MultiField::getLength().

Referenced by LinearSystem::createCoarse().

{

  const int xLen = coarseIndexField.getLength();
  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = coarseIndexField.getArrayIndex(i);

      const StorageSite& fineSite = *rowIndex.second;

      const StorageSite& coarseSite = *_coarseSites[rowIndex];

      const int nFineRows = fineSite.getCountLevel1();
      

      const Array<int>& coarseIndex =
        dynamic_cast<const Array<int>& >(coarseIndexField[rowIndex]);

      shared_ptr<CRConnectivity> coarseToFine(new CRConnectivity(coarseSite,fineSite));

      coarseToFine->initCount();


      
      for(int nr=0; nr<nFineRows; nr++){
        if (coarseIndex[nr]>=0)
          coarseToFine->addCount(coarseIndex[nr],1);
      } 
          

      coarseToFine->finishCount();

      for(int nr=0; nr<nFineRows; nr++)
        if (coarseIndex[nr]>=0)
          coarseToFine->add(coarseIndex[nr],nr);

      coarseToFine->finishAdd();
      _coarseToFineMappings[rowIndex]=coarseToFine;
  }
}

MultiFieldMatrix::DEFINE_TYPENAME

(

"MultiFieldMatrix"

)

void MultiFieldMatrix::forwardGS ( IContainer &  xB, const IContainer &  bB, IContainer &  temp  ) const

virtual

Definition at line 126 of file MultiFieldMatrix.cpp.

References ArrayBase::copyPartial(), Matrix::forwardGS(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), StorageSite::getSelfCount(), hasMatrix(), Matrix::multiplyAndAdd(), MultiField::sync(), MultiField::syncGather(), and MultiField::syncScatter().

Referenced by AMG::doSweeps().

{
  MultiField& x = dynamic_cast<MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
   MultiField& temp = dynamic_cast<MultiField&>(tempB);

  const int xLen = x.getLength();
  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const ArrayBase& bI = b[rowIndex];
          ArrayBase& r = temp[rowIndex];
          const StorageSite& rowSite = *rowIndex.second;
          r.copyPartial(bI,0,rowSite.getSelfCount());
          
          for(int j=0; j<xLen; j++)
          {
              const Index colIndex = x.getArrayIndex(j);
              if ((rowIndex!=colIndex) && hasMatrix(rowIndex,colIndex))
              {
                  const Matrix& mIJ = getMatrix(rowIndex,colIndex);
                  mIJ.multiplyAndAdd(r,x[colIndex]);
              }
          }
          
          const Matrix& mII = getMatrix(rowIndex,rowIndex);
#ifndef FVM_PARALLEL
          x.syncGather(rowIndex);
#endif
          mII.forwardGS(x[rowIndex],r,r);
#ifndef FVM_PARALLEL
          x.syncScatter(rowIndex);
#endif
      }
  }
  x.sync();
}

int MultiFieldMatrix::getLocalSize

(

)

const

Definition at line 870 of file MultiFieldMatrix.cpp.

References _matrices.

Referenced by AMG::createCoarseLevels().

{
  int size=0;
  foreach(const MatrixMap::value_type& pos, _matrices)
  {
      EntryIndex k = pos.first;

      if (k.first == k.second)
      {
          size += k.first.second->getCount();
      }
  }

  return size;

}

Matrix & MultiFieldMatrix::getMatrix	(	const Index &	rowIndex,
		const Index &	colIndex
	)

Definition at line 42 of file MultiFieldMatrix.cpp.

References _matrices.

Referenced by StructureModel< T >::Impl::applyBC(), computeResidual(), correctMassFluxInterior(), createCoarseConnectivity(), createCoarseMatrices(), createCoarsening(), Underrelaxer< X, Diag, OffDiag >::discretize(), CollisionTermDiscretization< X, Diag, OffDiag >::discretize(), PhononCollisionDiscretization< X, Diag, OffDiag >::discretize(), LinearizeInterfaceJumpUnconnected< X, Diag, OffDiag >::discretize(), TrapBandTunnelingDiscretization< X, Diag, OffDiag >::discretize(), LinearizeInterfaceJump< X, Diag, OffDiag >::discretize(), DriftDiscretization< X, Diag, OffDiag >::discretize(), TimeDerivativeDiscretization< X, Diag, OffDiag >::discretize(), SourceDiscretizationforFracture< T, Diag, OffDiag >::discretize(), BatteryPCTimeDerivativeDiscretization< X, Diag, OffDiag >::discretize(), GenericIBDiscretization< X, Diag, OffDiag >::discretize(), LinearizeDielectric< X, Diag, OffDiag >::discretize(), PhononConvectionDiscretization< X, Diag, OffDiag >::discretize(), IbmDiscretization< X, Diag, OffDiag >::discretize(), TimeDerivativeDiscretization_Kmodel< X, Diag, OffDiag >::discretize(), ConvectionDiscretization< X, Diag, OffDiag >::discretize(), BatteryLinearizePotentialInterface< X, Diag, OffDiag >::discretize(), BatteryLinearizeThermalInterface< X, Diag, OffDiag >::discretize(), LinearizePotentialInterface< X, Diag, OffDiag >::discretize(), BatteryPCLinearizeInterface_BV< X, Diag, OffDiag, otherMeshDiag >::discretize(), TimeDerivativeStructureDiscretization< X, Diag, OffDiag >::discretize(), BatteryLinearizeSpeciesInterface< X, Diag, OffDiag >::discretize(), GenericKineticIBDiscretization< X, Diag, OffDiag >::discretize(), LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize(), ConvectionDiscretization_Kmodel< X, Diag, OffDiag >::discretize(), MomentumPressureGradientDiscretization< X >::discretize(), TimeDerivativePlateDiscretization< X, Diag, OffDiag >::discretize(), ElecDiffusionDiscretization< X, Diag, OffDiag >::discretize(), CaptureDiscretization< X, Diag, OffDiag >::discretize(), DiffusionDiscretization< X, Diag, OffDiag >::discretize(), EmissionDiscretization< X, Diag, OffDiag >::discretize(), SourceDiscretizationdissi< X, Diag, OffDiag >::discretize(), TunnelingDiscretization< X, Diag, OffDiag >::discretize(), BatteryPCDiffusionDiscretization< X, Diag, OffDiag >::discretize(), WallDiscretization< X, Diag, OffDiag >::discretize(), PlateSourceDiscretization< T, Diag, OffDiag >::discretize(), StructureSourceDiscretization< T, Diag, OffDiag >::discretizeFaces(), StructurePlasticDiscretization< T, Diag, OffDiag >::discretizeFaces(), PlateSourceDiscretization< T, Diag, OffDiag >::discretizeFaces(), discretizeMassFluxInterior(), FlowModel< T >::Impl::dumpContinuityMatrix(), ThermalModel< T >::Impl::dumpMatrix(), VacancyModel< T >::Impl::dumpMatrix(), PlateModel< T >::Impl::dumpMatrix(), fixedFluxContinuityBC(), fixedPressureContinuityBC(), fixedPressureMomentumBC(), forwardGS(), iluSolve(), LinearSystem::initSolve(), Jacobi(), FlowModel< T >::Impl::linearizeContinuity(), PlateModel< T >::Impl::linearizeDeformation(), KeModel< T >::Impl::linearizedissipation(), multiply(), multiplyAndAdd(), BatteryModel< T >::Impl::printMatrixElementsOnFace(), quadProduct(), reverseGS(), FlowModel< T >::Impl::setDirichlet(), DirectSolver::solve(), solveBoundary(), FlowModel< T >::Impl::solveMomentum(), and spikeSolve().

{
  EntryIndex i(row,col);
  return *_matrices.find(i)->second;
}

const Matrix & MultiFieldMatrix::getMatrix	(	const Index &	rowIndex,
		const Index &	colIndex
	)		const

Definition at line 50 of file MultiFieldMatrix.cpp.

References _matrices.

{
  EntryIndex i(row,col);
  return *_matrices.find(i)->second;
}

int MultiFieldMatrix::getSize

(

)

const

Definition at line 845 of file MultiFieldMatrix.cpp.

References _matrices.

{
  int size=0;
  foreach(const MatrixMap::value_type& pos, _matrices)
  {
      EntryIndex k = pos.first;

      if (k.first == k.second)
      {
          size += k.first.second->getCount();
      }
  }

#ifdef FVM_PARALLEL
         int count = 1;
         MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &size, count, MPI::INT, MPI::MIN);

#endif

  return size;

}

bool MultiFieldMatrix::hasMatrix	(	const Index &	rowIndex,
		const Index &	colIndex
	)		const

Definition at line 34 of file MultiFieldMatrix.cpp.

References _matrices.

Referenced by computeResidual(), correctMassFluxInterior(), LinearSystem::createCoarse(), createCoarseConnectivity(), createCoarseMatrices(), createCoarsening(), MomentumPressureGradientDiscretization< X >::discretize(), discretizeMassFluxInterior(), fixedFluxContinuityBC(), forwardGS(), iluSolve(), LinearSystem::initSolve(), Jacobi(), multiply(), multiplyAndAdd(), FlowModel< T >::Impl::postContinuitySolve(), quadProduct(), reverseGS(), FlowModel< T >::Impl::setDirichlet(), solveBoundary(), and spikeSolve().

{
  EntryIndex i(row,col);
  return _matrices.find(i) != _matrices.end();
}

void MultiFieldMatrix::iluSolve ( IContainer &  xB, const IContainer &  bB, IContainer &  tempB  ) const

virtual

Definition at line 276 of file MultiFieldMatrix.cpp.

References MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), hasMatrix(), Matrix::iluSolve(), and MultiField::sync().

Referenced by ILU0Solver::doSweeps().

{
  MultiField& x = dynamic_cast<MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
  MultiField& temp = dynamic_cast<MultiField&>(tempB);
  
  
  const int xLen = x.getLength();
  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const Matrix& mII = getMatrix(rowIndex,rowIndex);
          mII.iluSolve(x[rowIndex],b[rowIndex],temp);
      }
  }
  x.sync();
}

void MultiFieldMatrix::initAssembly

(

)

Definition at line 58 of file MultiFieldMatrix.cpp.

References _matrices.

Referenced by LinearSystem::initAssembly().

{
  foreach(MatrixMap::value_type& pos, _matrices)
    pos.second->initAssembly();
}

void MultiFieldMatrix::injectResidual	(	const MultiField &	coarseIndex,
		const MultiField &	fineResidualField,
		MultiField &	coarseBField
	)

Definition at line 737 of file MultiFieldMatrix.cpp.

References _coarseSites, MultiField::getArrayIndex(), MultiField::getLength(), and ArrayBase::inject().

Referenced by AMG::cycle().

{
  const int xLen = fineResidualField.getLength();

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = fineResidualField.getArrayIndex(i);

      const StorageSite& coarseRowSite = *_coarseSites[rowIndex];

      const Index coarseRowIndex(rowIndex.first,&coarseRowSite);
      
      const ArrayBase& fineResidual =
        dynamic_cast<const ArrayBase&>(fineResidualField[rowIndex]);

      const ArrayBase& fineToCoarse =
        dynamic_cast<const ArrayBase&>(coarseIndex[rowIndex]);

      ArrayBase& coarseB =  dynamic_cast<ArrayBase&>(coarseBField[coarseRowIndex]);

      fineResidual.inject(coarseB,fineToCoarse,rowIndex.second->getSelfCount());
  }
}

void MultiFieldMatrix::Jacobi ( IContainer &  xB, const IContainer &  bB, IContainer &  tempB  ) const

virtual

Definition at line 168 of file MultiFieldMatrix.cpp.

References ArrayBase::copyPartial(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), StorageSite::getSelfCount(), hasMatrix(), Matrix::Jacobi(), Matrix::multiplyAndAdd(), MultiField::newClone(), and MultiField::sync().

Referenced by JacobiSolver::doSweeps(), and AMG::doSweeps().

{
  MultiField& x = dynamic_cast<MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
  MultiField& temp = dynamic_cast<MultiField&>(tempB);
  
  shared_ptr<MultiField> xnew = dynamic_pointer_cast<MultiField>(x.newClone());
  
  const int xLen = x.getLength();
  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const ArrayBase& bI = b[rowIndex];
          ArrayBase& r = temp[rowIndex];
          const StorageSite& rowSite = *rowIndex.second;
          r.copyPartial(bI,0,rowSite.getSelfCount());
          
          for(int j=0; j<xLen; j++)
          {
              const Index colIndex = x.getArrayIndex(j);
              if ((rowIndex!=colIndex) && hasMatrix(rowIndex,colIndex))
              {
                  const Matrix& mIJ = getMatrix(rowIndex,colIndex);
                  mIJ.multiplyAndAdd(r,x[colIndex]);
              }
          }
          
          const Matrix& mII = getMatrix(rowIndex,rowIndex);
          mII.Jacobi((*xnew)[rowIndex],x[rowIndex],r);

      }
  }

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const ArrayBase& xnewI = (*xnew)[rowIndex];
          ArrayBase& xI = x[rowIndex];
          const StorageSite& rowSite = *rowIndex.second;
          xI.copyPartial(xnewI,0,rowSite.getSelfCount());
      }
  }
  x.sync();
}

void MultiFieldMatrix::multiply ( IContainer &  yB, const IContainer &  xB  ) const

virtual

Definition at line 65 of file MultiFieldMatrix.cpp.

References MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), hasMatrix(), Matrix::multiplyAndAdd(), MultiField::sync(), and IContainer::zero().

Referenced by StructureModel< T >::Impl::advanceExplicit(), BCGStab::solve(), and CG::solve().

{
  MultiField& y = dynamic_cast<MultiField&>(yB);
  const MultiField& x = dynamic_cast<const MultiField&>(xB);

  const int yLen = y.getLength();
  const int xLen = x.getLength();

  //#pragma omp parallel for
  for(int i=0; i<yLen; i++)
  {
      const Index rowIndex = y.getArrayIndex(i);
      ArrayBase& yI = y[rowIndex];
      yI.zero();
      for(int j=0; j<xLen; j++)
      {
          const Index colIndex = x.getArrayIndex(j);
          if (hasMatrix(rowIndex,colIndex))
          {
              const Matrix& mIJ = getMatrix(rowIndex,colIndex);
              mIJ.multiplyAndAdd(yI,x[colIndex]);
          }
      }
  }
#ifdef FVM_PARALLEL
  y.sync();
#endif
}

void MultiFieldMatrix::multiplyAndAdd ( IContainer &  yB, const IContainer &  xB  ) const

virtual

Definition at line 95 of file MultiFieldMatrix.cpp.

References MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), hasMatrix(), Matrix::multiplyAndAdd(), and MultiField::sync().

{
  MultiField& y = dynamic_cast<MultiField&>(yB);
  const MultiField& x = dynamic_cast<const MultiField&>(xB);

  const int yLen = y.getLength();
  const int xLen = x.getLength();

  //#pragma omp parallel for
  for(int i=0; i<yLen; i++)
  {
      const Index rowIndex = y.getArrayIndex(i);
      ArrayBase& yI = y[rowIndex];
      for(int j=0; j<xLen; j++)
      {
          const Index colIndex = x.getArrayIndex(j);
          if (hasMatrix(rowIndex,colIndex))
          {
              const Matrix& mIJ = getMatrix(rowIndex,colIndex);
              mIJ.multiplyAndAdd(yI,x[colIndex]);
          }
      }
  }

#ifdef FVM_PARALLEL
  y.sync();
#endif

}

MFRPtr MultiFieldMatrix::quadProduct

(

const MultiField &

x

)

const

Definition at line 888 of file MultiFieldMatrix.cpp.

References MultiField::addArray(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), hasMatrix(), Matrix::quadProduct(), and MultiField::reduceSum().

Referenced by AMG::cycle().

{
  const int xLen = x.getLength();

  MultiField *p = new MultiField();
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      for(int j=0; j<xLen; j++)
      {
          const Index colIndex = x.getArrayIndex(j);
          if (hasMatrix(rowIndex,colIndex))
          {
              const Matrix& mIJ = getMatrix(rowIndex,colIndex);
              p->addArray(rowIndex,
                          mIJ.quadProduct(x[colIndex]));
          }
      }
  }

  MFRPtr r = p->reduceSum();
  delete p;
  return r;
}

void MultiFieldMatrix::removeMatrix	(	const Index &	rowIndex,
		const Index &	colIndex
	)

Definition at line 440 of file MultiFieldMatrix.cpp.

References _matrices.

{
  EntryIndex e(row,col);

  _matrices.erase(e);
}

void MultiFieldMatrix::reverseGS ( IContainer &  xB, const IContainer &  bB, IContainer &  temp  ) const

virtual

Definition at line 358 of file MultiFieldMatrix.cpp.

References ArrayBase::copyPartial(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), StorageSite::getSelfCount(), hasMatrix(), Matrix::multiplyAndAdd(), Matrix::reverseGS(), MultiField::sync(), MultiField::syncGather(), and MultiField::syncScatter().

Referenced by AMG::doSweeps().

{
  MultiField& x = dynamic_cast<MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
  MultiField& temp = dynamic_cast<MultiField&>(tempB);

  const int xLen = x.getLength();

  for(int i=xLen-1; i>=0; i--)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const ArrayBase& bI = b[rowIndex];
          ArrayBase& r = temp[rowIndex];
          const StorageSite& rowSite = *rowIndex.second;
          r.copyPartial(bI,0,rowSite.getSelfCount());
                    
          for(int j=0; j<xLen; j++)
          {
              const Index colIndex = x.getArrayIndex(j);
              if ((rowIndex!=colIndex) && hasMatrix(rowIndex,colIndex))
              {
                  const Matrix& mIJ = getMatrix(rowIndex,colIndex);
                  mIJ.multiplyAndAdd(r,x[colIndex]);
              }
          }

          const Matrix& mII = getMatrix(rowIndex,rowIndex);
#ifndef FVM_PARALLEL
          x.syncGather(rowIndex);
#endif
          mII.reverseGS(x[rowIndex],r,r);
#ifndef FVM_PARALLEL
         x.syncScatter(rowIndex);
#endif
      }
  }

  x.sync();

}

void MultiFieldMatrix::solveBoundary ( IContainer &  xB, const IContainer &  bB, IContainer &  temp  ) const

virtual

Definition at line 321 of file MultiFieldMatrix.cpp.

References ArrayBase::copyFrom(), MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), hasMatrix(), Matrix::multiplyAndAdd(), and Matrix::solveBoundary().

Referenced by LinearSystem::postSolve().

{
  MultiField& x = dynamic_cast<MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
   MultiField& temp = dynamic_cast<MultiField&>(tempB);

  const int xLen = x.getLength();

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const ArrayBase& bI = b[rowIndex];
          ArrayBase& r = temp[rowIndex];
          r.copyFrom(bI);
          
          for(int j=0; j<xLen; j++)
          {
              const Index colIndex = x.getArrayIndex(j);
              if ((rowIndex!=colIndex) && hasMatrix(rowIndex,colIndex))
              {
                  const Matrix& mIJ = getMatrix(rowIndex,colIndex);
                  mIJ.multiplyAndAdd(r,x[colIndex]);
              }
          }

          const Matrix& mII = getMatrix(rowIndex,rowIndex);
          mII.solveBoundary(x[rowIndex],r,r);
      }
  }
}

void MultiFieldMatrix::spikeSolve ( IContainer &  xB, const IContainer &  bB, IContainer &  tempB, const SpikeStorage &  spike_storage  ) const

virtual

Definition at line 299 of file MultiFieldMatrix.cpp.

References MultiField::getArrayIndex(), MultiField::getLength(), getMatrix(), hasMatrix(), Matrix::spikeSolve(), and MultiField::sync().

Referenced by SpikeSolver::doSweeps().

{
  MultiField& x = dynamic_cast<MultiField&>(xB);
  const MultiField& b = dynamic_cast<const MultiField&>(bB);
  MultiField& temp = dynamic_cast<MultiField&>(tempB);

  const int xLen = x.getLength();
  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = x.getArrayIndex(i);
      if (hasMatrix(rowIndex,rowIndex))
      {
          const Matrix& mII = getMatrix(rowIndex,rowIndex);
          mII.spikeSolve(x[rowIndex], b[rowIndex], temp, spike_storage);
      }
  }
  x.sync();
}

void MultiFieldMatrix::syncGhostCoarsening

(

MultiField &

coarseIndexField

)

Definition at line 476 of file MultiFieldMatrix.cpp.

References _coarseGatherMaps, _coarseGhostSizes, _coarseScatterMaps, _coarseSizes, MultiField::getArrayIndex(), UnorderedSet::getData(), StorageSite::getGatherMap(), StorageSite::getGatherMapLevel1(), MultiField::getLength(), Array< T >::getLength(), StorageSite::getScatterMap(), StorageSite::getScatterMapLevel1(), UnorderedSet::insert(), and UnorderedSet::size().

Referenced by LinearSystem::createCoarse().

{
  const int xLen = coarseIndexField.getLength();

  //#pragma omp parallel for
  for(int i=0; i<xLen; i++)
  {
      const Index rowIndex = coarseIndexField.getArrayIndex(i);

      Array<int>&  coarseIndex = dynamic_cast<Array<int>& >(coarseIndexField[rowIndex]);


      int coarseGhostSize=0;
      const int coarseSize = _coarseSizes.find(rowIndex)->second;

      const StorageSite& site = *rowIndex.second;

      const StorageSite::GatherMap&  gatherMap  = site.getGatherMap();
      const StorageSite::GatherMap&  gatherMapLevel1  = site.getGatherMapLevel1();

      // collect all the toIndices arrays for each storage site from
      // both gatherMap and gatherMapLevel1
      
      typedef map<const StorageSite*, vector<const Array<int>* > > IndicesMap;
      IndicesMap toIndicesMap;

      foreach(const StorageSite::GatherMap::value_type pos, gatherMap)
      {
          const StorageSite& oSite = *pos.first;
          const Array<int>& toIndices = *pos.second;

          toIndicesMap[&oSite].push_back(&toIndices);
      }
      
      foreach(const StorageSite::GatherMap::value_type pos, gatherMapLevel1)
      {
          const StorageSite& oSite = *pos.first;
          const Array<int>& toIndices = *pos.second;

          toIndicesMap[&oSite].push_back(&toIndices);
      }
      
      foreach(IndicesMap::value_type pos, toIndicesMap)
        {
          const StorageSite& oSite = *pos.first;
          const vector<const Array<int>* > toIndicesArrays = pos.second;

          
          map<int,int> otherToMyMapping;
          UnorderedSet  gatherSet;

          foreach(const Array<int>* toIndicesPtr, toIndicesArrays)
          {
              const Array<int>& toIndices = *toIndicesPtr;
              const int nGhostRows = toIndices.getLength();
              for(int ng=0; ng<nGhostRows; ng++)
              {
                  const int fineIndex = toIndices[ng];
                  const int coarseOtherIndex = coarseIndex[fineIndex];
                  
                  if (coarseOtherIndex < 0)
                    continue;
                  
                  if (otherToMyMapping.find(coarseOtherIndex) !=
                      otherToMyMapping.end())
                  {
                      coarseIndex[fineIndex] = otherToMyMapping[coarseOtherIndex];
                  }
                  else
                  {
                      coarseIndex[fineIndex] = coarseGhostSize+coarseSize;
                      otherToMyMapping[coarseOtherIndex] = coarseIndex[fineIndex];
                      gatherSet.insert( coarseIndex[fineIndex] );
                      coarseGhostSize++;
                  }
              }
          }

          const int coarseMappersSize = otherToMyMapping.size();

          shared_ptr<Array<int> > coarseToIndices(new Array<int>(coarseMappersSize));

          for(int n = 0; n < gatherSet.size(); n++)
          {
              (*coarseToIndices)[n]   = gatherSet.getData().at(n);
          }

          SSPair sskey(&site,&oSite);
          _coarseGatherMaps [sskey] = coarseToIndices;

        }

      const StorageSite::ScatterMap& scatterMap = site.getScatterMap();
      const StorageSite::ScatterMap& scatterMapLevel1 = site.getScatterMapLevel1();
      
      IndicesMap fromIndicesMap;

      foreach(const StorageSite::GatherMap::value_type pos, scatterMap)
      {
          const StorageSite& oSite = *pos.first;
          const Array<int>& fromIndices = *pos.second;

          fromIndicesMap[&oSite].push_back(&fromIndices);
      }
      
      foreach(const StorageSite::GatherMap::value_type pos, scatterMapLevel1)
      {
          const StorageSite& oSite = *pos.first;
          const Array<int>& fromIndices = *pos.second;

          fromIndicesMap[&oSite].push_back(&fromIndices);
      }

      foreach(IndicesMap::value_type pos, fromIndicesMap)
      {
          const StorageSite& oSite = *pos.first;
          const vector<const Array<int>* > fromIndicesArrays = pos.second;

          UnorderedSet  scatterSet;

          foreach(const Array<int>* fromIndicesPtr, fromIndicesArrays)
          {
              const Array<int>& fromIndices = *fromIndicesPtr;
              const int nGhostRows = fromIndices.getLength();
              for(int ng=0; ng<nGhostRows; ng++)
              {
                  const int fineIndex = fromIndices[ng];
                  const int coarseOtherIndex = coarseIndex[fineIndex];
                  if (coarseOtherIndex >= 0)
                    scatterSet.insert( coarseOtherIndex );
              }

          }

          const int coarseMappersSize = scatterSet.size();
          
          shared_ptr<Array<int> > coarseFromIndices(new Array<int>(coarseMappersSize));
          
          for(int n = 0; n < scatterSet.size(); n++ ) {
              (*coarseFromIndices)[n] = scatterSet.getData().at(n);
          }
          
          SSPair sskey(&site,&oSite);
          _coarseScatterMaps[sskey] = coarseFromIndices;
          
      }
      _coarseGhostSizes[rowIndex]=coarseGhostSize;
  }
}

void MultiFieldMatrix::transpose

(

)

Definition at line 220 of file MultiFieldMatrix.cpp.

References _matrices.

{
  // keep the transposed offdiagonals in this map since we can't modify _matrices while looping over it
  MatrixMap transposedMap;
  
  foreach(const MatrixMap::value_type& pos, _matrices)
  {
      EntryIndex k = pos.first;

      if (k.first == k.second)
      {
          pos.second->transpose();
      }
      else
      {
          EntryIndex kt(k.second,k.first);
          if ( (transposedMap.count(k) + transposedMap.count(kt)) == 0)
          {
              shared_ptr<Matrix> mIJ = pos.second;
              mIJ->transpose();

              if (_matrices.find(kt) != _matrices.end())
              {
                  shared_ptr<Matrix> mJI = _matrices[kt];
                  mJI->transpose();

                  transposedMap[k] = mJI;
                  transposedMap[kt] = mIJ;
              }
              else
              {
                  transposedMap[kt] = mIJ;
                  // store a null matrix so in the loop below we can
                  // erase the entry at this location in _matrices
                  transposedMap[k] = shared_ptr<Matrix>();
              }
          }
      }
  }

  // set the matrices from transposedMap in the transposed locations in _matrices
  foreach(const MatrixMap::value_type& pos, transposedMap)
  {
      EntryIndex k = pos.first;

      if (pos.second)
        _matrices[k] = pos.second;
      else
        _matrices.erase(k);
  }
}

Friends And Related Function Documentation

friend class LinearSystem

friend

Definition at line 53 of file MultiFieldMatrix.h.

friend class LinearSystemMerger

friend

Definition at line 54 of file MultiFieldMatrix.h.

Member Data Documentation

CoarseConnectivitiesMap MultiFieldMatrix::_coarseConnectivities

private

Definition at line 156 of file MultiFieldMatrix.h.

Referenced by createCoarseConnectivity(), and createCoarseMatrices().

MatrixMappersMap MultiFieldMatrix::_coarseGatherMaps

private

Definition at line 152 of file MultiFieldMatrix.h.

Referenced by LinearSystem::createCoarse(), and syncGhostCoarsening().

MatrixSizeMap MultiFieldMatrix::_coarseGhostSizes

private

Definition at line 150 of file MultiFieldMatrix.h.

Referenced by LinearSystem::createCoarse(), and syncGhostCoarsening().

MatrixMap MultiFieldMatrix::_coarseMatrices

private

Definition at line 157 of file MultiFieldMatrix.h.

Referenced by LinearSystem::createCoarse(), and createCoarseMatrices().

MatrixMappersMap MultiFieldMatrix::_coarseScatterMaps

private

Definition at line 151 of file MultiFieldMatrix.h.

Referenced by LinearSystem::createCoarse(), and syncGhostCoarsening().

StorageSiteMap MultiFieldMatrix::_coarseSites

private

Definition at line 154 of file MultiFieldMatrix.h.

Referenced by correctSolution(), LinearSystem::createCoarse(), createCoarseConnectivity(), createCoarseToFineMapping(), and injectResidual().

MatrixSizeMap MultiFieldMatrix::_coarseSizes

private

Definition at line 149 of file MultiFieldMatrix.h.

Referenced by LinearSystem::createCoarse(), createCoarsening(), and syncGhostCoarsening().

CoarseToFineMappingMap MultiFieldMatrix::_coarseToFineMappings

private

Definition at line 155 of file MultiFieldMatrix.h.

Referenced by LinearSystem::createCoarse(), createCoarseConnectivity(), createCoarseMatrices(), and createCoarseToFineMapping().

MatrixMap MultiFieldMatrix::_matrices

private

Definition at line 148 of file MultiFieldMatrix.h.

Referenced by addMatrix(), getLocalSize(), getMatrix(), getSize(), hasMatrix(), initAssembly(), removeMatrix(), and transpose().

---

The documentation for this class was generated from the following files:

• MultiFieldMatrix.h
• MultiFieldMatrix.cpp