EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
96 template<typename int_type>
97 int 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 
137 int 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
221 void 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_
229 template<typename int_type>
230 int 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 
489 int 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 //=======================================================================================================
506 template<typename int_type>
507 int 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 
538 int 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 //=======================================================================================================
555 int 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.
569 template<typename int_type>
570 Epetra_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 //=======================================================================================================
619 void 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 //=======================================================================================================
633 void 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 //=========================================================================
651 int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& /* A */, const Epetra_Import& /* Importer */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
652  return -1;
653 }
654 
655 //=========================================================================
656 int EpetraExt_PointToBlockDiagPermute::Import(const Epetra_SrcDistObject& /* A */, const Epetra_Export& /* Exporter */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
657  return -1;
658 }
659 
660 //=========================================================================
661 int EpetraExt_PointToBlockDiagPermute::Export(const Epetra_SrcDistObject& /* A */, const Epetra_Import & /* Importer */, Epetra_CombineMode /* CombineMode */, const Epetra_OffsetIndex * /* Indexor */){
662  return -1;
663 }
664 
665 //=========================================================================
666 int 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 
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...
const Epetra_Map & RangeMap() const
virtual int SetParameters(Teuchos::ParameterList &List)
Sets the parameter list.
#define MAX(x, y)
bool GlobalIndicesLongLong() const
bool SameAs(const Epetra_BlockMap &Map) const
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...
bool GlobalIndicesInt() const
int * FirstPointInElementList() const
virtual int MyPID() const =0
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
EpetraExt_BlockDiagMatrix: A class for storing distributed block matrices.
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().
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
const Epetra_Map & RowMap() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
int NumMyCols() const
virtual int CheckSizes(const Epetra_SrcDistObject &Source)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
int NumMyRows() const
int GlobalAssemble(bool callFillComplete=true, Epetra_CombineMode combineMode=Add, bool save_off_and_reuse_map_exporter=false)
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.
int LID(int GID) const
const Epetra_Comm & Comm() const
EpetraExt_PointToBlockDiagPermute(const Epetra_CrsMatrix &MAT)
@ Name Constructors
virtual int Compute()
Extracts the block-diagonal, builds maps, etc.
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.
const Epetra_Map & DomainMap() const
Epetra_CombineMode
virtual void Print(std::ostream &os) const
Print information about this object to the given output stream.
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.
const Epetra_Comm & Comm() const
virtual int ScanSum(double *MyVals, double *ScanSums, int Count) const =0
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 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.