45#include "Epetra_MultiVector.h"
46#include "Epetra_Comm.h"
47#include "Epetra_LAPACK.h"
48#include "Epetra_Distributor.h"
72 if(DataMap_)
delete DataMap_;
73 if(Pivots_)
delete [] Pivots_;
74 if(Values_)
delete [] Values_;
95void EpetraExt_BlockDiagMatrix::Allocate(){
99 int *ElementSize=
new int[NumBlocks];
101 for(
int i=0;i<NumBlocks;i++) {
103 dataSize+=ElementSize[i];
106#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
107 if(
Map().GlobalIndicesInt()) {
112#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
113 if(
Map().GlobalIndicesLongLong()) {
118 throw "EpetraExt_BlockDiagMatrix::Allocate: GlobalIndices type unknown";
120 Values_=
new double[dataSize];
121 delete [] ElementSize;
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);
145 for (
int i=0;i<MaxData;i++) Values_[i]=value;
157 if(!
Map().SameAs(Source.
Map()) || !DataMap_->
SameAs(*Source.DataMap_))
158 throw ReportError(
"Maps of DistBlockMatrices incompatible in assignment",-1);
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];
166 ApplyMode_=Source.ApplyMode_;
167 HasComputed_=Source.HasComputed_;
184 for(
int i=0;i<NumBlocks;i++){
193 double d=1/(v[0]*v[3]-v[1]*v[2]);
203 if(info) EPETRA_CHK_ERR(-2);
210 std::vector<double> work(lwork);
211 for(
int i=0;i<NumBlocks;i++){
220 if(info) EPETRA_CHK_ERR(-3);
250 for(
int i=0;i<NumBlocks;i++){
254 for(
int j=0;j<NumVectors;j++){
257 Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
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];
267 GEMV(
'N',Nb,Nb,1.0,&Values_[vidx0],Nb,&X[j][xidx0],0.0,&Y[j][xidx0]);
275 for(
int i=0;i<NumBlocks;i++){
279 for(
int j=0;j<NumVectors;j++){
282 Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
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];
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);
311 for (
int iproc=0; iproc < NumProc; iproc++) {
315 int * MyGlobalElements1_int = 0;
316 long long * MyGlobalElements1_LL = 0;
317#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
323#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
325 MyGlobalElements1_LL = DataMap_->MyGlobalElements64();
329 throw "EpetraExt_BlockDiagMatrix::Print: GlobalIndices type unknown";
331 int * FirstPointInElementList1 = (MaxElementSize1 == 1) ?
336 os <<
" MyPID"; os <<
" ";
338 if (MaxElementSize1==1)
346 for (
int i=0; i < NumMyElements1; i++) {
347 for (
int ii=0; ii < DataMap_->
ElementSize(i); ii++) {
350 os << MyPID; os <<
" ";
352 if (MaxElementSize1==1) {
353 os << (MyGlobalElements1_int ? MyGlobalElements1_int[i] : MyGlobalElements1_LL[i]) <<
" ";
357 os << (MyGlobalElements1_int ? MyGlobalElements1_int[i] : MyGlobalElements1_LL[i]) <<
"/" << ii <<
" ";
358 iii = FirstPointInElementList1[i]+ii;
380 return &
Map() == &Source.
Map();
390 int * PermuteFromLIDs,
398 double *To = Values_;
400 int * ToFirstPointInElementList = 0;
401 int * FromFirstPointInElementList = 0;
402 int * FromElementSizeList = 0;
406 if (!ConstantElementSize) {
418 if (MaxElementSize==1) {
420 NumSameEntries = NumSameIDs;
422 else if (ConstantElementSize) {
424 NumSameEntries = NumSameIDs * MaxElementSize;
428 NumSameEntries = FromFirstPointInElementList[NumSameIDs];
439 if (CombineMode==Epetra_AddLocalAlso) {
440 for (
int j=0; j<NumSameEntries; j++) {
444 for (
int j=0; j<NumSameEntries; j++) {
450 if (NumPermuteIDs>0) {
455 if (CombineMode==Epetra_AddLocalAlso) {
456 for (
int j=0; j<NumPermuteIDs; j++) {
457 To[PermuteToLIDs[j]] += From[PermuteFromLIDs[j]];
461 for (
int j=0; j<NumPermuteIDs; j++) {
462 To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
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];
477 for (
int k=0; k<MaxElementSize; k++) {
478 To[jj+k] = From[jjj+k];
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];
496 for (
int k=0; k<ElementSize; k++) {
497 To[jj+k] = From[jjj+k];
528 int * FromFirstPointInElementList = 0;
529 int * FromElementSizeList = 0;
531 if (!ConstantElementSize) {
536 SizeOfPacket = MaxElementSize * (int)
sizeof(
double);
538 if(NumExportIDs*SizeOfPacket>LenExports) {
539 if (LenExports>0)
delete [] Exports;
540 LenExports = NumExportIDs*SizeOfPacket;
541 Exports =
new char[LenExports];
546 if (NumExportIDs>0) {
547 ptr = (
double *) Exports;
550 if (MaxElementSize==1)
for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];
553 else if (ConstantElementSize) {
554 for (j=0; j<NumExportIDs; j++) {
555 jj = MaxElementSize*ExportLIDs[j];
556 for (k=0; k<MaxElementSize; k++)
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++)
595 if( CombineMode != Add
596 && CombineMode != Zero
597 && CombineMode != Insert
598 && CombineMode != Average
599 && CombineMode != AbsMax )
602 if (NumImportIDs<=0)
return(0);
604 double * To = Values_;
608 int * ToFirstPointInElementList = 0;
609 int * ToElementSizeList = 0;
611 if (!ConstantElementSize) {
619 ptr = (
double *) Imports;
622 if (MaxElementSize==1) {
624 if (CombineMode==Add)
625 for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += *ptr++;
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));
635 else if(CombineMode==Average)
636 for (j=0; j<NumImportIDs; j++) {To[ImportLIDs[j]] += *ptr++; To[ImportLIDs[j]] /= 2;}
641 else if (ConstantElementSize) {
643 if (CombineMode==Add) {
644 for (j=0; j<NumImportIDs; j++) {
645 jj = MaxElementSize*ImportLIDs[j];
646 for (k=0; k<MaxElementSize; k++)
650 else if(CombineMode==Insert) {
651 for (j=0; j<NumImportIDs; j++) {
652 jj = MaxElementSize*ImportLIDs[j];
653 for (k=0; k<MaxElementSize; k++)
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));
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;}
681 SizeOfPacket = MaxElementSize;
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++)
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++)
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));
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;}
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 * ElementSizeList() const
int * FirstPointInElementList() const
int MaxMyElementSize() const
bool SameAs(const Epetra_BlockMap &Map) const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
bool ConstantElementSize() const
int NumMyElements() const
int FirstPointInElement(int LID) const
bool GlobalIndicesLongLong() const
virtual int NumProc() const=0
virtual int MyPID() const=0
virtual void Barrier() const=0
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
virtual int ReportError(const std::string Message, int ErrorCode) const
virtual const Epetra_BlockMap & Map() const=0