EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
inout_LL/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//==============================================================================
52Poisson2dOperator::Poisson2dOperator(int nx, int ny, const Epetra_Comm & comm)
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_(0) {
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(-1LL, ((long long)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 long long[numImports_];
88 long long * ptr = importIDs_;
89 long long minGID = map_->MinMyGID64();
90 long long maxGID = map_->MaxMyGID64();
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
100 importMap_ = new Epetra_Map(-1LL, numImports_, importIDs_, 0LL, comm_);
101
102 importer_ = new Epetra_Import(*importMap_, *map_);
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)
128 else if (importX_->NumVectors()!=X.NumVectors()) {
129 delete importX_;
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
203 Epetra_CrsMatrix * A = new Epetra_CrsMatrix(Copy, *map_, 3);
204
205 int NumMyElements = map_->NumMyElements();
206 long long NumGlobalElements = map_->NumGlobalElements64();
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 long long GlobalRow = A->GRID64(i); long long RowLess1 = GlobalRow - 1; long long 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}
long long NumGlobalElements64() const
int NumMyElements() const
int FillComplete(bool OptimizeDataStorage=true)
long long GRID64(int LRID_in) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_BlockMap & Map() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int NumVectors() const
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
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.
Epetra_CrsMatrix * GeneratePrecMatrix() const
Generate a tridiagonal approximation to the 5-point Poisson as an Epetra_CrsMatrix.
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_Comm & comm_
Epetra_MultiVector * importX_