EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
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 //=========================================================================
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 
virtual const Epetra_BlockMap & DataMap() const
Returns the Epetra_BlockMap object with the distribution of underlying values.
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 ElementSize() const
bool GlobalIndicesLongLong() const
bool SameAs(const Epetra_BlockMap &Map) const
int NumMyBlocks() const
Returns the number of local blocks.
int MyGlobalElements(int *MyGlobalElementList) const
void GEMV(const char TRANS, const int M, const int N, const float ALPHA, const float *A, const int LDA, const float *X, const float BETA, float *Y, const int INCX=1, const int INCY=1) const
double * Values() const
Returns a pointer to the array containing the blocks.
bool ConstantElementSize() const
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
bool GlobalIndicesInt() const
int NumMyUnknowns() const
Returns the number of local unknowns.
EpetraExt_BlockDiagMatrix & operator=(const EpetraExt_BlockDiagMatrix &Source)
= Operator.
int * FirstPointInElementList() const
virtual int MyPID() const =0
void GETRF(const int M, const int N, float *A, const int LDA, int *IPIV, int *INFO) const
virtual ~EpetraExt_BlockDiagMatrix()
Destructor.
EpetraExt_BlockDiagMatrix: A class for storing distributed block matrices.
int NumMyElements() const
int BlockSize(int LID) const
Returns the size of the given block.
virtual void Print(std::ostream &os) const
Print method.
int NumData() const
Returns the size of the total Data block.
int FirstPointInElement(int LID) const
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.
int * ElementSizeList() const
void GETRS(const char TRANS, const int N, const int NRHS, const float *A, const int LDA, const int *IPIV, float *X, const int LDX, int *INFO) const
void GETRI(const int N, float *A, const int LDA, int *IPIV, float *WORK, const int *LWORK, int *INFO) const
virtual int SetParameters(Teuchos::ParameterList &List)
SetParameters.
const Epetra_Comm & Comm() const
virtual int NumProc() const =0
#define AM_INVERT
int MaxElementSize() const
Epetra_CombineMode
int MaxMyElementSize() const
int NumMyPoints() const
void PutScalar(double value)
PutScalar function.