Commit fbab6ab5 authored by Radu Popescu's avatar Radu Popescu
Browse files

Removing obsolete LinearEpetraOperator files

parent 5562eb50
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file
@brief Some linear operators which inherit from Epetra_Operator
@author Umberto Villa <uvilla.@emory.edu>
@date 03-09-2010
In this file we define a Lifev::Operators:LinearOperator which inherit from Epetra_Operator.
LinearOperator adds to the interface of Epetra_Operator a method "Apply" and "ApplyInverse"
which work with LifeV::VectorEpetra instead of Epetra_MultiVector.
LinearOperator is an abstract class with virtual public interface.
IdentityOperator and NullOperator are defined in this file and both inherit from LinearOperator.
*/
#ifndef LINEAREPETRAOPERATOR_H
#define LINEAREPETRAOPERATOR_H 1
#include <life/lifecore/LifeV.hpp>
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#pragma GCC diagnostic ignored "-Wextra"
#include <Epetra_Comm.h>
#include <Epetra_Map.h>
#include <Epetra_Operator.h>
#include <Epetra_MultiVector.h>
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"
#pragma GCC diagnostic warning "-Wextra"
#include <life/lifearray/VectorEpetra.hpp>
namespace LifeV
{
namespace Operators
{
//! @class LinearOperator
/*! @brief Abstract class which defines the interface of a Linear Operator.
*
* LinearOperator is an abstract which inherits from Epetra_Operator.
* LinearOperator should be the base class for all the LifeV class which implements a linear operator.
*
* LinearOperator ensures perfect compatibility with all the Trilinos solvers,
* plus it supports directly the LifeV::VectorEpetra data.
*
* Two concrete methods are implemented in LinearOperator
* int apply(const VectorEpetra &X, VectorEpetra & Y) const ;
* int applyInverse(const VectorEpetra &X, VectorEpetra &Y) const.
*
*
* Such methods extract a raw Epetra_MultiVector from the VectorEpetra and then call the virtual methods
* Int Apply(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
* or
* Int ApplyInverse(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
* respectively.
*
*
*/
class LinearOperator : public Epetra_Operator
{
public:
//! @name Destructor
//@{
//! Destructor
virtual ~LinearOperator() {};
//@}
//! @name Attribute set methods
//@{
//! If set true, transpose of this operator will be applied.
/*! This flag allows the transpose of the given operator to be used implicitly. Setting this flag
affects only the Apply() and ApplyInverse() methods. If the implementation of this interface
does not support transpose use, this method should return a value of -1.
\param In
UseTranspose -If true, multiply by the transpose of operator, otherwise just use operator.
\return Integer error code, set to 0 if successful. Set to -1 if this implementation does not support transpose.
*/
virtual int SetUseTranspose(bool UseTranspose) = 0;
//@}
//! @name Mathematical functions
//@{
//! Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
/*!
\param In
X - A Epetra_MultiVector of dimension NumVectors to multiply with matrix.
\param Out
Y -A Epetra_MultiVector of dimension NumVectors containing result.
\return Integer error code, set to 0 if successful.
*/
virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const = 0;
//! Returns the result of a LinearOperator applied to a VectorEpetra X in Y.
/*!
\param In
X - A VectorEpetra to multiply with matrix.
\param Out
Y -A VectorEpetra containing result.
\return Integer error code, set to 0 if successful.
*/
inline int apply(const VectorEpetra & X, VectorEpetra & Y) const
{
return Apply(X.epetraVector(), Y.epetraVector());
}
//! Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
/*!
\param In
X - A Epetra_MultiVector of dimension NumVectors to solve for.
\param Out
Y -A Epetra_MultiVector of dimension NumVectors containing result.
\return Integer error code, set to 0 if successful.
\warning In order to work with AztecOO, any implementation of this method must
support the case where X and Y are the same object.
*/
virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const = 0;
//! Returns the result of a LinearOperator inverse applied to an VectorEpetra X in Y.
/*!
\param In
X - A VectorEpetra to solve for.
\param Out
Y -A VectorEpetra containing result.
\return Integer error code, set to 0 if successful.
\warning In order to work with AztecOO, any implementation of this method must
support the case where X and Y are the same object.
*/
inline int applyInverse(const VectorEpetra & X, VectorEpetra Y)
{
return ApplyInverse(X.epetraVector(), Y.epetraVector());
}
//! Returns the infinity norm of the global matrix.
/* Returns the quantity \f$ \| A \|_\infty\f$ such that
\f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
\warning This method must not be called unless HasNormInf() returns true.
*/
virtual double NormInf() const = 0;
//@}
//! @name Attribute access functions
//@{
//! Returns a character string describing the operator
virtual const char * Label() const = 0;
//! Returns the current UseTranspose setting.
virtual bool UseTranspose() const = 0;
//! Returns true if the \e this object can provide an approximate Inf-norm, false otherwise.
virtual bool HasNormInf() const = 0;
//! Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual const Epetra_Comm & Comm() const = 0;
//! Returns the Epetra_Map object associated with the domain of this operator.
virtual const Epetra_Map & OperatorDomainMap() const = 0;
//! Returns the Epetra_Map object associated with the range of this operator.
virtual const Epetra_Map & OperatorRangeMap() const = 0;
//@}
};
//! @class IdentityOperator
/*! @brief Identity operator x = I*x. */
class IdentityOperator : public LinearOperator
{
public:
IdentityOperator():M_name("identity"), M_useTranspose(false) {};
void setUp(const boost::shared_ptr<Epetra_Map> & map) {M_map = map;}
int SetUseTranspose(bool useTranspose) {M_useTranspose = useTranspose; return 0; }
int Apply(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const {Y = X; return 0;}
int ApplyInverse(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const {Y=X; return 0;}
double NormInf() const {return 1.0;}
const char * Label() const {return M_name.c_str();}
bool UseTranspose() const {return M_useTranspose;}
bool HasNormInf() const {return true;}
const Epetra_Comm & Comm() const {return M_map->Comm();}
const Epetra_Map & OperatorDomainMap() const {return *M_map;}
const Epetra_Map & OperatorRangeMap() const {return *M_map;}
private:
std::string M_name;
boost::shared_ptr<Epetra_Map> M_map;
bool M_useTranspose;
};
//! @class NullOperator
/*! @brief Null operator 0 = Z*x. */
class NullOperator : public LinearOperator
{
public:
NullOperator():M_name("Null Operator"), M_useTranspose(false) {};
void setUp(const boost::shared_ptr<Epetra_Map> & domainMap,
const boost::shared_ptr<Epetra_Map> & rangeMap)
{
M_domainMap = domainMap;
M_rangeMap = rangeMap;
}
int SetUseTranspose(bool useTranspose) {M_useTranspose = useTranspose; return 0; }
int Apply(const Epetra_MultiVector & /*X*/, Epetra_MultiVector & Y) const {Y.PutScalar(0.0); return 0;}
int ApplyInverse(const Epetra_MultiVector & /*X*/, Epetra_MultiVector & Y) const
{Y.PutScalar(std::numeric_limits<Real>::quiet_NaN( )); return -1;}
double NormInf() const {return 0.0;}
const char * Label() const {return M_name.c_str();}
bool UseTranspose() const {return M_useTranspose;}
bool HasNormInf() const {return true;}
const Epetra_Comm & Comm() const {return M_rangeMap->Comm();}
const Epetra_Map & OperatorDomainMap() const {return *M_domainMap;}
const Epetra_Map & OperatorRangeMap() const {return *M_rangeMap;}
private:
std::string M_name;
boost::shared_ptr<Epetra_Map> M_domainMap;
boost::shared_ptr<Epetra_Map> M_rangeMap;
bool M_useTranspose;
};
} /*end namespace Operators*/
} /*end namespace */
#endif /* LINEAREPETRAOPERATOR_H */
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/*!
@file
@brief A class to manage block structured operators
@author Umberto Villa <uvilla@emory.edu>
@date 20-09-2010
LifeV::Operators::BlockOperator manages block structured operators in a monolithic approach.
BlockOperator inherit from LifeV::Operators::LinearOperator and assumes that each block is also a
LifeV::Operators::LinearOperator.
*/
#include <life/lifealg/LinearEpetraOperatorBlock.hpp>
namespace LifeV
{
namespace Operators
{
BlockOperator::BlockOperator():
M_name("BlockOperator"),
M_useTranspose(false),
M_structure(NoStructure)
{};
void BlockOperator::setUp(UInt nBlocks,
const std::vector<boost::shared_ptr<Epetra_Map> > domainMap,
const boost::shared_ptr<Epetra_Comm> & comm)
{
M_comm = comm;
M_nBlockRows = nBlocks;
M_nBlockCols = nBlocks;
M_domainBlockMaps = domainMap;
M_domainBlockMapsShift.resize(nBlocks);
for (UInt iblock=0; iblock < nBlocks; ++iblock)
M_domainBlockMapsShift[iblock].reset(new Epetra_Map(*M_domainBlockMaps[iblock]));
buildImporter(nBlocks, M_domainBlockMapsShift, M_domainMap,
M_block2monoDomain, M_mono2blockDomain);
//Since the operator is square I'll just copy the pointers not the content of the pointers
M_rangeBlockMaps = M_domainBlockMaps;
M_rangeBlockMapsShift = M_domainBlockMapsShift;
M_rangeMap = M_domainMap;
M_block2monoRange = M_block2monoDomain;
M_mono2blockRange = M_mono2blockDomain;
M_oper.resize(nBlocks, nBlocks);
}
//! SetUp for a "rectangular operator"
void BlockOperator::setUp(UInt nRowBlocks, UInt nColBlocks,
const std::vector<boost::shared_ptr<Epetra_Map> > & domainMap,
const std::vector<boost::shared_ptr<Epetra_Map> > & rangeMap,
const boost::shared_ptr<Epetra_Comm> & comm)
{
M_comm = comm;
M_nBlockRows = nRowBlocks;
M_nBlockCols = nColBlocks;
M_rangeBlockMaps = rangeMap;
M_domainBlockMaps = domainMap;
M_rangeBlockMapsShift.resize(M_nBlockRows);
for (UInt iblock=0; iblock < M_nBlockRows; ++iblock)
M_rangeBlockMapsShift[iblock].reset(new Epetra_Map(*M_rangeBlockMaps[iblock]));
M_domainBlockMapsShift.resize(M_nBlockCols);
for (UInt jblock=0; jblock < M_nBlockCols; ++jblock)
M_domainBlockMapsShift[jblock].reset(new Epetra_Map(*M_domainBlockMaps[jblock]));
buildImporter(M_nBlockRows, M_rangeBlockMapsShift, M_rangeMap,
M_block2monoRange, M_mono2blockRange);
buildImporter(M_nBlockCols, M_domainBlockMapsShift, M_domainMap,
M_block2monoDomain, M_mono2blockDomain);
M_oper.resize(M_nBlockRows, M_nBlockCols);
}
//! SetUp when the operator is given like a boost::matrix
void BlockOperator::setUp(const BlockOper & blockOper, const boost::shared_ptr<Epetra_Comm> & comm)
{
M_nBlockRows = blockOper.size1();
M_nBlockCols = blockOper.size2();
M_rangeBlockMaps.resize(M_nBlockRows);
M_domainBlockMaps.resize(M_nBlockCols);
for (UInt iblock=0; iblock < M_nBlockRows; ++iblock)
M_rangeBlockMaps[iblock].reset(new Epetra_Map(blockOper(iblock,0)->OperatorRangeMap() ));
for (UInt jblock=0; jblock < M_nBlockCols; ++jblock)
M_domainBlockMaps[jblock].reset(new Epetra_Map(blockOper(0,jblock)->OperatorDomainMap() ));
setUp(M_nBlockRows, M_nBlockCols, M_domainBlockMaps, M_rangeBlockMaps, comm);
M_oper = blockOper;
fillComplete();
}
//! set a component of the block operator
void BlockOperator::setBlock(UInt iblock, UInt jblock, const operator_ptr & operBlock)
{
ASSERT_PRE(M_rangeBlockMaps[iblock]->PointSameAs(operBlock->OperatorRangeMap()), "Wrong range map");
ASSERT_PRE(M_domainBlockMaps[jblock]->PointSameAs(operBlock->OperatorDomainMap()), "Wrong domain map");
M_oper(iblock, jblock) = operBlock;
}
void BlockOperator::fillComplete()
{
//Check for Structure
bool thereAreUpperDiagonalBlocks(false);
for (UInt iblock = 0; iblock < M_nBlockRows; ++iblock)
for (UInt jblock = iblock+1; jblock < M_nBlockCols; ++jblock)
if (M_oper(iblock,jblock).get() != 0 && M_oper(iblock, jblock)->HasNormInf() && M_oper(iblock, jblock)->NormInf()!=0 )
thereAreUpperDiagonalBlocks = true;
bool thereAreLowerDiagonalBlocks(false);
for (UInt iblock = 0; iblock < M_nBlockRows; ++iblock)
for (UInt jblock = 0; jblock < iblock; ++jblock)
if (M_oper(iblock,jblock).get() != 0 && M_oper(iblock, jblock)->HasNormInf() && M_oper(iblock, jblock)->NormInf()!=0 )
thereAreLowerDiagonalBlocks = true;
for (UInt iblock = 0; iblock < M_nBlockRows; ++iblock)
for (UInt jblock = 0; jblock < M_nBlockCols; ++jblock)
{
if (M_oper(iblock,jblock).get() == 0)
{
NullOperator * nullOp(new NullOperator);
nullOp->setUp(M_domainBlockMaps[jblock], M_rangeBlockMaps[iblock]);
M_oper(iblock,jblock).reset(nullOp);
}
ASSERT(M_rangeBlockMaps[iblock]->PointSameAs(M_oper(iblock,jblock)->OperatorRangeMap()), "Wrong range map");
ASSERT(M_domainBlockMaps[jblock]->PointSameAs(M_oper(iblock,jblock)->OperatorDomainMap()), "Wrong domain map");
}
if (M_nBlockRows != M_nBlockCols)
{
M_structure = Rectangular;
return;
}
if (!thereAreLowerDiagonalBlocks && !thereAreUpperDiagonalBlocks)
{
M_structure = Diagonal;
return;
}
if (thereAreLowerDiagonalBlocks && !thereAreUpperDiagonalBlocks)
{
M_structure = LowerTriangular;
return;
}
if (!thereAreLowerDiagonalBlocks && thereAreUpperDiagonalBlocks)
{
M_structure = UpperTriangular;
return;
}
if (thereAreLowerDiagonalBlocks && thereAreUpperDiagonalBlocks)
{
M_structure = NoStructure;
return;
}
}
int BlockOperator::Apply(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{
ASSERT_PRE(X.Map().SameAs(*M_domainMap),"The map of X is not conforming with domain map.");
ASSERT_PRE(Y.Map().SameAs(*M_rangeMap), "The map of Y is not conforming with range map.");
ASSERT_PRE(X.NumVectors() == Y.NumVectors(), "The number of vectors in X and Y is different" );
int nMultiVectors( X.NumVectors() );
std::vector< boost::shared_ptr<Epetra_MultiVector> > Xblock(M_nBlockCols);
std::vector< boost::shared_ptr<Epetra_MultiVector> > Yblock(M_nBlockRows);
std::vector< boost::shared_ptr<Epetra_MultiVector> > tmpYblock(M_nBlockRows);
// Split the vector
for (UInt jblock=0; jblock < M_nBlockCols; ++jblock)
{
Xblock[jblock].reset( new Epetra_MultiVector(*M_domainBlockMapsShift[jblock], nMultiVectors) );
Xblock[jblock]->PutScalar(0.0);
Xblock[jblock]->Import(X, *M_mono2blockDomain[jblock], Insert) ;
Xblock[jblock]->ReplaceMap(*M_domainBlockMaps[jblock]) ;
}
// Allocate Space for the Solution
for (UInt iblock=0; iblock < M_nBlockRows; ++iblock)
{
Yblock[iblock].reset( new Epetra_MultiVector(*M_rangeBlockMaps[iblock], nMultiVectors) );
Yblock[iblock]->PutScalar(0.0);
tmpYblock[iblock].reset( new Epetra_MultiVector(*M_rangeBlockMaps[iblock], nMultiVectors) );
}
//Perform the mat-vec multiplications
for (UInt iblock=0; iblock < M_nBlockRows; ++iblock)
for (UInt jblock=0; jblock < M_nBlockCols; ++jblock)
{
M_oper(iblock, jblock)->Apply(*Xblock[jblock], *tmpYblock[iblock]);
Yblock[iblock]->Update(1.0, *tmpYblock[iblock], 1.0);
}
//Reassemble the results
Y.PutScalar(0.0);
for (UInt iblock=0; iblock<M_nBlockRows; ++iblock)
{
Yblock[iblock]->ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
Y.Import(*Yblock[iblock], *M_block2monoRange[iblock], Insert) ;
}
return 0;
}
int BlockOperator::ApplyInverse(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{
switch (M_structure)
{
case Diagonal:
return blockJacobi(X,Y);
break;
case LowerTriangular:
return blockLowerTriangularSolve(X,Y);
break;
case UpperTriangular:
return blockUpperTriangularSolve(X,Y);
break;
case NoStructure:
Y.Scale(std::numeric_limits<Real>::quiet_NaN( ));
return -1;
break;
case Rectangular:
Y.Scale(std::numeric_limits<Real>::quiet_NaN( ));
return -1;
break;
default:
Y.Scale(std::numeric_limits<Real>::quiet_NaN( ));
return -1;
break;
}
}
//Merge two vectors using the domain map
int BlockOperator::merge(const Epetra_MultiVector & vBlock, Epetra_MultiVector & vMono, UInt jblock) const
{
ASSERT_PRE(vBlock.Map().SameAs(*M_domainBlockMaps[jblock]), "vBlock does not have the correct map" );
ASSERT_PRE(vMono.Map().SameAs(*M_domainMap), "vMono does not have the correct map");
ASSERT_PRE(vBlock.NumVectors() == vMono.NumVectors(), "The number of vectors in vBlock and vMono is different" );
Epetra_MultiVector tmp(vBlock);
tmp.ReplaceMap(*M_domainBlockMapsShift[jblock]) ;
return vMono.Import(tmp,*M_block2monoDomain[jblock], Insert);
}
//Extract vectors using the range map
int BlockOperator::extract(Epetra_MultiVector & vBlock, const Epetra_MultiVector & vMono, UInt iblock) const
{
ASSERT_PRE(M_rangeBlockMaps[iblock]->SameAs(vBlock.Map()), "vBlock does not have the correct map");
ASSERT_PRE(M_rangeMap->SameAs(vMono.Map()), "vMono does not have the correct map");
ASSERT_PRE(vBlock.NumVectors() == vMono.NumVectors(), "The number of vectors in vBlock and vMono is different" );
vBlock.ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
vBlock.Import(vMono, *M_mono2blockDomain[iblock], Insert) ;
vBlock.ReplaceMap(*M_rangeBlockMaps[iblock]) ;
return 0;
}
int BlockOperator::split(const Epetra_MultiVector & up,
vector_container & vi) const
{
UInt block(0);
for (vector_container::iterator it=vi.begin(); it != vi.end(); ++it, ++block)
extract(**it, up, block) ;
return 0;
}
int BlockOperator::merge( Epetra_MultiVector & up, const vector_container & vi) const
{
UInt block(0);
for (vector_container::const_iterator it=vi.begin(); it != vi.end(); ++it, ++block)
merge(**it, up, block) ;
return 0;
}
//===========================================================================//
//===========================================================================//
// Private Methods //
//===========================================================================//
//===========================================================================//
// Private Functions
void BlockOperator::buildImporter(UInt nblock,
std::vector<boost::shared_ptr<Epetra_Map> > & blockMaps,
boost::shared_ptr<Epetra_Map> & fullMap,
std::vector< boost::shared_ptr<Epetra_Import> > & block2mono,
std::vector< boost::shared_ptr<Epetra_Import> > & mono2block)
{
std::vector<int> numLocIdBlock(nblock), numGlobalElBlock(nblock);