LinearEpetraOperatorBlock.cpp 20.4 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
#include <life/lifealg/LinearEpetraOperatorBlock.hpp>
uvilla's avatar
uvilla committed
41
42
43
44
45
46

namespace LifeV
{
namespace Operators
{
BlockOperator::BlockOperator():
Tiziano Passerini's avatar
Tiziano Passerini committed
47
48
49
50
        M_name("BlockOperator"),
        M_useTranspose(false),
        M_structure(NoStructure)
{};
uvilla's avatar
uvilla committed
51
52

void BlockOperator::setUp(UInt nBlocks,
Tiziano Passerini's avatar
Tiziano Passerini committed
53
54
                          const std::vector<boost::shared_ptr<Epetra_Map> > domainMap,
                          const boost::shared_ptr<Epetra_Comm> & comm)
uvilla's avatar
uvilla committed
55
{
Tiziano Passerini's avatar
Tiziano Passerini committed
56
    M_comm = comm;
uvilla's avatar
uvilla committed
57

Tiziano Passerini's avatar
Tiziano Passerini committed
58
59
60
    M_nBlockRows = nBlocks;
    M_nBlockCols = nBlocks;
    M_domainBlockMaps = domainMap;
uvilla's avatar
uvilla committed
61

Tiziano Passerini's avatar
Tiziano Passerini committed
62
63
64
    M_domainBlockMapsShift.resize(nBlocks);
    for (UInt iblock=0; iblock < nBlocks; ++iblock)
        M_domainBlockMapsShift[iblock].reset(new Epetra_Map(*M_domainBlockMaps[iblock]));
uvilla's avatar
uvilla committed
65

Tiziano Passerini's avatar
Tiziano Passerini committed
66
67
    buildImporter(nBlocks, M_domainBlockMapsShift, M_domainMap,
                  M_block2monoDomain, M_mono2blockDomain);
uvilla's avatar
uvilla committed
68

Tiziano Passerini's avatar
Tiziano Passerini committed
69
70
71
72
73
74
    //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;
uvilla's avatar
uvilla committed
75

Tiziano Passerini's avatar
Tiziano Passerini committed
76
    M_oper.resize(nBlocks, nBlocks);
uvilla's avatar
uvilla committed
77
78
79
80
81
}


//! SetUp for a "rectangular operator"
void BlockOperator::setUp(UInt nRowBlocks, UInt nColBlocks,
Tiziano Passerini's avatar
Tiziano Passerini committed
82
83
84
                          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)
uvilla's avatar
uvilla committed
85
{
Tiziano Passerini's avatar
Tiziano Passerini committed
86
    M_comm = comm;
uvilla's avatar
uvilla committed
87

Tiziano Passerini's avatar
Tiziano Passerini committed
88
89
    M_nBlockRows = nRowBlocks;
    M_nBlockCols = nColBlocks;
uvilla's avatar
uvilla committed
90

Tiziano Passerini's avatar
Tiziano Passerini committed
91
92
    M_rangeBlockMaps = rangeMap;
    M_domainBlockMaps = domainMap;
uvilla's avatar
uvilla committed
93

Tiziano Passerini's avatar
Tiziano Passerini committed
94
95
96
    M_rangeBlockMapsShift.resize(M_nBlockRows);
    for (UInt iblock=0; iblock < M_nBlockRows; ++iblock)
        M_rangeBlockMapsShift[iblock].reset(new Epetra_Map(*M_rangeBlockMaps[iblock]));
uvilla's avatar
uvilla committed
97

Tiziano Passerini's avatar
Tiziano Passerini committed
98
99
100
    M_domainBlockMapsShift.resize(M_nBlockCols);
    for (UInt jblock=0; jblock < M_nBlockCols; ++jblock)
        M_domainBlockMapsShift[jblock].reset(new Epetra_Map(*M_domainBlockMaps[jblock]));
uvilla's avatar
uvilla committed
101

Tiziano Passerini's avatar
Tiziano Passerini committed
102
103
    buildImporter(M_nBlockRows, M_rangeBlockMapsShift, M_rangeMap,
                  M_block2monoRange, M_mono2blockRange);
uvilla's avatar
uvilla committed
104

Tiziano Passerini's avatar
Tiziano Passerini committed
105
106
    buildImporter(M_nBlockCols, M_domainBlockMapsShift, M_domainMap,
                  M_block2monoDomain, M_mono2blockDomain);
uvilla's avatar
uvilla committed
107

Tiziano Passerini's avatar
Tiziano Passerini committed
108
    M_oper.resize(M_nBlockRows, M_nBlockCols);
uvilla's avatar
uvilla committed
109
110
111
112
113
114

}

//! SetUp when the operator is given like a boost::matrix
void BlockOperator::setUp(const BlockOper & blockOper, const boost::shared_ptr<Epetra_Comm> & comm)
{
Tiziano Passerini's avatar
Tiziano Passerini committed
115
116
    M_nBlockRows = blockOper.size1();
    M_nBlockCols = blockOper.size2();
uvilla's avatar
uvilla committed
117

Tiziano Passerini's avatar
Tiziano Passerini committed
118
119
    M_rangeBlockMaps.resize(M_nBlockRows);
    M_domainBlockMaps.resize(M_nBlockCols);
uvilla's avatar
uvilla committed
120

Tiziano Passerini's avatar
Tiziano Passerini committed
121
122
    for (UInt iblock=0; iblock < M_nBlockRows; ++iblock)
        M_rangeBlockMaps[iblock].reset(new Epetra_Map(blockOper(iblock,0)->OperatorRangeMap() ));
uvilla's avatar
uvilla committed
123
124


Tiziano Passerini's avatar
Tiziano Passerini committed
125
126
    for (UInt jblock=0; jblock < M_nBlockCols; ++jblock)
        M_domainBlockMaps[jblock].reset(new Epetra_Map(blockOper(0,jblock)->OperatorDomainMap() ));
uvilla's avatar
uvilla committed
127

Tiziano Passerini's avatar
Tiziano Passerini committed
128
129
130
    setUp(M_nBlockRows, M_nBlockCols, M_domainBlockMaps, M_rangeBlockMaps, comm);
    M_oper = blockOper;
    fillComplete();
uvilla's avatar
uvilla committed
131
132
133
134
135
136

}

//! set a component of the block operator
void BlockOperator::setBlock(UInt iblock, UInt jblock, const operator_ptr & operBlock)
{
Tiziano Passerini's avatar
Tiziano Passerini committed
137
138
    ASSERT_PRE(M_rangeBlockMaps[iblock]->PointSameAs(operBlock->OperatorRangeMap()), "Wrong range map");
    ASSERT_PRE(M_domainBlockMaps[jblock]->PointSameAs(operBlock->OperatorDomainMap()), "Wrong domain map");
uvilla's avatar
uvilla committed
139

Tiziano Passerini's avatar
Tiziano Passerini committed
140
    M_oper(iblock, jblock) = operBlock;
uvilla's avatar
uvilla committed
141
142
143
144
145
146
}


void BlockOperator::fillComplete()
{
//Check for Structure
Tiziano Passerini's avatar
Tiziano Passerini committed
147
148
149
150
151
152
153
154
155
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
185
186
187
188
189
190
191
192
193
194
195
196
197
    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;
    }
uvilla's avatar
uvilla committed
198
199
200
201
202
203
204
205


}


int BlockOperator::Apply(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{

Tiziano Passerini's avatar
Tiziano Passerini committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
    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);
Radu Popescu's avatar
Radu Popescu committed
220
221
         Xblock[jblock]->Import(X, *M_mono2blockDomain[jblock], Insert) ;
        Xblock[jblock]->ReplaceMap(*M_domainBlockMaps[jblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
    }

    // 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)
        {
Radu Popescu's avatar
Radu Popescu committed
236
237
            M_oper(iblock, jblock)->Apply(*Xblock[jblock], *tmpYblock[iblock]);
            Yblock[iblock]->Update(1.0, *tmpYblock[iblock], 1.0);
Tiziano Passerini's avatar
Tiziano Passerini committed
238
239
240
241
242
243
        }

    //Reassemble the results
    Y.PutScalar(0.0);
    for (UInt iblock=0; iblock<M_nBlockRows; ++iblock)
    {
Radu Popescu's avatar
Radu Popescu committed
244
245
        Yblock[iblock]->ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
        Y.Import(*Yblock[iblock], *M_block2monoRange[iblock], Insert) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
246
247
    }
    return 0;
uvilla's avatar
uvilla committed
248
249
250
251
252

}

int BlockOperator::ApplyInverse(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{
Tiziano Passerini's avatar
Tiziano Passerini committed
253
254
255
256
257
258
259
260
261
262
263
264
    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:
265
        Y.Scale(std::numeric_limits<Real>::quiet_NaN( ));
Tiziano Passerini's avatar
Tiziano Passerini committed
266
267
268
        return -1;
        break;
    case Rectangular:
269
        Y.Scale(std::numeric_limits<Real>::quiet_NaN( ));
Tiziano Passerini's avatar
Tiziano Passerini committed
270
271
272
        return -1;
        break;
    default:
273
        Y.Scale(std::numeric_limits<Real>::quiet_NaN( ));
Tiziano Passerini's avatar
Tiziano Passerini committed
274
275
276
        return -1;
        break;
    }
uvilla's avatar
uvilla committed
277
278
279
280
281
282
}

//Merge two vectors using the domain map
int BlockOperator::merge(const Epetra_MultiVector & vBlock, Epetra_MultiVector & vMono, UInt jblock) const
{

Tiziano Passerini's avatar
Tiziano Passerini committed
283
284
285
    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" );
uvilla's avatar
uvilla committed
286

Tiziano Passerini's avatar
Tiziano Passerini committed
287
    Epetra_MultiVector tmp(vBlock);
Radu Popescu's avatar
Radu Popescu committed
288
     tmp.ReplaceMap(*M_domainBlockMapsShift[jblock]) ;
uvilla's avatar
uvilla committed
289

Tiziano Passerini's avatar
Tiziano Passerini committed
290
    return vMono.Import(tmp,*M_block2monoDomain[jblock], Insert);
uvilla's avatar
uvilla committed
291
292
293
294
295
296

}

//Extract vectors using the range map
int BlockOperator::extract(Epetra_MultiVector & vBlock, const Epetra_MultiVector & vMono, UInt iblock) const
{
Tiziano Passerini's avatar
Tiziano Passerini committed
297
298
299
300
    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" );

Radu Popescu's avatar
Radu Popescu committed
301
302
303
     vBlock.ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
     vBlock.Import(vMono, *M_mono2blockDomain[iblock], Insert) ;
     vBlock.ReplaceMap(*M_rangeBlockMaps[iblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
304
    return 0;
uvilla's avatar
uvilla committed
305
306
307
}

int BlockOperator::split(const Epetra_MultiVector & up,
Tiziano Passerini's avatar
Tiziano Passerini committed
308
                         vector_container & vi) const
uvilla's avatar
uvilla committed
309
{
Tiziano Passerini's avatar
Tiziano Passerini committed
310
311
    UInt block(0);
    for (vector_container::iterator it=vi.begin(); it != vi.end(); ++it, ++block)
Radu Popescu's avatar
Radu Popescu committed
312
         extract(**it, up, block) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
313
    return 0;
uvilla's avatar
uvilla committed
314
315
316
317
}

int BlockOperator::merge( Epetra_MultiVector & up, const vector_container & vi) const
{
Tiziano Passerini's avatar
Tiziano Passerini committed
318
319
    UInt block(0);
    for (vector_container::const_iterator it=vi.begin(); it != vi.end(); ++it, ++block)
Radu Popescu's avatar
Radu Popescu committed
320
         merge(**it, up, block) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
321
    return 0;
uvilla's avatar
uvilla committed
322
323
324
325
326
327
328
329
330
331

}

//===========================================================================//
//===========================================================================//
//  Private Methods                                                          //
//===========================================================================//
//===========================================================================//
// Private Functions
void BlockOperator::buildImporter(UInt nblock,
Tiziano Passerini's avatar
Tiziano Passerini committed
332
333
334
335
                                  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
336
337
{

Tiziano Passerini's avatar
Tiziano Passerini committed
338
339
    std::vector<int> numLocIdBlock(nblock), numGlobalElBlock(nblock);
    int numLocIdMono(0), numGlobalElMono(0);
uvilla's avatar
uvilla committed
340

Tiziano Passerini's avatar
Tiziano Passerini committed
341
342
343
344
345
    for (UInt iblock=0; iblock < nblock; ++iblock)
    {
        numLocIdBlock[iblock] = blockMaps[iblock]->NumMyElements();
        numGlobalElBlock[iblock] = blockMaps[iblock]->NumGlobalElements();
    }
uvilla's avatar
uvilla committed
346

Tiziano Passerini's avatar
Tiziano Passerini committed
347
348
    numLocIdMono = std::accumulate(numLocIdBlock.begin(), numLocIdBlock.end(), numLocIdMono);
    numGlobalElMono = std::accumulate(numGlobalElBlock.begin(), numGlobalElBlock.end(), numGlobalElMono);
uvilla's avatar
uvilla committed
349

Tiziano Passerini's avatar
Tiziano Passerini committed
350
351
352
353
354
    std::vector<int> shift(nblock+1);
    std::vector<int>::iterator itShift(shift.begin());
    *itShift = 0;
    itShift++;
    std::partial_sum(numGlobalElBlock.begin(), numGlobalElBlock.end(), itShift);
uvilla's avatar
uvilla committed
355

Tiziano Passerini's avatar
Tiziano Passerini committed
356
357
    std::vector<std::vector<int> > myGlobIdBlock(nblock);
    std::vector<int> myGlobIdMono(numGlobalElMono);
uvilla's avatar
uvilla committed
358

Tiziano Passerini's avatar
Tiziano Passerini committed
359
360
    std::vector<int>::iterator itMono;
    itMono = myGlobIdMono.begin();
uvilla's avatar
uvilla committed
361

Tiziano Passerini's avatar
Tiziano Passerini committed
362
363
364
365
366
367
368
    for (std::vector<std::vector<int> >::iterator myGlobIdBlockIt = myGlobIdBlock.begin();
            myGlobIdBlockIt != myGlobIdBlock.end();
            ++myGlobIdBlockIt)
    {
        UInt iblock = (UInt) (myGlobIdBlockIt - myGlobIdBlock.begin());
        myGlobIdBlockIt->resize(numLocIdBlock[iblock]);
        blockMaps[iblock]->MyGlobalElements( &(*myGlobIdBlockIt)[0] );
uvilla's avatar
uvilla committed
369

Tiziano Passerini's avatar
Tiziano Passerini committed
370
371
        for (std::vector<int>::iterator it=myGlobIdBlockIt->begin(); it!=myGlobIdBlockIt->end(); ++it)
            *it += shift[iblock];
uvilla's avatar
uvilla committed
372

Tiziano Passerini's avatar
Tiziano Passerini committed
373
        itMono = std::copy(myGlobIdBlockIt->begin(), myGlobIdBlockIt->end(), itMono);
uvilla's avatar
uvilla committed
374

Tiziano Passerini's avatar
Tiziano Passerini committed
375
376
        blockMaps[iblock].reset(new Epetra_Map(numGlobalElBlock[iblock], myGlobIdBlockIt->size(),
                                               &(*myGlobIdBlockIt)[0], shift[iblock]+1, *M_comm));
uvilla's avatar
uvilla committed
377

Tiziano Passerini's avatar
Tiziano Passerini committed
378
    }
uvilla's avatar
uvilla committed
379

Tiziano Passerini's avatar
Tiziano Passerini committed
380
    fullMap.reset(new Epetra_Map(numGlobalElMono, myGlobIdMono.size(), &myGlobIdMono[0], 1, *M_comm));
uvilla's avatar
uvilla committed
381

Tiziano Passerini's avatar
Tiziano Passerini committed
382
383
    block2mono.resize(nblock);
    mono2block.resize(nblock);
uvilla's avatar
uvilla committed
384

Tiziano Passerini's avatar
Tiziano Passerini committed
385
386
387
388
389
    for (UInt iblock=0; iblock<nblock; ++iblock)
    {
        block2mono[iblock].reset(new Epetra_Import(*fullMap, *blockMaps[iblock]));
        mono2block[iblock].reset(new Epetra_Import(*blockMaps[iblock], *fullMap));
    }
uvilla's avatar
uvilla committed
390
391
392
393
394
395
396

}


int BlockOperator::blockJacobi(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{

Tiziano Passerini's avatar
Tiziano Passerini committed
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
    ASSERT_PRE(M_nBlockRows == M_nBlockCols, "The operator must be squared");
    ASSERT_PRE(M_domainMap->SameAs(*M_rangeMap), "The operator must be squared");
    ASSERT_PRE(Y.Map().SameAs(*M_domainMap),"The map of Y is not conforming with domain map.");
    ASSERT_PRE(X.Map().SameAs(*M_rangeMap), "The map of X 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);

    // 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);
Radu Popescu's avatar
Radu Popescu committed
413
414
        Xblock[jblock]->Import(X, *M_mono2blockDomain[jblock], Insert) ;
        Xblock[jblock]->ReplaceMap(*M_domainBlockMaps[jblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
415
416
417
418
419
420
421
422
423
424
425
    }

    // 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);
    }

    //Perform the mat-vec multiplications
    for (UInt iblock = 0; iblock < M_nBlockRows; ++iblock)
Radu Popescu's avatar
Radu Popescu committed
426
        M_oper(iblock, iblock)->ApplyInverse(*Xblock[iblock], *Yblock[iblock]);
Tiziano Passerini's avatar
Tiziano Passerini committed
427
428
429
430
431

    //Reassemble the results
    Y.PutScalar(0.0);
    for (UInt iblock=0; iblock<M_nBlockRows; ++iblock)
    {
Radu Popescu's avatar
Radu Popescu committed
432
433
        Yblock[iblock]->ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
         Y.Import(*Yblock[iblock], *M_block2monoRange[iblock], Insert) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
434
435
436
    }

    return 0;
uvilla's avatar
uvilla committed
437
438
439
440
441
}

int BlockOperator::blockLowerTriangularSolve(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{

Tiziano Passerini's avatar
Tiziano Passerini committed
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
    ASSERT_PRE(M_nBlockRows == M_nBlockCols, "The operator must be squared");
    ASSERT_PRE(M_domainMap->SameAs(*M_rangeMap), "The operator must be squared");
    ASSERT_PRE(Y.Map().SameAs(*M_domainMap),"The map of Y is not conforming with domain map.");
    ASSERT_PRE(X.Map().SameAs(*M_rangeMap), "The map of X 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> > Zblock(M_nBlockCols);

    // 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);
Radu Popescu's avatar
Radu Popescu committed
459
460
        Xblock[jblock]->Import(X, *M_mono2blockDomain[jblock], Insert) ;
        Xblock[jblock]->ReplaceMap(*M_domainBlockMaps[jblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
461
462
463
464
465
466
467
468
469
470
471
472
473
474

        Zblock[jblock].reset( new Epetra_MultiVector(*M_domainBlockMaps[jblock], nMultiVectors) );
    }

    // 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);
    }

    //Perform the mat-vec multiplications
    for (UInt iblock = 0; iblock < M_nBlockRows; ++iblock)
    {
Radu Popescu's avatar
Radu Popescu committed
475
        M_oper(iblock, iblock)->ApplyInverse(*Xblock[iblock], *Yblock[iblock]);
Tiziano Passerini's avatar
Tiziano Passerini committed
476
477
        for (UInt kblock = iblock; kblock < M_nBlockRows; ++kblock)
        {
Radu Popescu's avatar
Radu Popescu committed
478
             M_oper(kblock, iblock)->Apply(*Yblock[iblock], *Zblock[kblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
479
480
481
482
483
484
485
486
            Xblock[kblock]->Update(-1.0, *Zblock[kblock], 1.0);
        }
    }

    //Reassemble the results
    Y.PutScalar(0.0);
    for (UInt iblock=0; iblock<M_nBlockRows; ++iblock)
    {
Radu Popescu's avatar
Radu Popescu committed
487
488
        Yblock[iblock]->ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
         Y.Import(*Yblock[iblock], *M_block2monoRange[iblock], Insert) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
489
490
491
    }

    return 0;
uvilla's avatar
uvilla committed
492
493
494
495
496
}

int BlockOperator::blockUpperTriangularSolve(const Epetra_MultiVector & X, Epetra_MultiVector & Y) const
{

Tiziano Passerini's avatar
Tiziano Passerini committed
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
    ASSERT_PRE(M_nBlockRows == M_nBlockCols, "The operator must be squared");
    ASSERT_PRE(M_domainMap->SameAs(*M_rangeMap), "The operator must be squared");
    ASSERT_PRE(Y.Map().SameAs(*M_domainMap),"The map of Y is not conforming with domain map.");
    ASSERT_PRE(X.Map().SameAs(*M_rangeMap), "The map of X 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> > Zblock(M_nBlockCols);

    // 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);
Radu Popescu's avatar
Radu Popescu committed
514
515
        Xblock[jblock]->Import(X, *M_mono2blockDomain[jblock], Insert) ;
        Xblock[jblock]->ReplaceMap(*M_domainBlockMaps[jblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
516
517
518
519
520
521
522
523
524
525
526
527
528
529

        Zblock[jblock].reset( new Epetra_MultiVector(*M_domainBlockMaps[jblock], nMultiVectors) );
    }

    // 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);
    }

    //Perform the mat-vec multiplications
    for (int iblock = M_nBlockRows - 1 ; iblock > -1 ; --iblock)
    {
Radu Popescu's avatar
Radu Popescu committed
530
        M_oper(iblock, iblock)->ApplyInverse(*Xblock[iblock], *Yblock[iblock]);
Tiziano Passerini's avatar
Tiziano Passerini committed
531
532
        for (int kblock = 0; kblock < iblock; ++kblock)
        {
Radu Popescu's avatar
Radu Popescu committed
533
             M_oper(kblock, iblock)->Apply(*Yblock[iblock], *Zblock[kblock]) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
534
535
536
537
538
539
540
541
            Xblock[kblock]->Update(-1.0, *Zblock[kblock], 1.0);
        }
    }

    //Reassemble the results
    Y.PutScalar(0.0);
    for (UInt iblock=0; iblock<M_nBlockRows; ++iblock)
    {
Radu Popescu's avatar
Radu Popescu committed
542
543
        Yblock[iblock]->ReplaceMap(*M_rangeBlockMapsShift[iblock]) ;
        Y.Import(*Yblock[iblock], *M_block2monoRange[iblock], Insert) ;
Tiziano Passerini's avatar
Tiziano Passerini committed
544
545
546
    }

    return 0;
uvilla's avatar
uvilla committed
547
548
549
550
551
}


} /* end namespace Operators */
} /*end namespace */