LinearEpetraOperatorBlock.hpp 9.83 KB
Newer Older
Tiziano Passerini's avatar
Tiziano Passerini committed
1
//@HEADER
uvilla's avatar
uvilla committed
2
/*
Tiziano Passerini's avatar
Tiziano Passerini committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
*******************************************************************************

    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.
uvilla's avatar
uvilla committed
38
39
 */

40
41
#ifndef LINEAREPETRAOPERATORBLOCK_H
#define LINEAREPETRAOPERATORBLOCK_H 1
uvilla's avatar
uvilla committed
42

43
#include <life/lifealg/LinearEpetraOperator.hpp>
44
#include <boost/numeric/ublas/matrix.hpp>
perego's avatar
perego committed
45
46
47
48
49
// 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"

uvilla's avatar
uvilla committed
50
51
#include <Epetra_Import.h>

uvilla's avatar
uvilla committed
52
53
54
55
// 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"
uvilla's avatar
uvilla committed
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
namespace LifeV
{
namespace Operators
{
//! @class BlockOperator
/*! @brief A abstract class for handling n-by-m block operators
 * This class inherits from LifeV::LinearOperator.
 *
 * The Transpose is not supported yet.
 */

class BlockOperator: public LinearOperator
{
public:

Tiziano Passerini's avatar
Tiziano Passerini committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    //! @name Public Typedefs
    //@{
    typedef LinearOperator super;
    typedef Epetra_Operator raw_operator;
    typedef boost::shared_ptr<raw_operator> operator_ptr;
    typedef boost::numeric::ublas::matrix<operator_ptr> BlockOper;

    typedef Epetra_MultiVector raw_vector;
    typedef boost::shared_ptr<raw_vector> vector_ptr;
    typedef std::vector<vector_ptr> vector_container;

    enum Structure {Diagonal, LowerTriangular, UpperTriangular, NoStructure, Rectangular};
    //@}

    //! Empty Constructor
    BlockOperator();

    //! @name Set Methods
    //@{
    //! SetUp for a "square operator"
    /*!
     * @param nBlocks:   number of blocks in a square n-by-n matrix
     * @param domainMap: the map of a vector in the domain of \e this
     *                   is obtained by concatenating the block maps in
     *                   domainMap.
     *                   rangeMap is assumed to be the same of domainMap.
     * @param comm:      the communicator.
     */
    void setUp(UInt nBlocks,
               const std::vector<boost::shared_ptr<Epetra_Map> > domainMap,
               const boost::shared_ptr<Epetra_Comm> & comm);

    //! SetUp for a "rectangular operator"
    /*!
     * @param nRowBlocks, nColBlocks:  \e this is a nRowBlocks-by-nColBlocks
     *                                  block operator
     * @param domainMap: the map of a vector in the domain of \e this
     *                   is obtained by concatenating the block maps in
     *                   domainMap.
     * @param rangeMap:  the map of a vector in the range of \e this
     *                   is obtained by concatenating the block maps in
     *                   rangeMap.
     * @param comm:      the communicator.
     */
    void 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
              );

    //! SetUp when the operator is given like a boost::matrix
    /*!
     * @param blockOper: a dense matrix to describe the block operator
     * @param comm:      the communicator
     */
    void setUp(const BlockOper & blockOper, const boost::shared_ptr<Epetra_Comm> & comm);

    //! set a component of the block operator
    /*!
     * @param iblock, jblock: The position of the block is (iblock, jblock).
     * @param operBlock     : an operator_ptr representing the block
     */
    void setBlock(UInt iblock, UInt jblock, const operator_ptr & operBlock);

    //! Complete the block matrix with null operators
    void fillComplete();

    //! If true the transpose of the operator will be computed.
    /*
     * Not Supported yet
     */
    int SetUseTranspose(bool useTranspose) {M_useTranspose = useTranspose; return -1;}

    //! Compute Y = Op*X;
    virtual int Apply(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const;
    //! Compute Y = Op\X;
    /*!
     * ApplyInverse is implemented for the matrices with block diagonal, lowerTriangular, upperTriangular form
     */
    virtual int ApplyInverse(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const;

    //! Compute the Inf norm of the operator
    /*!
     * Not implemented yet.
     */
uvilla's avatar
uvilla committed
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
    double NormInf() const {return -1;}

    //! Returns a character string describing the operator
    virtual const char * Label() const {return M_name.c_str();}

    //! Returns the current UseTranspose setting.
    bool UseTranspose() const {return M_useTranspose;}

    //! Returns true if the \e this object can provide an approximate Inf-norm, false otherwise.
    bool HasNormInf() const {return false;}

    //! Returns a pointer to the Epetra_Comm communicator associated with this operator.
    const Epetra_Comm & Comm() const {return *M_comm;}

    //! Returns the Epetra_Map object associated with the domain of this operator.
    const Epetra_Map & OperatorDomainMap() const {return *M_domainMap;}
    //! Returns the Epetra_Map object associated with the domain of this operator as a pointer
    const boost::shared_ptr<Epetra_Map> & OperatorDomainMap_ptr() const {return M_domainMap;}

    //! Returns the Epetra_Map object associated with the range of this operator.
    const Epetra_Map & OperatorRangeMap() const {return *M_rangeMap;}
    //! Returns the Epetra_Map object associated with the range of this operator as a pointer
    const boost::shared_ptr<Epetra_Map> & OperatorRangeMap_ptr() const {return M_rangeMap;}

    //! Merge two vectors using the domain map
    int merge(const Epetra_MultiVector & vBlock, Epetra_MultiVector & vMono, UInt jblock) const;
    //! Extract vectors using the range map
    int extract(Epetra_MultiVector & vBlock, const Epetra_MultiVector & vMono, UInt jblock) const;

Tiziano Passerini's avatar
Tiziano Passerini committed
185
186
    int split(const Epetra_MultiVector & up,
              vector_container & vi) const;
uvilla's avatar
uvilla committed
187

Tiziano Passerini's avatar
Tiziano Passerini committed
188
189
    int merge(      Epetra_MultiVector & up,
                    const vector_container & vi) const;
uvilla's avatar
uvilla committed
190
191
192
193
194
195
196
197
198
199


protected:
    //! Y = diag(block(i,i)^-1)*X
    int blockJacobi(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const;
    //! Y = backwardsubstitution(X)
    int blockUpperTriangularSolve(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const;
    //! Y = forwardsubstitution(X)
    int blockLowerTriangularSolve(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const;
    //! Change the name of the operator, (avaible for derivate classes).
Tiziano Passerini's avatar
Tiziano Passerini committed
200
    void setName(const std::string & name) {M_name =name;}
uvilla's avatar
uvilla committed
201
202
203
private:
    //! Construct the maps and the importers
    void buildImporter(UInt nblocks,
Tiziano Passerini's avatar
Tiziano Passerini committed
204
205
206
207
                       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 );
uvilla's avatar
uvilla committed
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226

    //! Number of blocks in each row
    UInt M_nBlockRows;
    //! Number of blocks in each column
    UInt M_nBlockCols;

    //! Name of the object
    std::string M_name;
    //! Communicator
    boost::shared_ptr<Epetra_Comm> M_comm;

    //! @name Maps
    //@{
    //! Domain Map
    boost::shared_ptr<Epetra_Map> M_domainMap;
    //! Range Map
    boost::shared_ptr<Epetra_Map> M_rangeMap;

    //! Domain Map of each block
Tiziano Passerini's avatar
Tiziano Passerini committed
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
    std::vector< boost::shared_ptr<Epetra_Map> > M_domainBlockMaps;
    //! Range Map of each block
    std::vector< boost::shared_ptr<Epetra_Map> > M_rangeBlockMaps;

    //! Shifted domain Map of each block
    std::vector< boost::shared_ptr<Epetra_Map> > M_domainBlockMapsShift;
    //! Shifted domain Map of each block
    std::vector< boost::shared_ptr<Epetra_Map> > M_rangeBlockMapsShift;
    //@}

    //! block operator represented like a dense matrix of pointers to Operators
    BlockOper M_oper;

    //! @name Importers
    //@{
    //! merge block domain vectors in the monolithic domain vector
    std::vector< boost::shared_ptr<Epetra_Import> > M_block2monoDomain;
    //! split the monolithic domain vector in the block domain vectors
    std::vector< boost::shared_ptr<Epetra_Import> > M_mono2blockDomain;
    //! merge block range vectors in the monolithic range vector
    std::vector< boost::shared_ptr<Epetra_Import> > M_block2monoRange;
    //! split the monolithic domain vector in the block domain vectors
    std::vector< boost::shared_ptr<Epetra_Import> > M_mono2blockRange;
    //@}

    //! whenever transpose should be used
    bool M_useTranspose;

    //! structure of the block operator (Diagonal, LowerDiagonal, UpperDiagonal, NoStructure)
    Structure M_structure;
uvilla's avatar
uvilla committed
257
258
259
260
};

} /*end namespace Operators*/
} /*end namespace */
261
#endif /* LINEAREPETRAOPERATORBLOCK_H */