43 #include "Epetra_CrsMatrix.h" 44 #include "Epetra_Map.h" 45 #include "Epetra_Import.h" 46 #include "Epetra_Vector.h" 47 #include "Epetra_MultiVector.h" 48 #include "Epetra_Comm.h" 49 #include "Epetra_Distributor.h" 65 Label_ =
"2D Poisson Operator";
66 int numProc = comm.NumProc();
67 int myPID = comm.MyPID();
69 ny = 2*numProc;
ny_ = ny;
70 std::cout <<
" Increasing ny to " << ny <<
" to avoid degenerate distribution on " << numProc <<
" processors." << std::endl;
73 int chunkSize = ny/numProc;
74 int remainder = ny%numProc;
76 if (myPID+1 <= remainder) chunkSize++;
80 map_ =
new Epetra_Map(-1LL, ((
long long)nx)*chunkSize, 0,
comm_);
89 long long minGID =
map_->MinMyGID64();
90 long long maxGID =
map_->MaxMyGID64();
92 if (myPID>0)
for (
int i=0; i< nx; i++) *ptr++ = minGID - nx + i;
93 if (myPID+1<numProc)
for (
int i=0; i< nx; i++) *ptr++ = maxGID + i +1;
123 if (Y.NumVectors()!=X.NumVectors()) abort();
125 if (
comm_.NumProc()>1) {
128 else if (
importX_->NumVectors()!=X.NumVectors()) {
135 double * importx1 = 0;
136 double * importx2 = 0;
140 for (
int j=0; j < X.NumVectors(); j++) {
142 const double * x = X[j];
143 if (
comm_.NumProc()>1) {
144 importx1 = (*importX_)[j];
145 importx2 = importx1+nx;
146 if (
comm_.MyPID()==0) importx2=importx1;
149 if (
comm_.MyPID()==0) {
150 y[0] = 4.0*x[0]-x[nx]-x[1];
151 y[nx-1] = 4.0*x[nx-1]-x[nx-2]-x[nx+nx-1];
152 for (
int ix=1; ix< nx-1; ix++)
153 y[ix] = 4.0*x[ix]-x[ix-1]-x[ix+1]-x[ix+nx];
156 y[0] = 4.0*x[0]-x[nx]-x[1]-importx1[0];
157 y[nx-1] = 4.0*x[nx-1]-x[nx-2]-x[nx+nx-1]-importx1[nx-1];
158 for (
int ix=1; ix< nx-1; ix++)
159 y[ix] = 4.0*x[ix]-x[ix-1]-x[ix+1]-x[ix+nx]-importx1[ix];
162 int curxy = nx*
myny_-1;
163 y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1];
165 y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1];
166 for (
int ix=1; ix< nx-1; ix++) {
168 y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx];
172 int curxy = nx*
myny_-1;
173 y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-importx2[nx-1];
175 y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1]-importx2[0];
176 for (
int ix=1; ix< nx-1; ix++) {
178 y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-importx2[ix];
181 for (
int iy=1; iy<
myny_-1; iy++) {
182 int curxy = nx*(iy+1)-1;
183 y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-x[curxy+nx];
185 y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1]-x[curxy+nx];
186 for (
int ix=1; ix< nx-1; ix++) {
188 y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-x[curxy+nx];
203 Epetra_CrsMatrix * A =
new Epetra_CrsMatrix(Copy, *
map_, 3);
205 int NumMyElements =
map_->NumMyElements();
206 long long NumGlobalElements =
map_->NumGlobalElements64();
209 double negOne = -1.0;
211 for (
int i=0; i<NumMyElements; i++) {
212 long long GlobalRow = A->GRID64(i);
long long RowLess1 = GlobalRow - 1;
long long RowPlus1 = GlobalRow + 1;
214 if (RowLess1!=-1) A->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
215 if (RowPlus1!=NumGlobalElements) A->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
216 A->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
Epetra_CrsMatrix * GeneratePrecMatrix() const
Generate a tridiagonal approximation to the 5-point Poisson as an Epetra_CrsMatrix.
~Poisson2dOperator()
Destructor.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Poisson2dOperator applied to a Epetra_MultiVector X in Y. ...
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
const Epetra_Comm & comm_
Epetra_Import * importer_
Epetra_MultiVector * importX_
Poisson2dOperator(int nx, int ny, const Epetra_Comm &comm)
Builds a 2 dimensional Poisson operator for a nx by ny grid, assuming zero Dirichlet BCs...
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.