EpetraExt Package Browser (Single Doxygen Collection) Development
Poisson2dOperator.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "Poisson2dOperator.h"
00043 #include "Epetra_CrsMatrix.h"
00044 #include "Epetra_Map.h"
00045 #include "Epetra_Import.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_MultiVector.h"
00048 #include "Epetra_Comm.h"
00049 #include "Epetra_Distributor.h"
00050 
00051 //==============================================================================
00052 Poisson2dOperator::Poisson2dOperator(int nx, int ny, const Epetra_Comm & comm) 
00053   : nx_(nx),
00054     ny_(ny),
00055     useTranspose_(false),
00056     comm_(comm),
00057     map_(0),
00058     numImports_(0),
00059     importIDs_(0),
00060     importMap_(0),
00061     importer_(0),
00062     importX_(0),
00063     Label_(0) {
00064 
00065   Label_ = "2D Poisson Operator";
00066   int numProc = comm.NumProc(); // Get number of processors
00067   int myPID = comm.MyPID(); // My rank
00068   if (2*numProc > ny) { // ny must be >= 2*numProc (to avoid degenerate cases)
00069     ny = 2*numProc; ny_ = ny;
00070     std::cout << " Increasing ny to " << ny << " to avoid degenerate distribution on " << numProc << " processors." << std::endl;
00071   }
00072   
00073   int chunkSize = ny/numProc;
00074   int remainder = ny%numProc;
00075 
00076   if (myPID+1 <= remainder) chunkSize++; // add on remainder
00077 
00078   myny_ = chunkSize;
00079 
00080   map_ = new Epetra_Map(-1, nx*chunkSize, 0, comm_);
00081 
00082   if (numProc>1) {
00083     // Build import GID list to build import map and importer
00084     if (myPID>0) numImports_ += nx;
00085     if (myPID+1<numProc) numImports_ += nx;
00086     
00087     if (numImports_>0) importIDs_ = new int[numImports_];
00088     int * ptr = importIDs_;
00089     int minGID = map_->MinMyGID();
00090     int maxGID = map_->MaxMyGID();
00091     
00092     if (myPID>0) for (int i=0; i< nx; i++) *ptr++ = minGID - nx + i;
00093     if (myPID+1<numProc) for (int i=0; i< nx; i++) *ptr++ = maxGID + i +1;
00094     
00095     // At the end of the above step importIDs_ will have a list of global IDs that are needed
00096     // to compute the matrix multiplication operation on this processor.  Now build import map 
00097     // and importer
00098     
00099     
00100     importMap_ = new Epetra_Map(-1, numImports_, importIDs_, 0, comm_);
00101     
00102     importer_ = new Epetra_Import(*importMap_, *map_);
00103 
00104   }
00105 }
00106 //==============================================================================
00107 Poisson2dOperator::~Poisson2dOperator() {
00108   if (importX_!=0) delete importX_;
00109   if (importer_!=0) delete importer_;
00110   if (importMap_!=0) delete importMap_;
00111   if (importIDs_!=0) delete [] importIDs_;
00112   if (map_!=0) delete map_;
00113 }
00114 //==============================================================================
00115 int Poisson2dOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00116 
00117 
00118   // This is a very brain-dead implementation of a 5-point finite difference stencil, but
00119   // it should illustrate the basic process for implementing the Epetra_Operator interface.
00120 
00121   if (!X.Map().SameAs(OperatorDomainMap())) abort();  // These aborts should be handled as int return codes.
00122   if (!Y.Map().SameAs(OperatorRangeMap())) abort();
00123   if (Y.NumVectors()!=X.NumVectors()) abort();
00124 
00125   if (comm_.NumProc()>1) {
00126     if (importX_==0)
00127       importX_ = new Epetra_MultiVector(*importMap_, X.NumVectors());
00128     else if (importX_->NumVectors()!=X.NumVectors()) {
00129       delete importX_;
00130       importX_ = new Epetra_MultiVector(*importMap_, X.NumVectors());
00131     }
00132     importX_->Import(X, *importer_, Insert); // Get x values we need
00133   }
00134 
00135   double * importx1 = 0;
00136   double * importx2 = 0;
00137   int nx = nx_;
00138   //int ny = ny_;
00139 
00140   for (int j=0; j < X.NumVectors(); j++) {
00141 
00142     const double * x = X[j];
00143     if (comm_.NumProc()>1) {
00144       importx1 = (*importX_)[j];
00145       importx2 = importx1+nx;
00146       if (comm_.MyPID()==0) importx2=importx1;
00147     }
00148     double * y = Y[j];
00149     if (comm_.MyPID()==0) {
00150       y[0] = 4.0*x[0]-x[nx]-x[1];
00151       y[nx-1] = 4.0*x[nx-1]-x[nx-2]-x[nx+nx-1];
00152       for (int ix=1; ix< nx-1; ix++)
00153   y[ix] = 4.0*x[ix]-x[ix-1]-x[ix+1]-x[ix+nx];
00154     }
00155     else {
00156       y[0] = 4.0*x[0]-x[nx]-x[1]-importx1[0];
00157       y[nx-1] = 4.0*x[nx-1]-x[nx-2]-x[nx+nx-1]-importx1[nx-1];
00158       for (int ix=1; ix< nx-1; ix++)
00159   y[ix] = 4.0*x[ix]-x[ix-1]-x[ix+1]-x[ix+nx]-importx1[ix];
00160     }
00161     if (comm_.MyPID()+1==comm_.NumProc()) {
00162       int curxy = nx*myny_-1;
00163       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1];
00164       curxy -= (nx-1);
00165       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1];
00166       for (int ix=1; ix< nx-1; ix++) {
00167   curxy++;
00168   y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx];
00169       }
00170     }
00171     else {
00172       int curxy = nx*myny_-1;
00173       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-importx2[nx-1];
00174       curxy -= (nx-1);
00175       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1]-importx2[0];
00176       for (int ix=1; ix< nx-1; ix++) {
00177   curxy++;
00178   y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-importx2[ix];
00179       }
00180     }
00181     for (int iy=1; iy< myny_-1; iy++) {
00182       int curxy = nx*(iy+1)-1;
00183       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-x[curxy+nx];
00184       curxy -= (nx-1);
00185       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1]-x[curxy+nx];
00186       for (int ix=1; ix< nx-1; ix++) {
00187   curxy++;
00188   y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-x[curxy+nx];
00189       }
00190     }
00191   }
00192   return(0);
00193 }
00194 //==============================================================================
00195 Epetra_CrsMatrix * Poisson2dOperator::GeneratePrecMatrix() const {
00196 
00197   // Generate a tridiagonal matrix as an Epetra_CrsMatrix
00198   // This method illustrates how to generate a matrix that is an approximation to the true
00199   // operator.  Given this matrix, we can use any of the Aztec or IFPACK preconditioners.
00200 
00201 
00202   // Create a Epetra_Matrix
00203   Epetra_CrsMatrix * A = new Epetra_CrsMatrix(Copy, *map_, 3);
00204 
00205   int NumMyElements = map_->NumMyElements();
00206   int NumGlobalElements = map_->NumGlobalElements();
00207 
00208   // Add  rows one-at-a-time
00209   double negOne = -1.0;
00210   double posTwo = 4.0;
00211   for (int i=0; i<NumMyElements; i++) {
00212     int GlobalRow = A->GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1;
00213 
00214     if (RowLess1!=-1) A->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
00215     if (RowPlus1!=NumGlobalElements) A->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
00216     A->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
00217   }
00218 
00219   // Finish up
00220   A->FillComplete();
00221 
00222   return(A);
00223 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines