45 #include "Epetra_MultiVector.h" 46 #include "Epetra_Comm.h" 47 #include "Epetra_LAPACK.h" 48 #include "Epetra_Distributor.h" 57 : Epetra_DistObject(map,
"EpetraExt::BlockDiagMatrix"),
72 if(DataMap_)
delete DataMap_;
73 if(Pivots_)
delete [] Pivots_;
74 if(Values_)
delete [] Values_;
81 : Epetra_DistObject(Source),
87 DataMap_=
new Epetra_BlockMap(*Source.DataMap_);
89 Values_=
new double[DataMap_->NumMyPoints()];
95 void 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()) {
108 DataMap_=
new Epetra_BlockMap(-1,Map().NumMyElements(),Map().MyGlobalElements(),ElementSize,0,Map().
Comm());
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());
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++){
188 Values_[DataMap_->FirstPointInElement(i)]=1.0/Values_[DataMap_->FirstPointInElement(i)];
192 double * v=&Values_[DataMap_->FirstPointInElement(i)];
193 double d=1/(v[0]*v[3]-v[1]*v[2]);
202 LAPACK.GETRF(Nb,Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&info);
203 if(info) EPETRA_CHK_ERR(-2);
209 int lwork=3*DataMap_->MaxMyElementSize();
210 std::vector<double> work(lwork);
211 for(
int i=0;i<NumBlocks;i++){
219 LAPACK.GETRI(Nb,&Values_[DataMap_->FirstPointInElement(i)],Nb,&Pivots_[Map().FirstPointInElement(i)],&work[0],&lwork,&info);
220 if(info) EPETRA_CHK_ERR(-3);
234 int NumVectors=X.NumVectors();
235 if(NumVectors!=Y.NumVectors())
243 const int *vlist=DataMap_->FirstPointInElementList();
244 const int *xlist=Map().FirstPointInElementList();
245 const int *blocksize=Map().ElementSizeList();
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);
308 int MyPID = DataMap_->Comm().MyPID();
309 int NumProc = DataMap_->Comm().NumProc();
311 for (
int iproc=0; iproc < NumProc; 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();
323 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 324 if (DataMap_->GlobalIndicesLongLong ()) {
325 MyGlobalElements1_LL = DataMap_->MyGlobalElements64();
329 throw "EpetraExt_BlockDiagMatrix::Print: GlobalIndices type unknown";
331 int * FirstPointInElementList1 = (MaxElementSize1 == 1) ?
332 NULL : DataMap_->FirstPointInElementList();
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;
369 Map().Comm().Barrier();
370 Map().Comm().Barrier();
371 Map().Comm().Barrier();
379 int EpetraExt_BlockDiagMatrix::CheckSizes(
const Epetra_SrcDistObject& Source){
380 return &Map() == &Source.Map();
386 int EpetraExt_BlockDiagMatrix::CopyAndPermute(
const Epetra_SrcDistObject& Source,
390 int * PermuteFromLIDs,
391 const Epetra_OffsetIndex * Indexor,
392 Epetra_CombineMode CombineMode){
398 double *To = Values_;
400 int * ToFirstPointInElementList = 0;
401 int * FromFirstPointInElementList = 0;
402 int * FromElementSizeList = 0;
403 int MaxElementSize =
DataMap().MaxElementSize();
404 bool ConstantElementSize =
DataMap().ConstantElementSize();
406 if (!ConstantElementSize) {
407 ToFirstPointInElementList =
DataMap().FirstPointInElementList();
408 FromFirstPointInElementList = A.
DataMap().FirstPointInElementList();
409 FromElementSizeList = A.
DataMap().ElementSizeList();
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];
508 int EpetraExt_BlockDiagMatrix::PackAndPrepare(
const Epetra_SrcDistObject& Source,
516 Epetra_Distributor& Distor){
525 int MaxElementSize =
DataMap().MaxElementSize();
526 bool ConstantElementSize =
DataMap().ConstantElementSize();
528 int * FromFirstPointInElementList = 0;
529 int * FromElementSizeList = 0;
531 if (!ConstantElementSize) {
532 FromFirstPointInElementList = A.
DataMap().FirstPointInElementList();
533 FromElementSizeList = A.
DataMap().ElementSizeList();
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++)
580 int EpetraExt_BlockDiagMatrix::UnpackAndCombine(
const Epetra_SrcDistObject& Source,
586 Epetra_Distributor& Distor,
587 Epetra_CombineMode CombineMode,
588 const Epetra_OffsetIndex * Indexor){
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_;
605 int MaxElementSize =
DataMap().MaxElementSize();
606 bool ConstantElementSize =
DataMap().ConstantElementSize();
608 int * ToFirstPointInElementList = 0;
609 int * ToElementSizeList = 0;
611 if (!ConstantElementSize) {
612 ToFirstPointInElementList =
DataMap().FirstPointInElementList();
613 ToElementSizeList =
DataMap().ElementSizeList();
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;}
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.
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.
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.