Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
BlockMatrix_def.hpp
1#ifndef BLOCKMATRIX_DEF_hpp
2#define BLOCKMATRIX_DEF_hpp
3#include "BlockMatrix_decl.hpp"
12
13namespace FEDD {
14
15template <class SC, class LO, class GO, class NO>
16BlockMatrix<SC,LO,GO,NO>::BlockMatrix():
17blockMatrix_(0),
18blockMap_(),
19mergedMatrix_(),
20mergedMap_(),
21globalBlockOffsets_(),
22localBlockOffsets_()
23{
24 blockMap_ = Teuchos::rcp( new BlockMap_Type( 0 ) );
25}
26
27template <class SC, class LO, class GO, class NO>
28BlockMatrix<SC,LO,GO,NO>::BlockMatrix(UN size):
29blockMatrix_(size),
30blockMap_(),
31mergedMatrix_(),
32mergedMap_(),
33globalBlockOffsets_(),
34localBlockOffsets_()
35{
36 blockMap_ = Teuchos::rcp( new BlockMap_Type( size ) );
37}
38
39
40template <class SC, class LO, class GO, class NO>
41BlockMatrix<SC,LO,GO,NO>::BlockMatrix(BlockMatrixPtr_Type bMatrixIn):
42blockMatrix_(bMatrixIn->size()),
43blockMap_(),
44mergedMatrix_(),
45mergedMap_(),
46globalBlockOffsets_(),
47localBlockOffsets_()
48{
49 blockMap_ = Teuchos::rcp( new BlockMap_Type( bMatrixIn->size() ) );
50
51 for(UN i = 0; i < bMatrixIn->size(); i++)
52 {
53 for(UN j = 0; j < bMatrixIn->size(); j++)
54 {
55 // copy matrix and then insert
56 if(bMatrixIn->blockExists(i,j))
57 {
58 MatrixPtr_Type matrixTmp = bMatrixIn->getBlock(i,j);
59 MatrixPtr_Type matrixCopy = Teuchos::rcp(new Matrix_Type( matrixTmp ) );
60 this->addBlock(matrixCopy, i, j);
61 }
62 }
63 }
64}
65
66
67template <class SC, class LO, class GO, class NO>
68BlockMatrix<SC,LO,GO,NO>::~BlockMatrix()
69{
70
71}
72
73template <class SC, class LO, class GO, class NO>
74int BlockMatrix<SC,LO,GO,NO>::size() const{
75 return blockMatrix_.size();
76}
77
78template <class SC, class LO, class GO, class NO>
79void BlockMatrix<SC,LO,GO,NO>::resize(UN size) {
80 blockMatrix_.resize( size );
81 blockMap_->resize( size );
82}
83
84template <class SC, class LO, class GO, class NO>
85typename BlockMatrix<SC,LO,GO,NO>::MatrixPtr_Type BlockMatrix<SC,LO,GO,NO>::getBlock(int i, int j){
86
87 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_[i][j].is_null(), std::runtime_error,"Block in BlockMatrix which you tried to access is null.");
88 return blockMatrix_[i][j];
89}
90
91template <class SC, class LO, class GO, class NO>
92typename BlockMatrix<SC,LO,GO,NO>::MatrixConstPtr_Type BlockMatrix<SC,LO,GO,NO>::getBlockConst(int i, int j) const{
93
94 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_[i][j].is_null(), std::runtime_error,"Block in BlockMatrix which you tried to access is null.");
95 return blockMatrix_[i][j];
96}
97
98template <class SC, class LO, class GO, class NO>
99bool BlockMatrix<SC,LO,GO,NO>::blockExists(int i, int j) const{
100 return ( !blockMatrix_[i][j].is_null() );
101}
102
103template <class SC, class LO, class GO, class NO>
104void BlockMatrix<SC,LO,GO,NO>::addBlock(const MatrixPtr_Type& matrix, int i, int j){
105 UN size = blockMatrix_.size();
106 if (i>size-1 || j>size-1)
107 blockMatrix_.resize( std::max(i,j)+1 );
108
109 if ( blockExists(i,j) )
110 blockMatrix_[i][j].reset();
111
112 blockMatrix_[i][j] = matrix;
113
114
115 blockMap_->addBlock( matrix->getMap(), i);
116}
117
118template <class SC, class LO, class GO, class NO>
119void BlockMatrix<SC,LO,GO,NO>::merge(){
120 if ( mergedMap_.is_null() ) {
121 blockMap_->merge();
122 mergedMap_ = Teuchos::rcp_const_cast<Map_Type>(blockMap_->getMergedMap());
123 }
124 LO maxNumEntries = 0;
125
126 for (UN i=0; i<blockMatrix_.size(); i++) {
127 LO maxNumEntriesBlockRow = 0;
128 for (UN j=0; j<blockMatrix_.size(); j++) {
129 if ( !blockMatrix_[i][j].is_null() )
130 maxNumEntriesBlockRow += blockMatrix_[i][j]->getGlobalMaxNumRowEntries();
131 }
132 if (maxNumEntriesBlockRow > maxNumEntries)
133 maxNumEntries = maxNumEntriesBlockRow;
134 }
135
136// if ( mergedMatrix_.is_null() ){
137 mergedMatrix_ = Teuchos::rcp( new Matrix_Type( mergedMap_, maxNumEntries ) ); // we should identify the number of entries per row.
138 this->determineLocalOffsets();
139 this->determineGlobalOffsets();
140
141 for (UN i=0; i<blockMatrix_.size(); i++) {
142 for (UN j=0; j<blockMatrix_.size(); j++) {
143 if ( !blockMatrix_[i][j].is_null() ){
144 mergeBlockNew( i, j );
145 }
146 }
147 }
148 mergedMatrix_->fillComplete( mergedMap_, mergedMap_ ); // only for square matrices!
149// }
150// else {
151// mergedMatrix_->resumeFill();
152// for (UN i=0; i<blockMatrix_.size(); i++) {
153// for (UN j=0; j<blockMatrix_.size(); j++) {
154// if ( !blockMatrix_[i][j].is_null() )
155// TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"We need to implement mergeBlockAgain");
156// }
157// }
158// }
160//
161// mergedMatrix_->fillComplete( mergedMap_, mergedMap_ ); // only for square matrices!
162
163}
164
165template <class SC, class LO, class GO, class NO>
166void BlockMatrix<SC,LO,GO,NO>::determineLocalOffsets(){
167
168
169 typedef Teuchos::ScalarTraits<LO> LOST;
170
171 Teuchos::Array<LO> offsetRows( blockMatrix_.size() - 1, LOST::zero() );
172 Teuchos::Array<LO> offsetCols( blockMatrix_.size() - 1, LOST::zero() );
173 for (UN row=0; row<blockMatrix_.size() - 1; row++) {
174 bool foundOffset = false;
175 for (UN col=0; col<blockMatrix_.size() && !foundOffset; col++) {
176 if ( !blockMatrix_[row][col].is_null() ){
177 LO mli = blockMatrix_[row][col]->getMap("row")->getMaxLocalIndex();
178 if (mli>=0)
179 offsetRows[row] = mli + 1;
180 }
181 }
182 }
183 for (UN col=0; col<blockMatrix_.size() - 1; col++) {
184 bool foundOffset = false;
185 for (UN row=0; row<blockMatrix_.size() && !foundOffset; row++) {
186 if ( !blockMatrix_[row][col].is_null() ){
187 LO mli = blockMatrix_[row][col]->getMap("col")->getMaxLocalIndex();
188 if (mli>=0)
189 offsetCols[col] = mli + 1;
190 }
191 }
192 }
193
194 localBlockOffsets_ =
195 Teuchos::rcp( new SmallMatrix< LOTuple_Type >( blockMatrix_.size() ) );
196
197 LO offsetRow = LOST::zero();
198 for (UN row=0; row<blockMatrix_.size(); row++) {
199 LO offsetCol = LOST::zero();
200 for (UN col=0; col<blockMatrix_.size(); col++) {
201 LOTuple_Type tuple = Teuchos::tuple( offsetRow, offsetCol );
202 (*localBlockOffsets_)[row][col] = tuple;
203 if (offsetCols.size() > col)
204 offsetCol += offsetCols[col];
205 }
206 if (offsetRows.size() > row)
207 offsetRow += offsetRows[row];
208 }
209}
210
211template <class SC, class LO, class GO, class NO>
212void BlockMatrix<SC,LO,GO,NO>::determineGlobalOffsets(){
213
214 typedef Teuchos::ScalarTraits<GO> GOST;
215
216 Teuchos::Array<GO> offsetRows( blockMatrix_.size() - 1, GOST::zero() );
217 Teuchos::Array<GO> offsetCols( blockMatrix_.size() - 1, GOST::zero() );
218 for (UN row=0; row<blockMatrix_.size() - 1; row++) {
219 bool foundOffset = false;
220 for (UN col=0; col<blockMatrix_.size() && !foundOffset; col++) {
221 if ( !blockMatrix_[row][col].is_null() )
222 {
223 offsetRows[row] = blockMatrix_[row][col]->getMap("row")->getMaxAllGlobalIndex() + 1;
224 foundOffset = true;
225 }
226 }
227 }
228 for (UN col=0; col<blockMatrix_.size() - 1; col++) {
229 bool foundOffset = false;
230 for (UN row=0; row<blockMatrix_.size() && !foundOffset; row++) {
231 if ( !blockMatrix_[row][col].is_null() )
232 {
233 offsetCols[col] = blockMatrix_[row][col]->getMap("col")->getMaxAllGlobalIndex() + 1;
234 foundOffset = true;
235 }
236 }
237 }
238
239 globalBlockOffsets_ =
240 Teuchos::rcp( new SmallMatrix< GOTuple_Type >( blockMatrix_.size() ) );
241
242 GO offsetRow = GOST::zero();
243 for (UN row=0; row<blockMatrix_.size(); row++) {
244 GO offsetCol = GOST::zero();
245 for (UN col=0; col<blockMatrix_.size(); col++) {
246 GOTuple_Type tuple = Teuchos::tuple( offsetRow, offsetCol );
247 (*globalBlockOffsets_)[row][col] = tuple;
248 if (offsetCols.size() > col)
249 offsetCol += offsetCols[col];
250 }
251 if (offsetRows.size() > row)
252 offsetRow += offsetRows[row];
253 }
254}
255
256template <class SC, class LO, class GO, class NO>
257void BlockMatrix<SC,LO,GO,NO>::mergeBlockNew(UN blockRow, UN blockCol){
258
259 MatrixPtr_Type matrix = blockMatrix_[blockRow][blockCol];
260 if ( matrix->isLocallyIndexed() ) {
261
262 Teuchos::ArrayView<const SC> values;
263 Teuchos::ArrayView<const LO> indices;
264 MapConstPtr_Type colMap = matrix->getMap("col");
265 MapConstPtr_Type rowMap = matrix->getMap("row");
266
267 for (UN i=0; i<matrix->getNodeNumRows(); i++) {
268 matrix->getLocalRowView( i, indices, values );
269 Teuchos::Array<GO> indicesGlobal( indices.size() );
270 for (UN j=0; j<indices.size(); j++) {
271 indicesGlobal[j] = colMap->getGlobalElement( indices[j] ) + (*globalBlockOffsets_)[blockRow][blockCol][1];
272// if (indicesGlobal[j]<0) {
273// std::cout << "blockRow:"<< blockRow << " blockCol:" << blockCol << " "<<i << " "<< j << " indicesGlobal[j]:"<< indicesGlobal[j]<<" indices.size():"<<indices.size() << " indices[j]:"<< indices[j] <<" gBO:" << (*globalBlockOffsets_)[blockRow][blockCol][1] << std::endl;
274// }
275 }
276 GO offset = (*globalBlockOffsets_)[blockRow][blockCol][0];
277 mergedMatrix_->insertGlobalValues( rowMap->getGlobalElement( i ) + offset, indicesGlobal(), values );
278 }
279
280 }
281// else if ( blockMatrix_[blockRow][blockCol]->isGloballyIndexed() ) {
282// TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"We need to implement mergeBlockNew for isGloballyIndexed()==true.");
283// }
284 else
285 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,"Call fillComplete() before you merge the matrix blocks.");
286
287}
288
289
290template <class SC, class LO, class GO, class NO>
291void BlockMatrix<SC,LO,GO,NO>::print(Teuchos::EVerbosityLevel verbLevel){
292
293 for (UN i=0; i<blockMatrix_.size(); i++) {
294 for (UN j=0; j<blockMatrix_.size(); j++) {
295 if ( !blockMatrix_[i][j].is_null() ) {
296 blockMatrix_[i][j]->print( verbLevel );
297 }
298 }
299 }
300}
301
302template <class SC, class LO, class GO, class NO>
303void BlockMatrix<SC,LO,GO,NO>::printMerge(Teuchos::EVerbosityLevel verbLevel){
304
305 if ( !mergedMatrix_.is_null() )
306 mergedMatrix_->print( verbLevel );
307
308}
309
310template <class SC, class LO, class GO, class NO>
311Teuchos::RCP<const Thyra::LinearOpBase<SC> > BlockMatrix<SC,LO,GO,NO>::getThyraLinOp(){
312
313 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() == 0, std::logic_error,"BlockMatrix size is 0.");
314 Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyraLinOp;
315 if (blockMatrix_.size() == 1){
316 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_[0][0].is_null(), std::runtime_error, "Block in BlockMatrix is null.");
317 thyraLinOp = blockMatrix_[0][0]->getThyraLinOp( );
318 }
319 else{
320 this->merge();
321 thyraLinOp = mergedMatrix_->getThyraLinOp( );
322 }
323 return thyraLinOp;
324}
325
326template <class SC, class LO, class GO, class NO>
327Teuchos::RCP< const Thyra::BlockedLinearOpBase<SC> > BlockMatrix<SC,LO,GO,NO>::getThyraLinBlockOp() const{
328
329 Teuchos::RCP< Thyra::DefaultBlockedLinearOp< SC > > dblOp = Thyra::defaultBlockedLinearOp<SC> ();
330 dblOp->beginBlockFill( this->size(), this->size() );
331
332 for (int i=0; i<this->size(); i++) {
333 for (int j=0; j<this->size(); j++) {
334 if ( this->blockExists(i,j) ) {
335 dblOp->setBlock( i, j, blockMatrix_[i][j]->getThyraLinOp( ) );
336 }
337 }
338 }
339
340 dblOp->endBlockFill();
341 Teuchos::RCP< const Thyra::BlockedLinearOpBase< SC > > blOpConst = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase< SC > > (dblOp);
342// Teuchos::RCP< const Thyra::BlockedLinearOpBase< SC > > blOpConst = Teuchos::rcp_const_cast<const Thyra::BlockedLinearOpBase< SC > > (blOp);
343 return blOpConst;
344}
345
346template <class SC, class LO, class GO, class NO>
347void BlockMatrix<SC,LO,GO,NO>::apply(const BlockMultiVector_Type& X,
348 BlockMultiVector_Type& Y ) const{
349
350 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() == 0, std::logic_error,"BlockMatrix size is 0.");
351 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() != X.size(), std::logic_error,"BlockMatrix and BlockMultiVector have different sizes.");
352 TEUCHOS_TEST_FOR_EXCEPTION( Y.size() != X.size(), std::logic_error,"BlockMultiVectors have different sizes.");
353
354 for (int i=0; i<Y.size(); i++)
355 Y[i]->putScalar( Teuchos::ScalarTraits<SC>::zero() );
356
357 typedef Teuchos::ScalarTraits<SC> ST;
358
359 for (int i=0; i<blockMatrix_.size(); i++) {
360 for (int j=0; j<blockMatrix_.size(); j++){
361 MultiVectorConstPtr_Type tmpX = X[j];
362 if (!blockMatrix_[i][j].is_null())
363 {
364 blockMatrix_[i][j]->apply( *tmpX, *Y[i], Teuchos::NO_TRANS, ST::one(), ST::one() );
365 }
366 }
367 }
368}
369
370template <class SC, class LO, class GO, class NO>
371void BlockMatrix<SC,LO,GO,NO>::apply(const BlockMultiVector_Type& X,
372 BlockMultiVector_Type& Y,
373 const SmallMatrix<SC>& coeff) const{
374
375 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() != coeff.size(), std::logic_error,"BlockMatrix size and coefficient size are different.");
376 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() == 0, std::logic_error,"BlockMatrix size is 0.");
377 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() != X.size(), std::logic_error,"BlockMatrix and BlockMultiVector have different sizes.");
378 TEUCHOS_TEST_FOR_EXCEPTION( Y.size() != X.size(), std::logic_error,"BlockMultiVectors have different sizes.");
379
380 for (int i=0; i<Y.size(); i++)
381 Y[i]->putScalar( Teuchos::ScalarTraits<SC>::zero() );
382
383 typedef Teuchos::ScalarTraits<SC> ST;
384
385 for (int i=0; i<blockMatrix_.size(); i++) {
386 for (int j=0; j<blockMatrix_.size(); j++){
387 MultiVectorConstPtr_Type tmpX = X[j];
388 if (!blockMatrix_[i][j].is_null())
389 blockMatrix_[i][j]->apply( *tmpX, *Y[i], Teuchos::NO_TRANS, coeff[i][j], ST::one() );
390 }
391 }
392}
393
394template <class SC, class LO, class GO, class NO>
395void BlockMatrix<SC,LO,GO,NO>::addMatrix(const SmallMatrix<SC>& coeffAlpha, const BlockMatrixPtr_Type &B, const SmallMatrix<SC>& coeffBeta){
396 //B = alpha*A + beta*B.
397 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() != coeffAlpha.size(), std::logic_error,"BlockMatrix size and coefficient size are different.");
398 TEUCHOS_TEST_FOR_EXCEPTION( coeffAlpha.size() != coeffBeta.size(), std::logic_error,"Coefficients have different sizes.");
399 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() == 0, std::logic_error,"BlockMatrix size is 0.");
400 TEUCHOS_TEST_FOR_EXCEPTION( blockMatrix_.size() != B->size(), std::logic_error,"BlockMatrices have different sizes.");
401
402 // CAREFUL TO NOT IGNORE ANYTHING
403 // B must have the correct blocks!
404 for (int i=0; i<blockMatrix_.size(); i++) {
405 for (int j=0; j<blockMatrix_.size(); j++){
406 if (blockMatrix_[i][j].is_null()){
407 if ( coeffBeta[i][j] != Teuchos::ScalarTraits<SC>::one() && B->blockExists(i,j) )
408 {
409 B->getBlock(i,j)->scale( coeffBeta[i][j] );
410 }
411 }
412 else
413 blockMatrix_[i][j]->addMatrix( coeffAlpha[i][j], B->getBlock(i,j), coeffBeta[i][j] );
414 }
415 }
416}
417
418template <class SC, class LO, class GO, class NO>
419void BlockMatrix<SC,LO,GO,NO>::writeMM(std::string fN) const{
420 for (int i=0; i<blockMatrix_.size(); i++) {
421 for (int j=0; j<blockMatrix_.size(); j++){
422 if (!blockMatrix_[i][j].is_null() ){
423 std::string fileName = fN + std::to_string(i) + std::to_string(j) + ".mm" ;
424 blockMatrix_[i][j]->writeMM(fileName);
425 }
426 }
427 }
428}
429
430template <class SC, class LO, class GO, class NO>
431typename BlockMatrix<SC,LO,GO,NO>::MatrixPtr_Type BlockMatrix<SC,LO,GO,NO>::getMergedMatrix() {
432 if (this->size()>1) {
433 this->merge();
434 return mergedMatrix_;
435 }
436 else
437 return this->getBlock(0,0);
438}
439
440template <class SC, class LO, class GO, class NO>
441typename BlockMatrix<SC,LO,GO,NO>::BlockMapPtr_Type BlockMatrix<SC,LO,GO,NO>::getMap() {
442
443 return blockMap_;
444}
445
446}
447#endif
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5