EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Poisson2dOperator.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 
42 #include "Poisson2dOperator.h"
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"
50 
51 //==============================================================================
53  : nx_(nx),
54  ny_(ny),
55  useTranspose_(false),
56  comm_(comm),
57  map_(0),
58  numImports_(0),
59  importIDs_(0),
60  importMap_(0),
61  importer_(0),
62  importX_(0),
63  Label_ ("") {
64 
65  Label_ = "2D Poisson Operator";
66  int numProc = comm.NumProc(); // Get number of processors
67  int myPID = comm.MyPID(); // My rank
68  if (2*numProc > ny) { // ny must be >= 2*numProc (to avoid degenerate cases)
69  ny = 2*numProc; ny_ = ny;
70  std::cout << " Increasing ny to " << ny << " to avoid degenerate distribution on " << numProc << " processors." << std::endl;
71  }
72 
73  int chunkSize = ny/numProc;
74  int remainder = ny%numProc;
75 
76  if (myPID+1 <= remainder) chunkSize++; // add on remainder
77 
78  myny_ = chunkSize;
79 
80  map_ = new Epetra_Map(-1, nx*chunkSize, 0, comm_);
81 
82  if (numProc>1) {
83  // Build import GID list to build import map and importer
84  if (myPID>0) numImports_ += nx;
85  if (myPID+1<numProc) numImports_ += nx;
86 
87  if (numImports_>0) importIDs_ = new int[numImports_];
88  int * ptr = importIDs_;
89  int minGID = map_->MinMyGID();
90  int maxGID = map_->MaxMyGID();
91 
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;
94 
95  // At the end of the above step importIDs_ will have a list of global IDs that are needed
96  // to compute the matrix multiplication operation on this processor. Now build import map
97  // and importer
98 
99 
101 
103 
104  }
105 }
106 //==============================================================================
108  if (importX_!=0) delete importX_;
109  if (importer_!=0) delete importer_;
110  if (importMap_!=0) delete importMap_;
111  if (importIDs_!=0) delete [] importIDs_;
112  if (map_!=0) delete map_;
113 }
114 //==============================================================================
116 
117 
118  // This is a very brain-dead implementation of a 5-point finite difference stencil, but
119  // it should illustrate the basic process for implementing the Epetra_Operator interface.
120 
121  if (!X.Map().SameAs(OperatorDomainMap())) abort(); // These aborts should be handled as int return codes.
122  if (!Y.Map().SameAs(OperatorRangeMap())) abort();
123  if (Y.NumVectors()!=X.NumVectors()) abort();
124 
125  if (comm_.NumProc()>1) {
126  if (importX_==0)
127  importX_ = new Epetra_MultiVector(*importMap_, X.NumVectors());
128  else if (importX_->NumVectors()!=X.NumVectors()) {
129  delete importX_;
130  importX_ = new Epetra_MultiVector(*importMap_, X.NumVectors());
131  }
132  importX_->Import(X, *importer_, Insert); // Get x values we need
133  }
134 
135  double * importx1 = 0;
136  double * importx2 = 0;
137  int nx = nx_;
138  //int ny = ny_;
139 
140  for (int j=0; j < X.NumVectors(); j++) {
141 
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;
147  }
148  double * y = Y[j];
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];
154  }
155  else {
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];
160  }
161  if (comm_.MyPID()+1==comm_.NumProc()) {
162  int curxy = nx*myny_-1;
163  y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1];
164  curxy -= (nx-1);
165  y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1];
166  for (int ix=1; ix< nx-1; ix++) {
167  curxy++;
168  y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx];
169  }
170  }
171  else {
172  int curxy = nx*myny_-1;
173  y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-importx2[nx-1];
174  curxy -= (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++) {
177  curxy++;
178  y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-importx2[ix];
179  }
180  }
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];
184  curxy -= (nx-1);
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++) {
187  curxy++;
188  y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-x[curxy+nx];
189  }
190  }
191  }
192  return(0);
193 }
194 //==============================================================================
196 
197  // Generate a tridiagonal matrix as an Epetra_CrsMatrix
198  // This method illustrates how to generate a matrix that is an approximation to the true
199  // operator. Given this matrix, we can use any of the Aztec or IFPACK preconditioners.
200 
201 
202  // Create a Epetra_Matrix
204 
205  int NumMyElements = map_->NumMyElements();
206  int NumGlobalElements = map_->NumGlobalElements();
207 
208  // Add rows one-at-a-time
209  double negOne = -1.0;
210  double posTwo = 4.0;
211  for (int i=0; i<NumMyElements; i++) {
212  int GlobalRow = A->GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1;
213 
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);
217  }
218 
219  // Finish up
220  A->FillComplete();
221 
222  return(A);
223 }
int NumGlobalElements() const
bool SameAs(const Epetra_BlockMap &Map) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_Map * importMap_
~Poisson2dOperator()
Destructor.
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Poisson2dOperator applied to a Epetra_MultiVector X in Y. ...
int NumMyElements() const
const Epetra_Comm & comm_
virtual const Epetra_BlockMap & Map() const =0
Epetra_CrsMatrix * GeneratePrecMatrix() const
Generate a tridiagonal approximation to the 5-point Poisson as an Epetra_CrsMatrix.
int MaxMyGID() const
int MinMyGID() const
Epetra_Import * importer_
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Epetra_MultiVector * importX_
virtual int NumProc() const =0
int GRID(int LRID_in) const
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 & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.