45 #include "Epetra_RowMatrix.h" 46 #include "Epetra_Map.h" 47 #include "Epetra_CrsMatrix.h" 48 #include "Epetra_Comm.h" 49 #include "Epetra_MultiVector.h" 51 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 52 #include "Epetra_IntVector.h" 53 #include "Epetra_MpiComm.h" 54 #include "Teuchos_Hashtable.hpp" 55 #include "Teuchos_Array.hpp" 56 #include "EpetraExt_OperatorOut.h" 57 #include "EpetraExt_BlockMapOut.h" 59 # ifdef IFPACK_NODE_AWARE_CODE 60 # include "Epetra_IntVector.h" 61 # include "Epetra_MpiComm.h" 62 # include "Teuchos_Hashtable.hpp" 63 # include "Teuchos_Array.hpp" 64 # include "EpetraExt_OperatorOut.h" 65 # include "EpetraExt_BlockMapOut.h" 66 extern int ML_NODE_ID;
75 template<
typename int_type>
84 const Epetra_Map *RowMap;
85 const Epetra_Map *ColMap;
87 std::vector<int_type> ExtElements;
89 for (
int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap) {
90 if (TmpMatrix != Teuchos::null) {
91 RowMap = &(TmpMatrix->RowMatrixRowMap());
92 ColMap = &(TmpMatrix->RowMatrixColMap());
95 RowMap = &(
A().RowMatrixRowMap());
96 ColMap = &(
A().RowMatrixColMap());
99 int size = ColMap->NumMyElements() - RowMap->NumMyElements();
100 std::vector<int_type> list(size);
105 for (
int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
106 int_type GID = (int_type) ColMap->GID64(i);
107 if (
A().RowMatrixRowMap().LID(GID) == -1) {
108 typename std::vector<int_type>::iterator pos
109 =
find(ExtElements.begin(),ExtElements.end(),GID);
110 if (pos == ExtElements.end()) {
111 ExtElements.push_back(GID);
118 const int_type *listptr = NULL;
119 if ( ! list.empty() ) listptr = &list[0];
120 TmpMap =
rcp(
new Epetra_Map(-1,count, listptr,0,
Comm()) );
122 TmpMatrix =
rcp(
new Epetra_CrsMatrix(
Copy,*TmpMap,0) );
124 TmpImporter =
rcp(
new Epetra_Import(*TmpMap,
A().RowMatrixRowMap()) );
126 TmpMatrix->Import(
A(),*TmpImporter,Insert);
127 TmpMatrix->FillComplete(
A().OperatorDomainMap(),*TmpMap);
133 std::vector<int_type> list(NumMyRowsA_ + ExtElements.size());
134 for (
int i = 0 ; i < NumMyRowsA_ ; ++i)
135 list[i] = (int_type)
A().RowMatrixRowMap().GID64(i);
136 for (
int i = 0 ; i < (int)ExtElements.size() ; ++i)
137 list[i + NumMyRowsA_] = ExtElements[i];
139 const int_type *listptr = NULL;
140 if ( ! list.empty() ) listptr = &list[0];
142 Map_ =
rcp(
new Epetra_Map((int_type) -1, NumMyRowsA_ + ExtElements.size(),
143 listptr, 0,
Comm()) );
145 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 148 # ifdef IFPACK_NODE_AWARE_CODE 155 const int_type * extelsptr = NULL;
156 if ( ! ExtElements.empty() ) extelsptr = &ExtElements[0];
157 ExtMap_ =
rcp(
new Epetra_Map((int_type) -1,ExtElements.size(),
158 extelsptr,0,
A().Comm()) );
166 int OverlapLevel_in) :
168 OverlapLevel_(OverlapLevel_in)
171 if (OverlapLevel_in == 0)
175 if (
Comm().NumProc() == 1)
182 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 184 BuildMap<int>(OverlapLevel_in);
188 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 189 if(
A().RowMatrixRowMap().GlobalIndicesLongLong()) {
190 BuildMap<long long>(OverlapLevel_in);
194 throw "Ifpack_OverlappingRowMatrix::Ifpack_OverlappingRowMatrix: GlobalIndices type unknown";
221 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 227 int OverlapLevel_in,
int subdomainID) :
229 OverlapLevel_(OverlapLevel_in)
233 #ifdef IFPACK_OVA_TIME_BUILD 234 double t0 = MPI_Wtime();
235 double t1, true_t0=t0;
239 const int votesMax = INT_MAX;
242 if (OverlapLevel_in == 0)
246 if (
Comm().NumProc() == 1)
253 const Epetra_MpiComm *pComm =
dynamic_cast<const Epetra_MpiComm*
>( &
Comm() );
254 assert(pComm != NULL);
255 MPI_Comm nodeMPIComm;
256 MPI_Comm_split(pComm->Comm(),subdomainID,pComm->MyPID(),&nodeMPIComm);
257 Epetra_MpiComm *nodeComm =
new Epetra_MpiComm(nodeMPIComm);
259 Epetra_SerialComm *nodeComm =
dynamic_cast<Epetra_MpiComm*
>( &(Matrix_in->RowMatrixRowMap().Comm()) );
270 const Epetra_Map *RowMap;
271 const Epetra_Map *ColMap;
272 const Epetra_Map *DomainMap;
278 #ifdef IFPACK_OVA_TIME_BUILD 280 fprintf(stderr,
"[node %d]: time for initialization %2.3e\n", subdomainID, t1-t0);
285 #ifdef IFPACK_OVA_TIME_BUILD 287 fprintf(stderr,
"[node %d]: overlap hash table allocs %2.3e\n", subdomainID, t1-t0);
300 for (
int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
305 if (TmpMatrix != Teuchos::null) {
306 RowMap = &(TmpMatrix->RowMatrixRowMap());
307 ColMap = &(TmpMatrix->RowMatrixColMap());
308 DomainMap = &(TmpMatrix->OperatorDomainMap());
311 RowMap = &(
A().RowMatrixRowMap());
312 ColMap = &(
A().RowMatrixColMap());
313 DomainMap = &(
A().OperatorDomainMap());
316 #ifdef IFPACK_OVA_TIME_BUILD 318 fprintf(stderr,
"[node %d]: overlap pointer pulls %2.3e\n", subdomainID, t1-t0);
324 Epetra_IntVector colIdList( *ColMap );
325 Epetra_IntVector rowIdList(*DomainMap);
326 rowIdList.PutValue(subdomainID);
329 #ifdef IFPACK_OVA_TIME_BUILD 331 fprintf(stderr,
"[node %d]: overlap intvector/importer allocs %2.3e\n", subdomainID, t1-t0);
335 colIdList.Import(rowIdList,*nodeIdImporter,Insert);
337 int size = ColMap->NumMyElements() - RowMap->NumMyElements();
340 #ifdef IFPACK_OVA_TIME_BUILD 342 fprintf(stderr,
"[node %d]: overlap col/row ID imports %2.3e\n", subdomainID, t1-t0);
349 for (
int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
350 int GID = ColMap->GID(i);
352 if (ghostTable.containsKey(GID)) myvotes = ghostTable.get(GID);
353 if ( colIdList[i] != subdomainID && myvotes < votesMax)
356 if (ghostTable.containsKey(GID)) {
357 votes = ghostTable.get(GID);
359 ghostTable.put(GID,votes);
361 ghostTable.put(GID,1);
367 ghostTable.arrayify(gidsarray,votesarray);
380 # ifdef HAVE_MPI //FIXME What if we build in serial?! This file will likely not compile. 386 #ifdef IFPACK_OVA_TIME_BUILD 388 fprintf(stderr,
"[node %d]: overlap before culling %2.3e\n", subdomainID, t1-t0);
394 if (nodeComm->MyPID() == 0)
398 int *lengths =
new int[nodeComm->NumProc()+1];
400 lengths[1] = ghostTable.size();
401 for (
int i=1; i<nodeComm->NumProc(); i++) {
403 MPI_Recv( &leng, 1, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
404 lengths[i+1] = lengths[i] + leng;
406 int total = lengths[nodeComm->NumProc()];
408 int* ghosts =
new int[total];
409 for (
int i=0; i<total; i++) ghosts[i] = -9;
410 int *round =
new int[total];
411 int *owningpid =
new int[total];
413 for (
int i=1; i<nodeComm->NumProc(); i++) {
414 int count = lengths[i+1] - lengths[i];
415 MPI_Recv( ghosts+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
416 MPI_Recv( round+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
420 for (
int i=0; i<lengths[1]; i++) {
427 for (
int j=1; j<nodeComm->NumProc(); j++) {
428 for (
int i=lengths[j]; i<lengths[j+1]; i++) {
435 roundpid[0] = round; roundpid[1] = owningpid;
436 Epetra_Util epetraUtil;
437 epetraUtil.Sort(
true,total,ghosts,0,0,2,roundpid);
440 int *nlosers =
new int[nodeComm->NumProc()];
441 int** losers =
new int*[nodeComm->NumProc()];
442 for (
int i=0; i<nodeComm->NumProc(); i++) {
444 losers[i] =
new int[ lengths[i+1]-lengths[i] ];
456 for (
int i=1; i<total; i++) {
457 if (ghosts[i] == ghosts[i-1]) {
459 if (round[i] > round[
max]) {
463 int j = owningpid[idx];
464 losers[j][nlosers[j]++] = ghosts[idx];
475 for (
int i=1; i<nodeComm->NumProc(); i++) {
476 MPI_Send( nlosers+i, 1, MPI_INT, i, 8675, nodeComm->Comm());
477 MPI_Send( losers[i], nlosers[i], MPI_INT, i, 8675, nodeComm->Comm());
484 cullList =
new int[ncull+1];
485 for (
int i=0; i<nlosers[0]; i++)
486 cullList[i] = losers[0][i];
488 for (
int i=0; i<nodeComm->NumProc(); i++)
498 int hashsize = ghostTable.size();
499 MPI_Send( &hashsize, 1, MPI_INT, 0, 8675, nodeComm->Comm());
500 MPI_Send( gids, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
501 MPI_Send( votes, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
505 MPI_Recv( &ncull, 1, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
506 cullList =
new int[ncull+1];
507 MPI_Recv( cullList, ncull, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
514 for (
int i=0; i<ncull; i++) {
515 try{ghostTable.remove(cullList[i]);}
518 printf(
"pid %d: In OverlappingRowMatrix ctr, problem removing ghost elt %d from ghostTable\n",
519 Comm().MyPID(),cullList[i]);
529 ghostTable.arrayify(gidsarray,votesarray);
531 std::vector<int> list(size);
533 for (
int i=0; i<ghostTable.size(); i++) {
535 if (votesarray[i] < votesMax) {
536 list[count] = gidsarray[i];
537 ghostTable.put(gidsarray[i],votesMax);
543 #ifdef IFPACK_OVA_TIME_BUILD 545 fprintf(stderr,
"[node %d]: overlap duplicate removal %2.3e\n", subdomainID, t1-t0);
550 # endif //ifdef HAVE_MPI 552 TmpMap =
rcp(
new Epetra_Map(-1,count, &list[0],0,
Comm()) );
554 TmpMatrix =
rcp(
new Epetra_CrsMatrix(Copy,*TmpMap,0) );
558 TmpMatrix->Import(
A(),*TmpImporter,Insert);
562 #ifdef IFPACK_OVA_TIME_BUILD 564 fprintf(stderr,
"[node %d]: overlap TmpMatrix fillcomplete %2.3e\n", subdomainID, t1-t0);
575 Epetra_IntVector ov_colIdList( TmpMatrix->ColMap() );
576 ov_colIdList.PutValue(-1);
577 Epetra_IntVector ov_rowIdList( TmpMatrix->RowMap() );
578 ov_rowIdList.PutValue(subdomainID);
580 ov_colIdList.Import(ov_rowIdList,*ov_nodeIdImporter,Insert);
582 for (
int i=0 ; i<ov_colIdList.MyLength(); i++) {
583 if (ov_colIdList[i] == subdomainID) {
584 int GID = (ov_colIdList.Map()).GID(i);
585 colMapTable.put(GID,1);
592 ov_colIdList.PutValue(-1);
594 ArowIdList.PutValue(subdomainID);
595 nodeIdImporter =
rcp(
new Epetra_Import( TmpMatrix->ColMap(),
A().RowMatrixRowMap() ));
596 ov_colIdList.Import(ArowIdList,*nodeIdImporter,Insert);
598 for (
int i=0 ; i<ov_colIdList.MyLength(); i++) {
599 if (ov_colIdList[i] == subdomainID) {
600 int GID = (ov_colIdList.Map()).GID(i);
601 colMapTable.put(GID,1);
606 #ifdef IFPACK_OVA_TIME_BUILD 608 fprintf(stderr,
"[node %d]: overlap 2 imports to fix up colmap %2.3e\n", subdomainID, t1-t0);
620 std::vector<int> ghostElements;
623 ghostTable.arrayify(gidsarray,votesarray);
624 for (
int i=0; i<ghostTable.size(); i++) {
628 for (
int i = 0 ; i <
A().RowMatrixColMap().NumMyElements() ; ++i) {
629 int GID =
A().RowMatrixColMap().GID(i);
631 if (colMapTable.containsKey(GID)) {
632 try{colMapTable.remove(GID);}
634 printf(
"pid %d: In OverlappingRowMatrix ctr, problem removing entry %d from colMapTable\n",
Comm().MyPID(),GID);
641 std::vector<int> colMapElements;
644 colMapTable.arrayify(gidsarray,votesarray);
645 for (
int i=0; i<colMapTable.size(); i++)
646 colMapElements.push_back(gidsarray[i]);
659 std::vector<int> rowList(
NumMyRowsA_ + ghostElements.size());
661 rowList[i] =
A().RowMatrixRowMap().GID(i);
662 for (
int i = 0 ; i < (int)ghostElements.size() ; ++i)
674 std::vector<int> colList(
A().
RowMatrixColMap().NumMyElements() + colMapElements.size());
675 int nc =
A().RowMatrixColMap().NumMyElements();
676 for (
int i = 0 ; i < nc; i++)
677 colList[i] =
A().RowMatrixColMap().GID(i);
678 for (
int i = 0 ; i < (int)colMapElements.size() ; i++)
679 colList[nc+i] = colMapElements[i];
683 colMap_ =
new Epetra_Map(-1,
A().
RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0,
Comm());
687 #ifdef IFPACK_OVA_TIME_BUILD 689 fprintf(stderr,
"[node %d]: time to init B() col/row maps %2.3e\n", subdomainID, t1-t0);
701 Epetra_Map* nodeMap =
new Epetra_Map(-1,
NumMyRowsA_ + ghostElements.size(),&rowList[0],0,*nodeComm);
702 int numMyElts = colMap_->NumMyElements();
703 assert(numMyElts!=0);
709 int* myGlobalElts =
new int[numMyElts];
710 colMap_->MyGlobalElements(myGlobalElts);
711 int* pidList =
new int[numMyElts];
712 nodeMap->RemoteIDList(numMyElts, myGlobalElts, pidList, 0);
717 tt[0] = myGlobalElts;
718 Util.Sort(
true, numMyElts, pidList, 0, (
double**)0, 1, tt);
723 while (localStart<numMyElts) {
724 int currPID = (pidList)[localStart];
726 while (i<numMyElts && (pidList)[i] == currPID) i++;
727 if (currPID != nodeComm->MyPID())
728 Util.Sort(
true, i-localStart, (myGlobalElts)+localStart, 0, 0, 0, 0);
734 while (localStart<numMyElts && (pidList)[localStart] != nodeComm->MyPID()) localStart++;
735 assert(localStart != numMyElts);
736 int localEnd=localStart;
737 while (localEnd<numMyElts && (pidList)[localEnd] == nodeComm->MyPID()) localEnd++;
738 int* mySortedGlobalElts =
new int[numMyElts];
739 int leng = localEnd - localStart;
743 int *rowGlobalElts = nodeMap->MyGlobalElements();
763 for (
int i=0; i<leng; i++) {
765 (mySortedGlobalElts)[i] = rowGlobalElts[i];
769 for (
int i=leng; i< localEnd; i++) {
770 (mySortedGlobalElts)[i] = (myGlobalElts)[i-leng];
772 for (
int i=localEnd; i<numMyElts; i++) {
773 (mySortedGlobalElts)[i] = (myGlobalElts)[i];
777 #ifdef IFPACK_OVA_TIME_BUILD 779 fprintf(stderr,
"[node %d]: time to sort col map arrays (cp=1) %2.3e\n", subdomainID, t1-t0);
784 int indexBase = colMap_->IndexBase();
785 if (colMap_)
delete colMap_;
786 colMap_ =
new Epetra_Map(-1,numMyElts,mySortedGlobalElts,indexBase,
Comm());
789 delete [] myGlobalElts;
791 delete [] mySortedGlobalElts;
795 printf(
"** * Ifpack_OverlappingRowmatrix ctor: problem creating column map * **\n\n");
803 #ifdef IFPACK_OVA_TIME_BUILD 805 fprintf(stderr,
"[node %d]: time to sort col map arrays (cp=2) %2.3e\n", subdomainID, t1-t0);
835 ExtMap_ =
rcp(
new Epetra_Map(-1,ghostElements.size(), &ghostElements[0],0,
Comm()) );
842 #ifdef IFPACK_OVA_TIME_BUILD 844 fprintf(stderr,
"[node %d]: time to import and setup ExtMap_ %2.3e\n", subdomainID, t1-t0);
860 #ifdef IFPACK_OVA_TIME_BUILD 862 fprintf(stderr,
"[node %d]: time to FillComplete B() %2.3e\n", subdomainID, t1-t0);
896 #ifdef IFPACK_OVA_TIME_BUILD 898 fprintf(stderr,
"[node %d]: time for final calcs %2.3e\n", subdomainID, t1-t0);
899 fprintf(stderr,
"[node %d]: total IORM ctor time %2.3e\n", subdomainID, t1-true_t0);
907 # ifdef IFPACK_NODE_AWARE_CODE 913 int OverlapLevel_in,
int nodeID) :
915 OverlapLevel_(OverlapLevel_in)
919 #ifdef IFPACK_OVA_TIME_BUILD 920 double t0 = MPI_Wtime();
921 double t1, true_t0=t0;
925 const int votesMax = INT_MAX;
928 if (OverlapLevel_in == 0)
932 if (
Comm().NumProc() == 1)
939 const Epetra_MpiComm *pComm =
dynamic_cast<const Epetra_MpiComm*
>( &
Comm() );
940 assert(pComm != NULL);
941 MPI_Comm nodeMPIComm;
942 MPI_Comm_split(pComm->Comm(),nodeID,pComm->MyPID(),&nodeMPIComm);
943 Epetra_MpiComm *nodeComm =
new Epetra_MpiComm(nodeMPIComm);
945 Epetra_SerialComm *nodeComm =
dynamic_cast<Epetra_MpiComm*
>( &(Matrix_in->RowMatrixRowMap().Comm()) );
956 const Epetra_Map *RowMap;
957 const Epetra_Map *ColMap;
958 const Epetra_Map *DomainMap;
964 #ifdef IFPACK_OVA_TIME_BUILD 966 fprintf(stderr,
"[node %d]: time for initialization %2.3e\n", nodeID, t1-t0);
971 #ifdef IFPACK_OVA_TIME_BUILD 973 fprintf(stderr,
"[node %d]: overlap hash table allocs %2.3e\n", nodeID, t1-t0);
986 for (
int overlap = 0 ; overlap < OverlapLevel_in ; ++overlap)
991 if (TmpMatrix != Teuchos::null) {
992 RowMap = &(TmpMatrix->RowMatrixRowMap());
993 ColMap = &(TmpMatrix->RowMatrixColMap());
994 DomainMap = &(TmpMatrix->OperatorDomainMap());
997 RowMap = &(
A().RowMatrixRowMap());
998 ColMap = &(
A().RowMatrixColMap());
999 DomainMap = &(
A().OperatorDomainMap());
1002 #ifdef IFPACK_OVA_TIME_BUILD 1004 fprintf(stderr,
"[node %d]: overlap pointer pulls %2.3e\n", nodeID, t1-t0);
1010 Epetra_IntVector colIdList( *ColMap );
1011 Epetra_IntVector rowIdList(*DomainMap);
1012 rowIdList.PutValue(nodeID);
1015 #ifdef IFPACK_OVA_TIME_BUILD 1017 fprintf(stderr,
"[node %d]: overlap intvector/importer allocs %2.3e\n", nodeID, t1-t0);
1021 colIdList.Import(rowIdList,*nodeIdImporter,Insert);
1023 int size = ColMap->NumMyElements() - RowMap->NumMyElements();
1026 #ifdef IFPACK_OVA_TIME_BUILD 1028 fprintf(stderr,
"[node %d]: overlap col/row ID imports %2.3e\n", nodeID, t1-t0);
1035 for (
int i = 0 ; i < ColMap->NumMyElements() ; ++i) {
1036 int GID = ColMap->GID(i);
1038 if (ghostTable.containsKey(GID)) myvotes = ghostTable.get(GID);
1039 if ( colIdList[i] != nodeID && myvotes < votesMax)
1042 if (ghostTable.containsKey(GID)) {
1043 votes = ghostTable.get(GID);
1045 ghostTable.put(GID,votes);
1047 ghostTable.put(GID,1);
1053 ghostTable.arrayify(gidsarray,votesarray);
1066 # ifdef HAVE_MPI //FIXME What if we build in serial?! This file will likely not compile. 1072 #ifdef IFPACK_OVA_TIME_BUILD 1074 fprintf(stderr,
"[node %d]: overlap before culling %2.3e\n", nodeID, t1-t0);
1080 if (nodeComm->MyPID() == 0)
1084 int *lengths =
new int[nodeComm->NumProc()+1];
1086 lengths[1] = ghostTable.size();
1087 for (
int i=1; i<nodeComm->NumProc(); i++) {
1089 MPI_Recv( &leng, 1, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
1090 lengths[i+1] = lengths[i] + leng;
1092 int total = lengths[nodeComm->NumProc()];
1094 int* ghosts =
new int[total];
1095 for (
int i=0; i<total; i++) ghosts[i] = -9;
1096 int *round =
new int[total];
1097 int *owningpid =
new int[total];
1099 for (
int i=1; i<nodeComm->NumProc(); i++) {
1100 int count = lengths[i+1] - lengths[i];
1101 MPI_Recv( ghosts+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
1102 MPI_Recv( round+lengths[i], count, MPI_INT, i, MPI_ANY_TAG, nodeComm->Comm(), &status);
1106 for (
int i=0; i<lengths[1]; i++) {
1107 ghosts[i] = gids[i];
1108 round[i] = votes[i];
1113 for (
int j=1; j<nodeComm->NumProc(); j++) {
1114 for (
int i=lengths[j]; i<lengths[j+1]; i++) {
1121 roundpid[0] = round; roundpid[1] = owningpid;
1122 Epetra_Util epetraUtil;
1123 epetraUtil.Sort(
true,total,ghosts,0,0,2,roundpid);
1126 int *nlosers =
new int[nodeComm->NumProc()];
1127 int** losers =
new int*[nodeComm->NumProc()];
1128 for (
int i=0; i<nodeComm->NumProc(); i++) {
1130 losers[i] =
new int[ lengths[i+1]-lengths[i] ];
1142 for (
int i=1; i<total; i++) {
1143 if (ghosts[i] == ghosts[i-1]) {
1145 if (round[i] > round[
max]) {
1149 int j = owningpid[idx];
1150 losers[j][nlosers[j]++] = ghosts[idx];
1158 delete [] owningpid;
1161 for (
int i=1; i<nodeComm->NumProc(); i++) {
1162 MPI_Send( nlosers+i, 1, MPI_INT, i, 8675, nodeComm->Comm());
1163 MPI_Send( losers[i], nlosers[i], MPI_INT, i, 8675, nodeComm->Comm());
1170 cullList =
new int[ncull+1];
1171 for (
int i=0; i<nlosers[0]; i++)
1172 cullList[i] = losers[0][i];
1174 for (
int i=0; i<nodeComm->NumProc(); i++)
1175 delete [] losers[i];
1184 int hashsize = ghostTable.size();
1185 MPI_Send( &hashsize, 1, MPI_INT, 0, 8675, nodeComm->Comm());
1186 MPI_Send( gids, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
1187 MPI_Send( votes, hashsize, MPI_INT, 0, 8675, nodeComm->Comm());
1191 MPI_Recv( &ncull, 1, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
1192 cullList =
new int[ncull+1];
1193 MPI_Recv( cullList, ncull, MPI_INT, 0, 8675, nodeComm->Comm(), &status);
1200 for (
int i=0; i<ncull; i++) {
1201 try{ghostTable.remove(cullList[i]);}
1204 printf(
"pid %d: In OverlappingRowMatrix ctr, problem removing ghost elt %d from ghostTable\n",
1205 Comm().MyPID(),cullList[i]);
1215 ghostTable.arrayify(gidsarray,votesarray);
1217 std::vector<int> list(size);
1219 for (
int i=0; i<ghostTable.size(); i++) {
1221 if (votesarray[i] < votesMax) {
1222 list[count] = gidsarray[i];
1223 ghostTable.put(gidsarray[i],votesMax);
1229 #ifdef IFPACK_OVA_TIME_BUILD 1231 fprintf(stderr,
"[node %d]: overlap duplicate removal %2.3e\n", nodeID, t1-t0);
1236 # endif //ifdef HAVE_MPI 1238 TmpMap =
rcp(
new Epetra_Map(-1,count, &list[0],0,
Comm()) );
1240 TmpMatrix =
rcp(
new Epetra_CrsMatrix(Copy,*TmpMap,0) );
1244 TmpMatrix->Import(
A(),*TmpImporter,Insert);
1248 #ifdef IFPACK_OVA_TIME_BUILD 1250 fprintf(stderr,
"[node %d]: overlap TmpMatrix fillcomplete %2.3e\n", nodeID, t1-t0);
1261 Epetra_IntVector ov_colIdList( TmpMatrix->ColMap() );
1262 ov_colIdList.PutValue(-1);
1263 Epetra_IntVector ov_rowIdList( TmpMatrix->RowMap() );
1264 ov_rowIdList.PutValue(nodeID);
1266 ov_colIdList.Import(ov_rowIdList,*ov_nodeIdImporter,Insert);
1268 for (
int i=0 ; i<ov_colIdList.MyLength(); i++) {
1269 if (ov_colIdList[i] == nodeID) {
1270 int GID = (ov_colIdList.Map()).GID(i);
1271 colMapTable.put(GID,1);
1278 ov_colIdList.PutValue(-1);
1280 ArowIdList.PutValue(nodeID);
1281 nodeIdImporter =
rcp(
new Epetra_Import( TmpMatrix->ColMap(),
A().RowMatrixRowMap() ));
1282 ov_colIdList.Import(ArowIdList,*nodeIdImporter,Insert);
1284 for (
int i=0 ; i<ov_colIdList.MyLength(); i++) {
1285 if (ov_colIdList[i] == nodeID) {
1286 int GID = (ov_colIdList.Map()).GID(i);
1287 colMapTable.put(GID,1);
1292 #ifdef IFPACK_OVA_TIME_BUILD 1294 fprintf(stderr,
"[node %d]: overlap 2 imports to fix up colmap %2.3e\n", nodeID, t1-t0);
1306 std::vector<int> ghostElements;
1309 ghostTable.arrayify(gidsarray,votesarray);
1310 for (
int i=0; i<ghostTable.size(); i++) {
1314 for (
int i = 0 ; i <
A().RowMatrixColMap().NumMyElements() ; ++i) {
1315 int GID =
A().RowMatrixColMap().GID(i);
1317 if (colMapTable.containsKey(GID)) {
1318 try{colMapTable.remove(GID);}
1320 printf(
"pid %d: In OverlappingRowMatrix ctr, problem removing entry %d from colMapTable\n",
Comm().MyPID(),GID);
1327 std::vector<int> colMapElements;
1330 colMapTable.arrayify(gidsarray,votesarray);
1331 for (
int i=0; i<colMapTable.size(); i++)
1332 colMapElements.push_back(gidsarray[i]);
1345 std::vector<int> rowList(
NumMyRowsA_ + ghostElements.size());
1347 rowList[i] =
A().RowMatrixRowMap().GID(i);
1348 for (
int i = 0 ; i < (int)ghostElements.size() ; ++i)
1360 std::vector<int> colList(
A().
RowMatrixColMap().NumMyElements() + colMapElements.size());
1361 int nc =
A().RowMatrixColMap().NumMyElements();
1362 for (
int i = 0 ; i < nc; i++)
1363 colList[i] =
A().RowMatrixColMap().GID(i);
1364 for (
int i = 0 ; i < (int)colMapElements.size() ; i++)
1365 colList[nc+i] = colMapElements[i];
1369 colMap_ =
new Epetra_Map(-1,
A().
RowMatrixColMap().NumMyElements() + colMapElements.size(), &colList[0], 0,
Comm());
1373 #ifdef IFPACK_OVA_TIME_BUILD 1375 fprintf(stderr,
"[node %d]: time to init B() col/row maps %2.3e\n", nodeID, t1-t0);
1387 Epetra_Map* nodeMap =
new Epetra_Map(-1,
NumMyRowsA_ + ghostElements.size(),&rowList[0],0,*nodeComm);
1388 int numMyElts = colMap_->NumMyElements();
1389 assert(numMyElts!=0);
1395 int* myGlobalElts =
new int[numMyElts];
1396 colMap_->MyGlobalElements(myGlobalElts);
1397 int* pidList =
new int[numMyElts];
1398 nodeMap->RemoteIDList(numMyElts, myGlobalElts, pidList, 0);
1403 tt[0] = myGlobalElts;
1404 Util.Sort(
true, numMyElts, pidList, 0, (
double**)0, 1, tt);
1409 while (localStart<numMyElts) {
1410 int currPID = (pidList)[localStart];
1412 while (i<numMyElts && (pidList)[i] == currPID) i++;
1413 if (currPID != nodeComm->MyPID())
1414 Util.Sort(
true, i-localStart, (myGlobalElts)+localStart, 0, 0, 0, 0);
1420 while (localStart<numMyElts && (pidList)[localStart] != nodeComm->MyPID()) localStart++;
1421 assert(localStart != numMyElts);
1422 int localEnd=localStart;
1423 while (localEnd<numMyElts && (pidList)[localEnd] == nodeComm->MyPID()) localEnd++;
1424 int* mySortedGlobalElts =
new int[numMyElts];
1425 int leng = localEnd - localStart;
1429 int *rowGlobalElts = nodeMap->MyGlobalElements();
1449 for (
int i=0; i<leng; i++) {
1451 (mySortedGlobalElts)[i] = rowGlobalElts[i];
1455 for (
int i=leng; i< localEnd; i++) {
1456 (mySortedGlobalElts)[i] = (myGlobalElts)[i-leng];
1458 for (
int i=localEnd; i<numMyElts; i++) {
1459 (mySortedGlobalElts)[i] = (myGlobalElts)[i];
1463 #ifdef IFPACK_OVA_TIME_BUILD 1465 fprintf(stderr,
"[node %d]: time to sort col map arrays (cp=1) %2.3e\n", nodeID, t1-t0);
1470 int indexBase = colMap_->IndexBase();
1471 if (colMap_)
delete colMap_;
1472 colMap_ =
new Epetra_Map(-1,numMyElts,mySortedGlobalElts,indexBase,
Comm());
1475 delete [] myGlobalElts;
1477 delete [] mySortedGlobalElts;
1481 printf(
"** * Ifpack_OverlappingRowmatrix ctor: problem creating column map * **\n\n");
1489 #ifdef IFPACK_OVA_TIME_BUILD 1491 fprintf(stderr,
"[node %d]: time to sort col map arrays (cp=2) %2.3e\n", nodeID, t1-t0);
1521 ExtMap_ =
rcp(
new Epetra_Map(-1,ghostElements.size(), &ghostElements[0],0,
Comm()) );
1528 #ifdef IFPACK_OVA_TIME_BUILD 1530 fprintf(stderr,
"[node %d]: time to import and setup ExtMap_ %2.3e\n", nodeID, t1-t0);
1546 #ifdef IFPACK_OVA_TIME_BUILD 1548 fprintf(stderr,
"[node %d]: time to FillComplete B() %2.3e\n", nodeID, t1-t0);
1581 #ifdef IFPACK_OVA_TIME_BUILD 1583 fprintf(stderr,
"[node %d]: time for final calcs %2.3e\n", nodeID, t1-t0);
1584 fprintf(stderr,
"[node %d]: total IORM ctor time %2.3e\n", nodeID, t1-true_t0);
1596 # endif //ifdef IFPACK_NODE_AWARE_CODE 1597 #endif // ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 1610 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 1614 int * Indices)
const 1618 const Epetra_Map *Themap;
1620 ierr =
A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
1621 Themap=&
A().RowMatrixColMap();
1624 ierr =
B().ExtractMyRowCopy(LocRow-
NumMyRowsA_,Length,NumEntries,Values,Indices);
1625 Themap=&
B().RowMatrixColMap();
1632 int Ifpack_OverlappingRowMatrix::
1633 ExtractGlobalRowCopy(
int GlobRow,
int Length,
int & NumEntries,
double *Values,
1634 int * Indices)
const 1637 const Epetra_Map *Themap;
1638 int LocRow =
A().RowMatrixRowMap().LID(GlobRow);
1640 ierr =
A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
1641 Themap=&
A().RowMatrixColMap();
1644 LocRow =
B().RowMatrixRowMap().LID(GlobRow);
1647 ierr =
B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
1648 Themap=&
B().RowMatrixColMap();
1651 for (
int i=0; i<NumEntries; i++) {
1652 Indices[i]=Themap->GID(Indices[i]);
1653 assert(Indices[i] != -1);
1659 # ifdef IFPACK_NODE_AWARE_CODE 1663 int * Indices)
const 1667 const Epetra_Map *Themap;
1669 ierr =
A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
1670 Themap=&
A().RowMatrixColMap();
1673 ierr =
B().ExtractMyRowCopy(LocRow-
NumMyRowsA_,Length,NumEntries,Values,Indices);
1674 Themap=&
B().RowMatrixColMap();
1681 int Ifpack_OverlappingRowMatrix::
1682 ExtractGlobalRowCopy(
int GlobRow,
int Length,
int & NumEntries,
double *Values,
1683 int * Indices)
const 1686 const Epetra_Map *Themap;
1687 int LocRow =
A().RowMatrixRowMap().LID(GlobRow);
1689 ierr =
A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
1690 Themap=&
A().RowMatrixColMap();
1693 LocRow =
B().RowMatrixRowMap().LID(GlobRow);
1696 ierr =
B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
1697 Themap=&
B().RowMatrixColMap();
1700 for (
int i=0; i<NumEntries; i++) {
1701 Indices[i]=Themap->GID(Indices[i]);
1702 assert(Indices[i] != -1);
1712 int * Indices)
const 1716 ierr =
A().ExtractMyRowCopy(MyRow,Length,NumEntries,Values,Indices);
1718 ierr =
B().ExtractMyRowCopy(MyRow -
NumMyRowsA_,Length,NumEntries,
1723 # endif // ifdef IFPACK_NODE_AWARE_CODE 1724 #endif // ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 1736 Multiply(
bool TransA,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1750 for (
int j = 0 ; j < Nnz ; ++j) {
1751 Y[k][i] += Val[j] * X[k][Ind[j]];
1762 for (
int j = 0 ; j < Nnz ; ++j) {
1772 Apply(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const 1786 #ifndef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS 1787 # ifndef IFPACK_NODE_AWARE_CODE 1803 Epetra_CombineMode CM)
1812 Epetra_CombineMode CM)
void BuildMap(int OverlapLevel_in)
Teuchos::RefCountPtr< Epetra_CrsMatrix > ExtMatrix_
~Ifpack_OverlappingRowMatrix()
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
Epetra_RowMatrix & B() const
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Teuchos::RefCountPtr< const Epetra_Map > Map_
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
const Epetra_BlockMap & Map() const
Teuchos::RefCountPtr< Epetra_Import > ExtImporter_
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
long long NumGlobalNonzeros_
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
#define IFPACK_CHK_ERRV(ifpack_err)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual const Epetra_Map & RowMatrixRowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
const Epetra_RowMatrix & A() const
virtual const Epetra_Map & RowMatrixColMap() const
Returns the Epetra_Map object associated with the columns of this matrix.
#define IFPACK_RETURN(ifpack_err)
void push_back(const value_type &x)
Teuchos::RefCountPtr< Epetra_Map > ExtMap_
int ImportMultiVector(const Epetra_MultiVector &X, Epetra_MultiVector &OvX, Epetra_CombineMode CM=Insert)
virtual int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
int ExportMultiVector(const Epetra_MultiVector &OvX, Epetra_MultiVector &X, Epetra_CombineMode CM=Add)
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const
Returns the number of nonzero entries in MyRow.
#define IFPACK_CHK_ERR(ifpack_err)
bool UseTranspose() const
Returns the current UseTranspose setting.
Ifpack_OverlappingRowMatrix(const Teuchos::RefCountPtr< const Epetra_RowMatrix > &Matrix_in, int OverlapLevel_in)
Teuchos::RefCountPtr< const Epetra_Import > Importer_
static bool find(Parser_dh p, char *option, OptionsNode **ptr)
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const