EpetraExt Development
Loading...
Searching...
No Matches
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
95void 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//=========================================================================
127int 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//=========================================================================
155int 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
307void 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.
379int 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.
386int 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().
508int 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().
580int 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
#define AM_MULTIPLY
#define AM_FACTOR
#define AM_INVERT
Epetra_CombineMode
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.
int NumMyUnknowns() const
Returns the number of local unknowns.
virtual int Compute()
Computes the inverse / factorization if such is set on the list.
int BlockSize(int LID) const
Returns the size of the given block.
virtual void Print(std::ostream &os) const
Print method.
EpetraExt_BlockDiagMatrix & operator=(const EpetraExt_BlockDiagMatrix &Source)
= Operator.
virtual ~EpetraExt_BlockDiagMatrix()
Destructor.
int NumData() const
Returns the size of the total Data block.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
void PutScalar(double value)
PutScalar function.
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.
EpetraExt_BlockDiagMatrix(const Epetra_BlockMap &Map, bool zero_out=true)
Constructor - This map is the map of the vector this can be applied to.
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.
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
int MaxElementSize() const
int MyGlobalElements(int *MyGlobalElementList) const
int NumMyPoints() const
int ElementSize() const
int * FirstPointInElementList() const
int MaxMyElementSize() const
bool SameAs(const Epetra_BlockMap &Map) const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
int NumMyElements() const
int FirstPointInElement(int LID) const
bool GlobalIndicesLongLong() const
const Epetra_BlockMap & Map() const
void GETRF(const int M, const int N, float *A, const int LDA, int *IPIV, int *INFO) const
void GETRI(const int N, float *A, const int LDA, int *IPIV, float *WORK, const int *LWORK, int *INFO) 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
int NumVectors() const
virtual int ReportError(const std::string Message, int ErrorCode) const