EpetraExt Package Browser (Single Doxygen Collection) 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
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
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
154//=========================================================================
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
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
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
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);
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
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
325 MyGlobalElements1_LL = DataMap_->MyGlobalElements64();
326 }
327 else
328#endif
329 throw "EpetraExt_BlockDiagMatrix::Print: GlobalIndices type unknown";
330
331 int * FirstPointInElementList1 = (MaxElementSize1 == 1) ?
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.
380 return &Map() == &Source.Map();
381}
382
383
384 //=========================================================================
385// Perform ID copies and permutations that are on processor.
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().
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().
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
Insert
Add
AbsMax
Average
Zero
Epetra_AddLocalAlso
#define EPETRA_MAX(x, y)
#define EPETRA_CHK_ERR(a)
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 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)
int NumMyUnknowns() const
Returns the number of local unknowns.
int * Pivots_
Pivots for factorization.
bool HasComputed_
Has Computed? Needed for Inverse/Factorization modes.
int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
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.
int PackAndPrepare(const Epetra_SrcDistObject &Source, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
virtual void Print(std::ostream &os) const
Print method.
EpetraExt_BlockDiagMatrix & operator=(const EpetraExt_BlockDiagMatrix &Source)
= Operator.
Epetra_BlockMap * DataMap_
Map for the data.
int CheckSizes(const Epetra_SrcDistObject &Source)
int DoCopy(const EpetraExt_BlockDiagMatrix &Source)
virtual ~EpetraExt_BlockDiagMatrix()
Destructor.
int NumData() const
Returns the size of the total Data block.
double * Values_
Actual Data values.
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 ApplyMode_
Which Apply Mode to use.
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
long long * MyGlobalElements64() 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