EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_PointToBlockDiagPermute.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// EpetraExt: Epetra Extended - Linear Algebra Services Package
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44#include "Epetra_ConfigDefs.h"
47#include "Epetra_MultiVector.h"
48#include "Epetra_Import.h"
49#include "Epetra_Export.h"
50#include "Epetra_Comm.h"
51
52#include <stdio.h>
53#include <fstream>
54
55#define MAX(x,y) ((x)>(y)?(x):(y))
56
57//=========================================================================
59 :Epetra_DistObject(MAT.RowMap()),
60 Matrix_(&MAT),
61 PurelyLocalMode_(true),
62 ContiguousBlockMode_(false),
63 ContiguousBlockSize_(0),
64 NumBlocks_(0),
65 Blockstart_(0),
66 Blockids_int_(0),
67#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
68 Blockids_LL_(0),
69#endif
70 BDMap_(0),
71 CompatibleMap_(0),
72 BDMat_(0),
73 Importer_(0),
74 Exporter_(0),
75 ImportVector_(0),
76 ExportVector_(0)
77{
78
79}
80
81//=========================================================================
82// Destructor
84{
85 if(BDMap_) delete BDMap_;
86 if(CompatibleMap_) delete CompatibleMap_;
87 if(BDMat_) delete BDMat_;
88 if(Importer_) delete Importer_;
89 if(Exporter_) delete Exporter_;
90 if(ImportVector_) delete ImportVector_;
91 if(ExportVector_) delete ExportVector_;
92}
93
94//=========================================================================
95// Set list
96template<typename int_type>
97int EpetraExt_PointToBlockDiagPermute::TSetParameters(Teuchos::ParameterList & List){
98 List_=List;
99
100 // Check for contiguous blocking first
101 ContiguousBlockSize_=List_.get("contiguous block size",0);
102 if(ContiguousBlockSize_!=0){
103 ContiguousBlockMode_=true;
104 PurelyLocalMode_=false;
105 }
106
107 // Local vs. global ids & mode
108 NumBlocks_=List_.get("number of local blocks",0);
109 Blockstart_=List_.get("block start index",(int*)0);
110 Blockids_ref<int>()=List_.get("block entry lids",(int*)0);
111 if(Blockids_const_ptr<int>())
112 PurelyLocalMode_=true;
113 else{
114 Blockids_ref<int_type>()=List_.get("block entry gids",(int_type*)0);
115 if(Blockids_const_ptr<int_type>()) {
116 PurelyLocalMode_=false;
117 if(sizeof(int) != sizeof(int_type)) {
118 Blockids_ref<int>() = new int[Matrix_->RowMap().NumMyElements()];
119 }
120 }
121 }
122
123 // Sanity checks
124 if(ContiguousBlockMode_){
125 // Can't use contiguous at the same time as the other modes
126 if(NumBlocks_ || Blockstart_ || Blockids_const_ptr<int>() || Blockids_const_ptr<int_type>()) EPETRA_CHK_ERR(-4);
127 }
128 else {
129 if(NumBlocks_ <= 0) EPETRA_CHK_ERR(-1);
130 if(!Blockstart_) EPETRA_CHK_ERR(-2);
131 if(!Blockids_const_ptr<int>() && !Blockids_const_ptr<int_type>()) EPETRA_CHK_ERR(-3);
132 }
133
134 return 0;
135}
136
137int EpetraExt_PointToBlockDiagPermute::SetParameters(Teuchos::ParameterList & List){
138#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
139 if(Matrix_->RowMap().GlobalIndicesInt()) {
140 return TSetParameters<int>(List);
141 }
142 else
143#endif
144#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
145 if(Matrix_->RowMap().GlobalIndicesLongLong()) {
146 return TSetParameters<long long>(List);
147 }
148 else
149#endif
150 throw "EpetraExt_PointToBlockDiagPermute::SetParameters: ERROR, GlobalIndices type unknown.";
151}
152
153//=========================================================================
154// Extracts the block-diagonal, builds maps, etc.
156 int rv=ExtractBlockDiagonal();
157 return rv;
158}
159
160//=========================================================================
161// Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
163 return -1;
164
165}
166
167//=========================================================================
168// Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
170 // Stuff borrowed from Epetra_CrsMatrix
171 int NumVectors = X.NumVectors();
172 if (NumVectors!=Y.NumVectors()) {
173 EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
174 }
175
176 const Epetra_MultiVector *Xp=&X;
177 Epetra_MultiVector *Yp=&Y;
178
179 // Allocate temp workspace if X==Y and there are no imports or exports
180 Epetra_MultiVector * Xcopy = 0;
181 if (&X==&Y && Importer_==0 && Exporter_==0) {
182 Xcopy = new Epetra_MultiVector(X);
183 Xp=Xcopy;
184 }
185
186 UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
187 UpdateExportVector(NumVectors);
188
189 // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
190 if (Importer_){
191 EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer_, Insert));
192 Xp=ImportVector_;
193 }
194
195 // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
196 if (Exporter_) {
197 Yp=ExportVector_;
198 }
199
200 // Do the matvec
201 BDMat_->ApplyInverse(*Xp,*Yp);
202
203 // Export if needed
204 if (Exporter_) {
205 Y.PutScalar(0.0); // Make sure target is zero
206 Y.Export(*ExportVector_, *Exporter_, Add); // Fill Y with Values from export vector
207 }
208
209 // Cleanup
210 if(Xcopy) {
211 delete Xcopy;
212 EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of X
213 return 1;
214 }
215
216 return 0;
217}
218
219//=========================================================================
220// Print method
221void EpetraExt_PointToBlockDiagPermute::Print(std::ostream& /* os */) const{
222 if(Importer_) std::cout<<*Importer_<<std::endl;
223 if(Exporter_) std::cout<<*Exporter_<<std::endl;
224 if(BDMat_) std::cout<<*BDMat_<<std::endl;
225}
226
227//=========================================================================
228// Pulls the block diagonal of the matrix and then builds the BDMat_
229template<typename int_type>
230int EpetraExt_PointToBlockDiagPermute::TExtractBlockDiagonal(){
231 int i,j;
232 std::vector<int_type> ExtRows;
233 int* LocalColIDS=0;
234 int ExtSize=0,ExtCols=0,MainCols=0;
235 int Nrows=Matrix_->NumMyRows();
236 int *l2b,*block_offset;
237 int_type *l2blockid;
238 const Epetra_Map &RowMap=Matrix_->RowMap();
239 int index,col,row_in_block,col_in_block,length,*colind;
240 double *values;
241
242 bool verbose=(bool)(List_.get("output",0) > 0);
243
244 // Contiguous Setup
245 SetupContiguousMode();
246
247 // Compute block size lists
248 int *bsize=new int[NumBlocks_];
249 for(i=0;i<NumBlocks_;i++)
250 bsize[i]=Blockstart_[i+1]-Blockstart_[i];
251
252 // Use the ScanSum function to compute a prefix sum of the number of points
253 int_type MyMinGID;
254 int_type NumBlocks_int_type = NumBlocks_;
255 Matrix_->Comm().ScanSum(&NumBlocks_int_type,&MyMinGID, 1);
256 MyMinGID-=NumBlocks_;
257 int_type *MyBlockGIDs=new int_type[NumBlocks_];
258 for(i=0;i<NumBlocks_;i++)
259 MyBlockGIDs[i]=MyMinGID+i;
260
261 BDMap_=new Epetra_BlockMap((int_type) -1,NumBlocks_,MyBlockGIDs,bsize,0,Matrix_->Comm());
262
263 // Allocate the Epetra_DistBlockMatrix
264 BDMat_=new EpetraExt_BlockDiagMatrix(*BDMap_,true);
265
266 // First check to see if we can switch back to PurelyLocalMode_ if we're not
267 // in it
268 if(!PurelyLocalMode_){
269 // Find the non-local row IDs
270 ExtSize=MAX(0,Blockstart_[NumBlocks_]-RowMap.NumMyElements());// estimate
271 ExtRows.reserve(ExtSize);
272
273 for(i=0;i<Blockstart_[NumBlocks_];i++){
274 if(RowMap.LID(Blockids_const_ptr<int_type>()[i])==-1)
275 ExtRows.push_back(Blockids_const_ptr<int_type>()[i]);
276 }
277
278 // Switch to PurelyLocalMode_ if we never need GIDs
279 int size_sum;
280 ExtSize=ExtRows.size();
281 RowMap.Comm().SumAll(&ExtSize,&size_sum,1);
282 if(size_sum==0){
283 if(verbose && !Matrix_->Comm().MyPID()) printf("EpetraExt_PointToBlockDiagPermute: Switching to purely local mode\n");
284 PurelyLocalMode_=true;
285 for(i=0;i<Blockstart_[NumBlocks_];i++){
286 Blockids_ref<int>()[i]=RowMap.LID(Blockids_const_ptr<int_type>()[i]);
287 }
288 }
289 }
290
291 if(PurelyLocalMode_){
292 /*******************************************************/
293 // Allocations
294 l2b=new int[Nrows];
295 block_offset=new int[Nrows];
296
297 // Build the local-id-to-block-id list, block offset list
298 for(i=0;i<NumBlocks_;i++) {
299 for(j=Blockstart_[i];j<Blockstart_[i+1];j++){
300 block_offset[Blockids_const_ptr<int>()[j]]=j-Blockstart_[i];
301 l2b[Blockids_const_ptr<int>()[j]]=i;
302 }
303 }
304
305 // Copy the data to the EpetraExt_BlockDiagMatrix
306 for(i=0;i<Nrows;i++) {
307 int block_num = l2b[i];
308 if(block_num>=0 && block_num<NumBlocks_) {
309 row_in_block=block_offset[i];
310 Matrix_->ExtractMyRowView(i,length,values,colind);
311 int Nnz=0;
312 for (j = 0; j < length; j++) {
313 col = colind[j];
314 if(col < Nrows && l2b[col]==block_num){
315 Nnz++;
316 col_in_block = block_offset[col];
317 index = col_in_block * bsize[block_num] + row_in_block;
318 (*BDMat_)[block_num][index]=values[j];
319 }
320 }
321 /* Handle the case of a zero row. */
322 /* By just putting a 1 on the diagonal. */
323 if (Nnz == 0) {
324 index = row_in_block * bsize[block_num] + row_in_block;
325 (*BDMat_)[block_num][index]=1.0;
326 }
327 }
328 }
329
330 // Build the compatible map for import/export
331 l2blockid=new int_type[Nrows];
332 for(i=0;i<Nrows;i++)
333 l2blockid[Blockstart_[l2b[i]]+block_offset[i]]= (int_type) Matrix_->RowMap().GID64(i);
334
335 // Build the Compatible Map, Import/Export Objects
336 CompatibleMap_=new Epetra_Map(-1,Nrows,l2blockid,0,Matrix_->Comm());
337
338 }
339 else{
340 /*******************************************************/
341 // Do the import to grab matrix entries
342 // Allocate temporaries for import
343 int_type* ExtRowsPtr = ExtRows.size()>0 ? &ExtRows[0] : NULL;
344 Epetra_Map TmpMap((int_type) -1,ExtSize, ExtRowsPtr,0,Matrix_->Comm());;
345 Epetra_CrsMatrix TmpMatrix(Copy,TmpMap,0);
346 Epetra_Import TmpImporter(TmpMap,RowMap);
347
348 TmpMatrix.Import(*Matrix_,TmpImporter,Insert);
349 TmpMatrix.FillComplete();
350
351 ExtCols=TmpMatrix.NumMyCols();
352 MainCols=Matrix_->NumMyCols();
353
354
355 // Build the column reidex - main matrix
356 LocalColIDS=new int[MainCols+ExtCols];
357 for(i=0;i<MainCols;i++){
358 int_type GID= (int_type) Matrix_->GCID64(i);
359 int MainLID=RowMap.LID(GID);
360 if(MainLID!=-1) LocalColIDS[i]=MainLID;
361 else{
362 int ExtLID=TmpMatrix.LRID(GID);
363 if(ExtLID!=-1) LocalColIDS[i]=Nrows+ExtLID;
364 else LocalColIDS[i]=-1;
365 }
366 }
367
368 // Build the column reidex - ext matrix
369 for(i=0;i<ExtCols;i++){
370 int_type GID= (int_type) TmpMatrix.GCID64(i);
371 int MainLID=RowMap.LID(GID);
372 if(MainLID!=-1) LocalColIDS[MainCols+i]=MainLID;
373 else{
374 int ExtLID=TmpMatrix.LRID(GID);
375 if(ExtLID!=-1) LocalColIDS[MainCols+i]=Nrows+ExtLID;
376 else LocalColIDS[MainCols+i]=-1;
377 }
378 }
379
380 // Allocations
381 l2b=new int[Nrows+ExtSize];
382 block_offset=new int[Nrows+ExtSize];
383
384 // Build l2b/block_offset with the expanded local index
385 //NTS: Uses the same ordering of operation as the above routine, which is why
386 //it works.
387 for(i=0;i<Nrows+ExtSize;i++) block_offset[i]=l2b[i]=-1;
388 int ext_idx=0;
389 for(i=0;i<NumBlocks_;i++) {
390 for(j=Blockstart_[i];j<Blockstart_[i+1];j++){
391 int LID=RowMap.LID(Blockids_const_ptr<int_type>()[j]);
392 if(LID==-1) {LID=Nrows+ext_idx;ext_idx++;}
393 block_offset[LID]=j-Blockstart_[i];
394 l2b[LID]=i;
395 }
396 }
397
398 // Copy the data to the EpetraExt_BlockDiagMatrix from Matrix_
399 for(i=0;i<Nrows;i++) {
400 int block_num = l2b[i];
401 if(block_num>=0 && block_num<NumBlocks_) {
402 row_in_block=block_offset[i];
403 Matrix_->ExtractMyRowView(i,length,values,colind);
404 int Nnz=0;
405 for (j = 0; j < length; j++) {
406 col = LocalColIDS[colind[j]];
407 if(col!=-1 && l2b[col]==block_num){
408 Nnz++;
409 col_in_block = block_offset[col];
410 index = col_in_block * bsize[block_num] + row_in_block;
411 (*BDMat_)[block_num][index]=values[j];
412 }
413 }
414 /* Handle the case of a zero row. */
415 /* By just putting a 1 on the diagonal. */
416 if (Nnz == 0) {
417 index = row_in_block * bsize[block_num] + row_in_block;
418 (*BDMat_)[block_num][index]=values[j];
419 }
420 }
421 }
422
423 // Copy the data to the EpetraExt_BlockDiagMatrix from TmpMatrix
424 for(i=0;i<ExtSize;i++) {
425 int block_num = l2b[Nrows+i];
426 if(block_num>=0 && block_num<NumBlocks_) {
427 row_in_block=block_offset[Nrows+i];
428 TmpMatrix.ExtractMyRowView(i,length,values,colind);
429 int Nnz=0;
430 for (j = 0; j < length; j++) {
431 col = LocalColIDS[MainCols+colind[j]];
432 if(col!=-1 && l2b[col]==block_num){
433 Nnz++;
434 col_in_block = block_offset[col];
435 index = col_in_block * bsize[block_num] + row_in_block;
436 (*BDMat_)[block_num][index]=values[j];
437 }
438 }
439 /* Handle the case of a zero row. */
440 /* By just putting a 1 on the diagonal. */
441 if (Nnz == 0) {
442 index = row_in_block * bsize[block_num] + row_in_block;
443 (*BDMat_)[block_num][index]=1.0;
444 }
445 }
446 }
447
448 // Build the compatible map for import/export
449 l2blockid=new int_type[Blockstart_[NumBlocks_]];
450 for(i=0;i<Nrows;i++){
451 int bkid=l2b[i];
452 if(bkid>-1) l2blockid[Blockstart_[bkid]+block_offset[i]]= (int_type) RowMap.GID64(i);
453 }
454 // NTS: This is easier - if we imported it, we know we need it
455 for(i=0;i<ExtSize;i++)
456 l2blockid[Blockstart_[l2b[Nrows+i]]+block_offset[Nrows+i]]= (int_type) TmpMatrix.GRID64(i);
457
458 // Build the Compatible Map, Import/Export Objects
459 CompatibleMap_=new Epetra_Map(-1,Blockstart_[NumBlocks_],l2blockid,0,Matrix_->Comm());
460
461 }//end else
462
463 // Set BDMat_'s parameter list and compute!
464 Teuchos::ParameterList dummy,inList;
465 inList=List_.get("blockdiagmatrix: list",dummy);
466 BDMat_->SetParameters(inList);
467 BDMat_->Compute();
468
469 // Build importer/exporter
470 if(!CompatibleMap_->SameAs(Matrix_->DomainMap()))
471 Importer_ = new Epetra_Import(*CompatibleMap_,Matrix_->DomainMap());
472 if(!CompatibleMap_->SameAs(Matrix_->RangeMap()))
473 Exporter_ = new Epetra_Export(*CompatibleMap_,Matrix_->RangeMap());
474
475 // Cleanup
476 delete [] LocalColIDS;
477 delete [] block_offset;
478 delete [] l2b;
479 delete [] l2blockid;
480 delete [] bsize;
481 delete [] MyBlockGIDs;
482
483 // Contiguous Cleanup
484 CleanupContiguousMode();
485
486 return 0;
487}
488
489int EpetraExt_PointToBlockDiagPermute::ExtractBlockDiagonal(){
490#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
491 if(Matrix_->RowMap().GlobalIndicesInt()) {
492 return TExtractBlockDiagonal<int>();
493 }
494 else
495#endif
496#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
497 if(Matrix_->RowMap().GlobalIndicesLongLong()) {
498 return TExtractBlockDiagonal<long long>();
499 }
500 else
501#endif
502 throw "EpetraExt_PointToBlockDiagPermute::ExtractBlockDiagonal: ERROR, GlobalIndices type unknown.";
503}
504
505//=======================================================================================================
506template<typename int_type>
507int EpetraExt_PointToBlockDiagPermute::TSetupContiguousMode(){
508 if(!ContiguousBlockMode_) return 0;
509 // NTS: In case of processor-crossing blocks, the lowest PID always gets the block;
510 const Epetra_Map &RowMap=Matrix_->RowMap();
511
512 int_type MinMyGID= (int_type) RowMap.MinMyGID64();
513 int_type MaxMyGID= (int_type) RowMap.MaxMyGID64();
514 int_type Base= (int_type) Matrix_->IndexBase64();
515
516 // Find the GID that begins my first block
517 int_type MyFirstBlockGID=ContiguousBlockSize_*(int_type)ceil(((double)(MinMyGID - Base))/ContiguousBlockSize_)+Base;
518 NumBlocks_=(int)ceil((double)((MaxMyGID-MyFirstBlockGID+1.0)) / ContiguousBlockSize_);
519
520 // Allocate memory
521 Blockstart_=new int[NumBlocks_+1];
522 Blockids_ref<int_type>()=new int_type[NumBlocks_*ContiguousBlockSize_];
523 if(sizeof(int) != sizeof(int_type))
524 Blockids_ref<int>()=new int[NumBlocks_*ContiguousBlockSize_];
525 Blockstart_[NumBlocks_]=NumBlocks_*ContiguousBlockSize_;
526
527 // Fill the arrays
528 for(int i=0,ct=0;i<NumBlocks_;i++){
529 Blockstart_[i]=ct;
530 for(int j=0;j<ContiguousBlockSize_;j++,ct++){
531 Blockids_ref<int_type>()[ct]=MyFirstBlockGID+ct;
532 }
533 }
534
535 return 0;
536}
537
538int EpetraExt_PointToBlockDiagPermute::SetupContiguousMode(){
539#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
540 if(Matrix_->RowMap().GlobalIndicesInt()) {
541 return TSetupContiguousMode<int>();
542 }
543 else
544#endif
545#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
546 if(Matrix_->RowMap().GlobalIndicesLongLong()) {
547 return TSetupContiguousMode<long long>();
548 }
549 else
550#endif
551 throw "EpetraExt_PointToBlockDiagPermute::SetupContiguousMode: ERROR, GlobalIndices type unknown.";
552}
553
554//=======================================================================================================
555int EpetraExt_PointToBlockDiagPermute::CleanupContiguousMode(){
556 if(!ContiguousBlockMode_) return 0;
557 NumBlocks_=0;
558 if(Blockstart_) {delete [] Blockstart_; Blockstart_=0;}
559 if(Blockids_int_) {delete [] Blockids_int_; Blockids_int_=0;}
560#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
561 if(Blockids_LL_) {delete [] Blockids_LL_; Blockids_LL_=0;}
562#endif
563 return 0;
564}
565
566
567//=======================================================================================================
568// Creates an Epetra_CrsMatrix from the BlockDiagMatrix. This is generally only useful if you want to do a matrix-matrix multiply.
569template<typename int_type>
570Epetra_FECrsMatrix * EpetraExt_PointToBlockDiagPermute::TCreateFECrsMatrix(){
571 Epetra_FECrsMatrix * NewMat=new Epetra_FECrsMatrix(Copy,Matrix_->RowMap(),0);
572
573 const Epetra_BlockMap &BlockMap=BDMat_->BlockMap();
574 const Epetra_BlockMap &DataMap=BDMat_->DataMap();
575 const int *vlist=DataMap.FirstPointInElementList();
576 const int *xlist=BlockMap.FirstPointInElementList();
577 const int *blocksize=BlockMap.ElementSizeList();
578 const double *values=BDMat_->Values();
579 int NumBlocks=BDMat_->NumMyBlocks();
580
581 // Maximum size vector for stashing GIDs
582 std::vector<int_type> GIDs;
583 GIDs.resize(BlockMap.MaxMyElementSize());
584
585
586 for(int i=0;i<NumBlocks;i++){
587 int Nb=blocksize[i];
588 int vidx0=vlist[i];
589 int xidx0=xlist[i];
590 // Get global indices
591 for(int j=0;j<Nb;j++)
592 GIDs[j]= (int_type) CompatibleMap_->GID64(xidx0+j);
593
594 // Remember: We're using column-major storage for LAPACK's benefit
595 int ierr=NewMat->InsertGlobalValues(Nb,&GIDs[0],&values[vidx0],Epetra_FECrsMatrix::COLUMN_MAJOR);
596 if(ierr < 0) throw "CreateFECrsMatrix: ERROR in InsertGlobalValues";
597 }
598 NewMat->GlobalAssemble();
599 return NewMat;
600}
601
603#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
604 if(Matrix_->RowMap().GlobalIndicesInt()) {
605 return TCreateFECrsMatrix<int>();
606 }
607 else
608#endif
609#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
610 if(Matrix_->RowMap().GlobalIndicesLongLong()) {
611 return TCreateFECrsMatrix<long long>();
612 }
613 else
614#endif
615 throw "EpetraExt_PointToBlockDiagPermute::CreateFECrsMatrix: ERROR, GlobalIndices type unknown.";
616}
617
618//=======================================================================================================
619void EpetraExt_PointToBlockDiagPermute::UpdateImportVector(int NumVectors) const {
620 if(Importer_ != 0) {
621 if(ImportVector_ != 0) {
622 if(ImportVector_->NumVectors() != NumVectors) {
623 delete ImportVector_;
624 ImportVector_= 0;
625 }
626 }
627 if(ImportVector_ == 0)
628 ImportVector_ = new Epetra_MultiVector(*CompatibleMap_,NumVectors); // Create import vector if needed
629 }
630 return;
631}
632//=======================================================================================================
633void EpetraExt_PointToBlockDiagPermute::UpdateExportVector(int NumVectors) const {
634 if(Exporter_ != 0) {
635 if(ExportVector_ != 0) {
636 if(ExportVector_->NumVectors() != NumVectors) {
637 delete ExportVector_;
638 ExportVector_= 0;
639 }
640 }
641 if(ExportVector_ == 0)
642 ExportVector_ = new Epetra_MultiVector(*CompatibleMap_,NumVectors); // Create Export vector if needed
643 }
644 return;
645}
646
647
648
649
650//=========================================================================
651int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& /* A */, const Epetra_Import& /* Importer */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
652 return -1;
653}
654
655//=========================================================================
656int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& /* A */, const Epetra_Export& /* Exporter */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
657 return -1;
658}
659
660//=========================================================================
661int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& /* A */, const Epetra_Import & /* Importer */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
662 return -1;
663}
664
665//=========================================================================
666int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& /* A */, const Epetra_Export& /* Exporter */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
667 return -1;
668}
669
670//=========================================================================
671// Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
673 return -1;
674}
675
676//=========================================================================
677// Perform ID copies and permutations that are on processor.
679 int /* NumSameIDs */,
680 int /* NumPermuteIDs */,
681 int * /* PermuteToLIDs */,
682 int * /* PermuteFromLIDs */,
683 const Epetra_OffsetIndex * /* Indexor */,
684 Epetra_CombineMode /* CombineMode */){
685 return -1;
686}
687
688//=========================================================================
689// Perform any packing or preparation required for call to DoTransfer().
691 int /* NumExportIDs */,
692 int* /* ExportLIDs */,
693 int& /* LenExports */,
694 char*& /* Exports */,
695 int& /* SizeOfPacket */,
696 int* /* Sizes */,
697 bool & /* VarSizes */,
698 Epetra_Distributor& /* Distor */){
699 return -1;
700}
701
702//=========================================================================
703// Perform any unpacking and combining after call to DoTransfer().
705 int /* NumImportIDs */,
706 int* /* ImportLIDs */,
707 int /* LenImports */,
708 char* /* Imports */,
709 int& /* SizeOfPacket */,
710 Epetra_Distributor& /* Distor */,
711 Epetra_CombineMode /* CombineMode */,
712 const Epetra_OffsetIndex * /* Indexor */){
713 return -1;
714}
715
716
Epetra_CombineMode
Insert
Add
EpetraExt_BlockDiagMatrix: A class for storing distributed block matrices.
double * Values() const
Returns a pointer to the array containing the blocks.
virtual int SetParameters(Teuchos::ParameterList &List)
SetParameters.
virtual int Compute()
Computes the inverse / factorization if such is set on the list.
virtual const Epetra_BlockMap & BlockMap() const
Returns the Epetra_BlockMap object associated with the range of this operator.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
virtual const Epetra_BlockMap & DataMap() const
Returns the Epetra_BlockMap object with the distribution of underlying values.
int NumMyBlocks() const
Returns the number of local blocks.
virtual void Print(std::ostream &os) const
Print information about this object to the given output stream.
virtual int PackAndPrepare(const Epetra_SrcDistObject &Source, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
Perform any packing or preparation required for call to DoTransfer().
virtual Epetra_FECrsMatrix * CreateFECrsMatrix()
Create an Epetra_FECrsMatrix from the BlockDiagMatrix.
EpetraExt_PointToBlockDiagPermute(const Epetra_CrsMatrix &MAT)
@ Name Constructors
virtual int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
Perform ID copies and permutations that are on processor.
virtual int UnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
Perform any unpacking and combining after call to DoTransfer().
virtual int SetParameters(Teuchos::ParameterList &List)
Sets the parameter list.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
virtual int Compute()
Extracts the block-diagonal, builds maps, etc.
virtual int CheckSizes(const Epetra_SrcDistObject &Source)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
int * ElementSizeList() const
int LID(int GID) const
int * FirstPointInElementList() const
int MaxMyElementSize() const
bool SameAs(const Epetra_BlockMap &Map) const
const Epetra_Comm & Comm() const
int NumMyElements() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int NumVectors() const