EpetraExt  Development
EpetraExt_BlockDiagMatrix.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 
45 #include "Epetra_MultiVector.h"
46 #include "Epetra_Comm.h"
47 #include "Epetra_LAPACK.h"
48 #include "Epetra_Distributor.h"
49 
50 #define AM_MULTIPLY 0
51 #define AM_INVERT 1
52 #define AM_FACTOR 2
53 
54 //=========================================================================
55 // Constructor
56 EpetraExt_BlockDiagMatrix::EpetraExt_BlockDiagMatrix(const Epetra_BlockMap& map,bool zero_out)
57  : Epetra_DistObject(map, "EpetraExt::BlockDiagMatrix"),
58  HasComputed_(false),
59  ApplyMode_(AM_MULTIPLY),
60  DataMap_(0),
61  Values_(0),
62  Pivots_(0)
63 {
64  Allocate();
65  if(zero_out) PutScalar(0.0);
66 }
67 
68 
69 //=========================================================================
70 // Destructor
72  if(DataMap_) delete DataMap_;
73  if(Pivots_) delete [] Pivots_;
74  if(Values_) delete [] Values_;
75 }
76 
77 
78 //=========================================================================
79 // Copy constructor.
81  : Epetra_DistObject(Source),
82  HasComputed_(false),
83  ApplyMode_(AM_MULTIPLY),
84  Values_(0),
85  Pivots_(0)
86 {
87  DataMap_=new Epetra_BlockMap(*Source.DataMap_);
88  Pivots_=new int[NumMyUnknowns()];
89  Values_=new double[DataMap_->NumMyPoints()];
90  DoCopy(Source);
91 }
92 
93 //=========================================================================
94 // Allocate
95 void EpetraExt_BlockDiagMatrix::Allocate(){
96 
97  int dataSize=0, NumBlocks=NumMyBlocks();
98  Pivots_=new int[NumMyUnknowns()];
99  int *ElementSize=new int[NumBlocks];
100 
101  for(int i=0;i<NumBlocks;i++) {
102  ElementSize[i]=BlockSize(i)*BlockSize(i);
103  dataSize+=ElementSize[i];
104  }
105 
106 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
107  if(Map().GlobalIndicesInt()) {
108  DataMap_=new Epetra_BlockMap(-1,Map().NumMyElements(),Map().MyGlobalElements(),ElementSize,0,Map().Comm());
109  }
110  else
111 #endif
112 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
113  if(Map().GlobalIndicesLongLong()) {
114  DataMap_=new Epetra_BlockMap((long long) -1,Map().NumMyElements(),Map().MyGlobalElements64(),ElementSize,0,Map().Comm());
115  }
116  else
117 #endif
118  throw "EpetraExt_BlockDiagMatrix::Allocate: GlobalIndices type unknown";
119 
120  Values_=new double[dataSize];
121  delete [] ElementSize;
122 }
123 
124 
125 //=========================================================================
127 int EpetraExt_BlockDiagMatrix::SetParameters(Teuchos::ParameterList & List){
128  List_=List;
129 
130  // Inverse storage mode
131  std::string dummy("multiply");
132  std::string ApplyMode=List_.get("apply mode",dummy);
133  if(ApplyMode==std::string("multiply")) ApplyMode_=AM_MULTIPLY;
134  else if(ApplyMode==std::string("invert")) ApplyMode_=AM_INVERT;
135  else if(ApplyMode==std::string("factor")) ApplyMode_=AM_FACTOR;
136  else EPETRA_CHK_ERR(-1);
137 
138  return 0;
139 }
140 
141 
142 //=========================================================================
144  int MaxData=NumData();
145  for (int i=0;i<MaxData;i++) Values_[i]=value;
146 }
147 
148 //=========================================================================
149 // Assignment operator: Needs the same maps
151  DoCopy(Source);
152  return(*this);
153 }
154 //=========================================================================
155 int EpetraExt_BlockDiagMatrix::DoCopy(const EpetraExt_BlockDiagMatrix& Source){
156  // Need the same map
157  if(!Map().SameAs(Source.Map()) || !DataMap_->SameAs(*Source.DataMap_))
158  throw ReportError("Maps of DistBlockMatrices incompatible in assignment",-1);
159 
160  int MaxData=Source.NumData();
161 
162  for(int i=0;i<MaxData;i++) Values_[i]=Source.Values_[i];
163  for(int i=0;i<Source.NumMyUnknowns();i++) Pivots_[i]=Source.Pivots_[i];
164 
165  List_=Source.List_;
166  ApplyMode_=Source.ApplyMode_;
167  HasComputed_=Source.HasComputed_;
168 
169  return 0;
170 }
171 
172 
173 //=========================================================================
176  int info;
177 
178  if(ApplyMode_==AM_MULTIPLY)
179  // Multiply mode - noop
180  return 0;
181  else {
182  // Factorization - Needed for both 'factor' and 'invert' modes
183  int NumBlocks=NumMyBlocks();
184  for(int i=0;i<NumBlocks;i++){
185  int Nb=BlockSize(i);
186  if(Nb==1) {
187  // Optimize for Block Size 1
188  Values_[DataMap_->FirstPointInElement(i)]=1.0/Values_[DataMap_->FirstPointInElement(i)];
189  }
190  else if(Nb==2) {
191  // Optimize for Block Size 2
192  double * v=&Values_[DataMap_->FirstPointInElement(i)];
193  double d=1/(v[0]*v[3]-v[1]*v[2]);
194  double v0old=v[0];
195  v[0]=v[3]*d;
196  v[1]=-v[1]*d;
197  v[2]=-v[2]*d;
198  v[3]=v0old*d;
199  }
200  else{
201  // "Large" Block - Use LAPACK
202  LAPACK.GETRF(Nb,Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&info);
203  if(info) EPETRA_CHK_ERR(-2);
204  }
205  }
206 
207  if(ApplyMode_==AM_INVERT){
208  // Invert - Use the factorization and invert the blocks in-place
209  int lwork=3*DataMap_->MaxMyElementSize();
210  std::vector<double> work(lwork);
211  for(int i=0;i<NumBlocks;i++){
212  int Nb=BlockSize(i);
213  if(Nb==1 || Nb==2){
214  // Optimize for Block Size 1 and 2
215  // No need to do anything - factorization has already inverted the value
216  }
217  else{
218  // "Large" Block - Use LAPACK
219  LAPACK.GETRI(Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&work[0],&lwork,&info);
220  if(info) EPETRA_CHK_ERR(-3);
221  }
222  }
223  }
224  }
225  HasComputed_=true;
226  return 0;
227 }
228 
229 
230 //=========================================================================
231 int EpetraExt_BlockDiagMatrix::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
232  int info;
233  // Sanity Checks
234  int NumVectors=X.NumVectors();
235  if(NumVectors!=Y.NumVectors())
236  EPETRA_CHK_ERR(-1);
237  if(!HasComputed_ && (ApplyMode_==AM_INVERT || ApplyMode_==AM_FACTOR))
238  EPETRA_CHK_ERR(-2);
239 
240  //NTS: MultiVector's MyLength and [] Operators are "points" level operators
241  //not a "block/element" level operators.
242 
243  const int *vlist=DataMap_->FirstPointInElementList();
244  const int *xlist=Map().FirstPointInElementList();
245  const int *blocksize=Map().ElementSizeList();
246 
247  if(ApplyMode_==AM_MULTIPLY || ApplyMode_==AM_INVERT){
248  // Multiply & Invert mode have the same apply
249  int NumBlocks=NumMyBlocks();
250  for(int i=0;i<NumBlocks;i++){
251  int Nb=blocksize[i];
252  int vidx0=vlist[i];
253  int xidx0=xlist[i];
254  for(int j=0;j<NumVectors;j++){
255  if(Nb==1) {
256  // Optimize for size = 1
257  Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
258  }
259  else if(Nb==2){
260  // Optimize for size = 2
261  Y[j][xidx0 ]=Values_[vidx0 ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
262  Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
263  }
264  else{
265  // "Large" Block - Use BLAS
266  //void GEMV (const char TRANS, const int M, const int N, const double ALPHA, const double *A, const int LDA, const double *X, const double BETA, double *Y, const int INCX=1, const int INCY=1) const
267  GEMV('N',Nb,Nb,1.0,&Values_[vidx0],Nb,&X[j][xidx0],0.0,&Y[j][xidx0]);
268  }
269  }
270  }
271  }
272  else{
273  // Factorization mode has a different apply
274  int NumBlocks=NumMyBlocks();
275  for(int i=0;i<NumBlocks;i++){
276  int Nb=blocksize[i];
277  int vidx0=vlist[i];
278  int xidx0=xlist[i];
279  for(int j=0;j<NumVectors;j++){
280  if(Nb==1) {
281  // Optimize for size = 1 - use the inverse
282  Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
283  }
284  else if(Nb==2){
285  // Optimize for size = 2 - use the inverse
286  Y[j][xidx0 ]=Values_[vidx0 ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
287  Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
288  }
289  else{
290  // "Large" Block - use LAPACK
291  // void GETRS (const char TRANS, const int N, const int NRHS, const double *A, const int LDA, const int *IPIV, double *X, const int LDX, int *INFO) const
292  for(int k=0;k<Nb;k++) Y[j][xidx0+k]=X[j][xidx0+k];
293  LAPACK.GETRS('N',Nb,1,&Values_[vidx0],Nb,&Pivots_[xidx0],&Y[j][xidx0],Nb,&info);
294  if(info) EPETRA_CHK_ERR(info);
295  }
296  }
297  }
298  }
299  return 0;
300 }
301 
302 
303 
304 
305 //=========================================================================
306 // Print method
307 void EpetraExt_BlockDiagMatrix::Print(std::ostream & os) const{
308  int MyPID = DataMap_->Comm().MyPID();
309  int NumProc = DataMap_->Comm().NumProc();
310 
311  for (int iproc=0; iproc < NumProc; iproc++) {
312  if (MyPID==iproc) {
313  int NumMyElements1 =DataMap_->NumMyElements();
314  int MaxElementSize1 = DataMap_->MaxElementSize();
315  int * MyGlobalElements1_int = 0;
316  long long * MyGlobalElements1_LL = 0;
317 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
318  if (DataMap_->GlobalIndicesInt ()) {
319  MyGlobalElements1_int = DataMap_->MyGlobalElements();
320  }
321  else
322 #endif
323 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
324  if (DataMap_->GlobalIndicesLongLong ()) {
325  MyGlobalElements1_LL = DataMap_->MyGlobalElements64();
326  }
327  else
328 #endif
329  throw "EpetraExt_BlockDiagMatrix::Print: GlobalIndices type unknown";
330 
331  int * FirstPointInElementList1 = (MaxElementSize1 == 1) ?
332  NULL : DataMap_->FirstPointInElementList();
333 
334  if (MyPID==0) {
335  os.width(8);
336  os << " MyPID"; os << " ";
337  os.width(12);
338  if (MaxElementSize1==1)
339  os << "GID ";
340  else
341  os << " GID/Point";
342  os.width(20);
343  os << "Values ";
344  os << std::endl;
345  }
346  for (int i=0; i < NumMyElements1; i++) {
347  for (int ii=0; ii < DataMap_->ElementSize(i); ii++) {
348  int iii;
349  os.width(10);
350  os << MyPID; os << " ";
351  os.width(10);
352  if (MaxElementSize1==1) {
353  os << (MyGlobalElements1_int ? MyGlobalElements1_int[i] : MyGlobalElements1_LL[i]) << " ";
354  iii = i;
355  }
356  else {
357  os << (MyGlobalElements1_int ? MyGlobalElements1_int[i] : MyGlobalElements1_LL[i]) << "/" << ii << " ";
358  iii = FirstPointInElementList1[i]+ii;
359  }
360  os.width(20);
361  os << Values_[iii];
362  os << std::endl;
363  }
364  }
365  os << std::flush;
366  }
367 
368  // Do a few global ops to give I/O a chance to complete
369  Map().Comm().Barrier();
370  Map().Comm().Barrier();
371  Map().Comm().Barrier();
372  }
373  return;
374 }
375 
376 
377 //=========================================================================
378 // Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
379 int EpetraExt_BlockDiagMatrix::CheckSizes(const Epetra_SrcDistObject& Source){
380  return &Map() == &Source.Map();
381 }
382 
383 
384  //=========================================================================
385 // Perform ID copies and permutations that are on processor.
386 int EpetraExt_BlockDiagMatrix::CopyAndPermute(const Epetra_SrcDistObject& Source,
387  int NumSameIDs,
388  int NumPermuteIDs,
389  int * PermuteToLIDs,
390  int * PermuteFromLIDs,
391  const Epetra_OffsetIndex * Indexor,
392  Epetra_CombineMode CombineMode){
393  (void)Indexor;
394 
395  const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
396 
397  double *From=A.Values();
398  double *To = Values_;
399 
400  int * ToFirstPointInElementList = 0;
401  int * FromFirstPointInElementList = 0;
402  int * FromElementSizeList = 0;
403  int MaxElementSize = DataMap().MaxElementSize();
404  bool ConstantElementSize = DataMap().ConstantElementSize();
405 
406  if (!ConstantElementSize) {
407  ToFirstPointInElementList = DataMap().FirstPointInElementList();
408  FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
409  FromElementSizeList = A.DataMap().ElementSizeList();
410  }
411 
412  int NumSameEntries;
413 
414  bool Case1 = false;
415  bool Case2 = false;
416  // bool Case3 = false;
417 
418  if (MaxElementSize==1) {
419  Case1 = true;
420  NumSameEntries = NumSameIDs;
421  }
422  else if (ConstantElementSize) {
423  Case2 = true;
424  NumSameEntries = NumSameIDs * MaxElementSize;
425  }
426  else {
427  // Case3 = true;
428  NumSameEntries = FromFirstPointInElementList[NumSameIDs];
429  }
430 
431  // Short circuit for the case where the source and target vector is the same.
432  if (To==From) {
433  NumSameEntries = 0;
434  }
435 
436  // Do copy first
437  if (NumSameIDs>0)
438  if (To!=From) {
439  if (CombineMode==Epetra_AddLocalAlso) {
440  for (int j=0; j<NumSameEntries; j++) {
441  To[j] += From[j]; // Add to existing value
442  }
443  } else {
444  for (int j=0; j<NumSameEntries; j++) {
445  To[j] = From[j];
446  }
447  }
448  }
449  // Do local permutation next
450  if (NumPermuteIDs>0) {
451 
452  // Point entry case
453  if (Case1) {
454 
455  if (CombineMode==Epetra_AddLocalAlso) {
456  for (int j=0; j<NumPermuteIDs; j++) {
457  To[PermuteToLIDs[j]] += From[PermuteFromLIDs[j]]; // Add to existing value
458  }
459  }
460  else {
461  for (int j=0; j<NumPermuteIDs; j++) {
462  To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
463  }
464  }
465  }
466  // constant element size case
467  else if (Case2) {
468  for (int j=0; j<NumPermuteIDs; j++) {
469  int jj = MaxElementSize*PermuteToLIDs[j];
470  int jjj = MaxElementSize*PermuteFromLIDs[j];
471  if (CombineMode==Epetra_AddLocalAlso) {
472  for (int k=0; k<MaxElementSize; k++) {
473  To[jj+k] += From[jjj+k]; // Add to existing value
474  }
475  }
476  else {
477  for (int k=0; k<MaxElementSize; k++) {
478  To[jj+k] = From[jjj+k];
479  }
480  }
481  }
482  }
483 
484  // variable element size case
485  else {
486  for (int j=0; j<NumPermuteIDs; j++) {
487  int jj = ToFirstPointInElementList[PermuteToLIDs[j]];
488  int jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
489  int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
490  if (CombineMode==Epetra_AddLocalAlso) {
491  for (int k=0; k<ElementSize; k++) {
492  To[jj+k] += From[jjj+k]; // Add to existing value
493  }
494  }
495  else {
496  for (int k=0; k<ElementSize; k++) {
497  To[jj+k] = From[jjj+k];
498  }
499  }
500  }
501  }
502  }
503  return(0);
504 }
505 
506 //=========================================================================
507 // Perform any packing or preparation required for call to DoTransfer().
508 int EpetraExt_BlockDiagMatrix::PackAndPrepare(const Epetra_SrcDistObject& Source,
509  int NumExportIDs,
510  int* ExportLIDs,
511  int& LenExports,
512  char*& Exports,
513  int& SizeOfPacket,
514  int* Sizes,
515  bool & VarSizes,
516  Epetra_Distributor& Distor){
517  (void)Sizes;
518  (void)VarSizes;
519  (void)Distor;
520  const EpetraExt_BlockDiagMatrix & A = dynamic_cast<const EpetraExt_BlockDiagMatrix &>(Source);
521 
522  int j, jj, k;
523 
524  double *From=A.Values();
525  int MaxElementSize = DataMap().MaxElementSize();
526  bool ConstantElementSize = DataMap().ConstantElementSize();
527 
528  int * FromFirstPointInElementList = 0;
529  int * FromElementSizeList = 0;
530 
531  if (!ConstantElementSize) {
532  FromFirstPointInElementList = A.DataMap().FirstPointInElementList();
533  FromElementSizeList = A.DataMap().ElementSizeList();
534  }
535 
536  SizeOfPacket = MaxElementSize * (int)sizeof(double);
537 
538  if(NumExportIDs*SizeOfPacket>LenExports) {
539  if (LenExports>0) delete [] Exports;
540  LenExports = NumExportIDs*SizeOfPacket;
541  Exports = new char[LenExports];
542  }
543 
544  double * ptr;
545 
546  if (NumExportIDs>0) {
547  ptr = (double *) Exports;
548 
549  // Point entry case
550  if (MaxElementSize==1) for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];
551 
552  // constant element size case
553  else if (ConstantElementSize) {
554  for (j=0; j<NumExportIDs; j++) {
555  jj = MaxElementSize*ExportLIDs[j];
556  for (k=0; k<MaxElementSize; k++)
557  *ptr++ = From[jj+k];
558  }
559  }
560 
561  // variable element size case
562  else {
563  SizeOfPacket = MaxElementSize;
564  for (j=0; j<NumExportIDs; j++) {
565  ptr = (double *) Exports + j*SizeOfPacket;
566  jj = FromFirstPointInElementList[ExportLIDs[j]];
567  int ElementSize = FromElementSizeList[ExportLIDs[j]];
568  for (k=0; k<ElementSize; k++)
569  *ptr++ = From[jj+k];
570  }
571  }
572  }
573 
574  return(0);
575 }
576 
577 
578 //=========================================================================
579 // Perform any unpacking and combining after call to DoTransfer().
580 int EpetraExt_BlockDiagMatrix::UnpackAndCombine(const Epetra_SrcDistObject& Source,
581  int NumImportIDs,
582  int* ImportLIDs,
583  int LenImports,
584  char* Imports,
585  int& SizeOfPacket,
586  Epetra_Distributor& Distor,
587  Epetra_CombineMode CombineMode,
588  const Epetra_OffsetIndex * Indexor){
589  (void)Source;
590  (void)LenImports;
591  (void)Distor;
592  (void)Indexor;
593  int j, jj, k;
594 
595  if( CombineMode != Add
596  && CombineMode != Zero
597  && CombineMode != Insert
598  && CombineMode != Average
599  && CombineMode != AbsMax )
600  EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero
601 
602  if (NumImportIDs<=0) return(0);
603 
604  double * To = Values_;
605  int MaxElementSize = DataMap().MaxElementSize();
606  bool ConstantElementSize = DataMap().ConstantElementSize();
607 
608  int * ToFirstPointInElementList = 0;
609  int * ToElementSizeList = 0;
610 
611  if (!ConstantElementSize) {
612  ToFirstPointInElementList = DataMap().FirstPointInElementList();
613  ToElementSizeList = DataMap().ElementSizeList();
614  }
615 
616  double * ptr;
617  // Unpack it...
618 
619  ptr = (double *) Imports;
620 
621  // Point entry case
622  if (MaxElementSize==1) {
623 
624  if (CombineMode==Add)
625  for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += *ptr++; // Add to existing value
626  else if(CombineMode==Insert)
627  for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = *ptr++;
628  else if(CombineMode==AbsMax)
629  for (j=0; j<NumImportIDs; j++) {
630  To[ImportLIDs[j]] = EPETRA_MAX( To[ImportLIDs[j]],std::abs(*ptr));
631  ptr++;
632  }
633  // Note: The following form of averaging is not a true average if more that one value is combined.
634  // This might be an issue in the future, but we leave this way for now.
635  else if(CombineMode==Average)
636  for (j=0; j<NumImportIDs; j++) {To[ImportLIDs[j]] += *ptr++; To[ImportLIDs[j]] /= 2;}
637  }
638 
639  // constant element size case
640 
641  else if (ConstantElementSize) {
642 
643  if (CombineMode==Add) {
644  for (j=0; j<NumImportIDs; j++) {
645  jj = MaxElementSize*ImportLIDs[j];
646  for (k=0; k<MaxElementSize; k++)
647  To[jj+k] += *ptr++; // Add to existing value
648  }
649  }
650  else if(CombineMode==Insert) {
651  for (j=0; j<NumImportIDs; j++) {
652  jj = MaxElementSize*ImportLIDs[j];
653  for (k=0; k<MaxElementSize; k++)
654  To[jj+k] = *ptr++;
655  }
656  }
657  else if(CombineMode==AbsMax) {
658  for (j=0; j<NumImportIDs; j++) {
659  jj = MaxElementSize*ImportLIDs[j];
660  for (k=0; k<MaxElementSize; k++) {
661  To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
662  ptr++;
663  }
664  }
665  }
666  // Note: The following form of averaging is not a true average if more that one value is combined.
667  // This might be an issue in the future, but we leave this way for now.
668  else if(CombineMode==Average) {
669  for (j=0; j<NumImportIDs; j++) {
670  jj = MaxElementSize*ImportLIDs[j];
671  for (k=0; k<MaxElementSize; k++)
672  { To[jj+k] += *ptr++; To[jj+k] /= 2;}
673  }
674  }
675  }
676 
677  // variable element size case
678 
679  else {
680 
681  SizeOfPacket = MaxElementSize;
682 
683  if (CombineMode==Add) {
684  for (j=0; j<NumImportIDs; j++) {
685  ptr = (double *) Imports + j*SizeOfPacket;
686  jj = ToFirstPointInElementList[ImportLIDs[j]];
687  int ElementSize = ToElementSizeList[ImportLIDs[j]];
688  for (k=0; k<ElementSize; k++)
689  To[jj+k] += *ptr++; // Add to existing value
690  }
691  }
692  else if(CombineMode==Insert){
693  for (j=0; j<NumImportIDs; j++) {
694  ptr = (double *) Imports + j*SizeOfPacket;
695  jj = ToFirstPointInElementList[ImportLIDs[j]];
696  int ElementSize = ToElementSizeList[ImportLIDs[j]];
697  for (k=0; k<ElementSize; k++)
698  To[jj+k] = *ptr++;
699  }
700  }
701  else if(CombineMode==AbsMax){
702  for (j=0; j<NumImportIDs; j++) {
703  ptr = (double *) Imports + j*SizeOfPacket;
704  jj = ToFirstPointInElementList[ImportLIDs[j]];
705  int ElementSize = ToElementSizeList[ImportLIDs[j]];
706  for (k=0; k<ElementSize; k++) {
707  To[jj+k] = EPETRA_MAX( To[jj+k], std::abs(*ptr));
708  ptr++;
709  }
710  }
711  }
712  // Note: The following form of averaging is not a true average if more that one value is combined.
713  // This might be an issue in the future, but we leave this way for now.
714  else if(CombineMode==Average) {
715  for (j=0; j<NumImportIDs; j++) {
716  ptr = (double *) Imports + j*SizeOfPacket;
717  jj = ToFirstPointInElementList[ImportLIDs[j]];
718  int ElementSize = ToElementSizeList[ImportLIDs[j]];
719  for (k=0; k<ElementSize; k++)
720  { To[jj+k] += *ptr++; To[jj+k] /= 2;}
721  }
722  }
723  }
724 
725  return(0);
726 }
727 
int NumMyUnknowns() const
Returns the number of local unknowns.
int BlockSize(int LID) const
Returns the size of the given block.
int NumMyBlocks() const
Returns the number of local blocks.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual const Epetra_BlockMap & DataMap() const
Returns the Epetra_BlockMap object with the distribution of underlying values.
EpetraExt_BlockDiagMatrix & operator=(const EpetraExt_BlockDiagMatrix &Source)
= Operator.
virtual ~EpetraExt_BlockDiagMatrix()
Destructor.
EpetraExt_BlockDiagMatrix: A class for storing distributed block matrices.
EpetraExt_BlockDiagMatrix(const Epetra_BlockMap &Map, bool zero_out=true)
Constructor - This map is the map of the vector this can be applied to.
#define AM_MULTIPLY
#define AM_FACTOR
virtual int Compute()
Computes the inverse / factorization if such is set on the list.
virtual void Print(std::ostream &os) const
Print method.
virtual int SetParameters(Teuchos::ParameterList &List)
SetParameters.
#define AM_INVERT
double * Values() const
Returns a pointer to the array containing the blocks.
int NumData() const
Returns the size of the total Data block.
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...
void PutScalar(double value)
PutScalar function.