EpetraExt Package Browser (Single Doxygen Collection) Development
Poisson2dOperator.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //        AztecOO: An Object-Oriented Aztec Linear Solver Package 
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #include "Poisson2dOperator.h"
00030 #include "Epetra_CrsMatrix.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_Import.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_MultiVector.h"
00035 #include "Epetra_Comm.h"
00036 #include "Epetra_Distributor.h"
00037 
00038 //==============================================================================
00039 Poisson2dOperator::Poisson2dOperator(int nx, int ny, const Epetra_Comm & comm) 
00040   : nx_(nx),
00041     ny_(ny),
00042     useTranspose_(false),
00043     comm_(comm),
00044     map_(0),
00045     numImports_(0),
00046     importIDs_(0),
00047     importMap_(0),
00048     importer_(0),
00049     importX_(0),
00050     Label_(0) {
00051 
00052   Label_ = "2D Poisson Operator";
00053   int numProc = comm.NumProc(); // Get number of processors
00054   int myPID = comm.MyPID(); // My rank
00055   if (2*numProc > ny) { // ny must be >= 2*numProc (to avoid degenerate cases)
00056     ny = 2*numProc; ny_ = ny;
00057     cout << " Increasing ny to " << ny << " to avoid degenerate distribution on " << numProc << " processors." << endl;
00058   }
00059   
00060   int chunkSize = ny/numProc;
00061   int remainder = ny%numProc;
00062 
00063   if (myPID+1 <= remainder) chunkSize++; // add on remainder
00064 
00065   myny_ = chunkSize;
00066 
00067   map_ = new Epetra_Map(-1, nx*chunkSize, 0, comm_);
00068 
00069   if (numProc>1) {
00070     // Build import GID list to build import map and importer
00071     if (myPID>0) numImports_ += nx;
00072     if (myPID+1<numProc) numImports_ += nx;
00073     
00074     if (numImports_>0) importIDs_ = new int[numImports_];
00075     int * ptr = importIDs_;
00076     int minGID = map_->MinMyGID();
00077     int maxGID = map_->MaxMyGID();
00078     
00079     if (myPID>0) for (int i=0; i< nx; i++) *ptr++ = minGID - nx + i;
00080     if (myPID+1<numProc) for (int i=0; i< nx; i++) *ptr++ = maxGID + i +1;
00081     
00082     // At the end of the above step importIDs_ will have a list of global IDs that are needed
00083     // to compute the matrix multiplication operation on this processor.  Now build import map 
00084     // and importer
00085     
00086     
00087     importMap_ = new Epetra_Map(-1, numImports_, importIDs_, 0, comm_);
00088     
00089     importer_ = new Epetra_Import(*importMap_, *map_);
00090 
00091   }
00092 }
00093 //==============================================================================
00094 Poisson2dOperator::~Poisson2dOperator() {
00095   if (importX_!=0) delete importX_;
00096   if (importer_!=0) delete importer_;
00097   if (importMap_!=0) delete importMap_;
00098   if (importIDs_!=0) delete [] importIDs_;
00099   if (map_!=0) delete map_;
00100 }
00101 //==============================================================================
00102 int Poisson2dOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00103 
00104 
00105   // This is a very brain-dead implementation of a 5-point finite difference stencil, but
00106   // it should illustrate the basic process for implementing the Epetra_Operator interface.
00107 
00108   if (!X.Map().SameAs(OperatorDomainMap())) abort();  // These aborts should be handled as int return codes.
00109   if (!Y.Map().SameAs(OperatorRangeMap())) abort();
00110   if (Y.NumVectors()!=X.NumVectors()) abort();
00111 
00112   if (comm_.NumProc()>1) {
00113     if (importX_==0)
00114       importX_ = new Epetra_MultiVector(*importMap_, X.NumVectors());
00115     else if (importX_->NumVectors()!=X.NumVectors()) {
00116       delete importX_;
00117       importX_ = new Epetra_MultiVector(*importMap_, X.NumVectors());
00118     }
00119     importX_->Import(X, *importer_, Insert); // Get x values we need
00120   }
00121 
00122   double * importx1 = 0;
00123   double * importx2 = 0;
00124   int nx = nx_;
00125   //int ny = ny_;
00126 
00127   for (int j=0; j < X.NumVectors(); j++) {
00128 
00129     const double * x = X[j];
00130     if (comm_.NumProc()>1) {
00131       importx1 = (*importX_)[j];
00132       importx2 = importx1+nx;
00133       if (comm_.MyPID()==0) importx2=importx1;
00134     }
00135     double * y = Y[j];
00136     if (comm_.MyPID()==0) {
00137       y[0] = 4.0*x[0]-x[nx]-x[1];
00138       y[nx-1] = 4.0*x[nx-1]-x[nx-2]-x[nx+nx-1];
00139       for (int ix=1; ix< nx-1; ix++)
00140   y[ix] = 4.0*x[ix]-x[ix-1]-x[ix+1]-x[ix+nx];
00141     }
00142     else {
00143       y[0] = 4.0*x[0]-x[nx]-x[1]-importx1[0];
00144       y[nx-1] = 4.0*x[nx-1]-x[nx-2]-x[nx+nx-1]-importx1[nx-1];
00145       for (int ix=1; ix< nx-1; ix++)
00146   y[ix] = 4.0*x[ix]-x[ix-1]-x[ix+1]-x[ix+nx]-importx1[ix];
00147     }
00148     if (comm_.MyPID()+1==comm_.NumProc()) {
00149       int curxy = nx*myny_-1;
00150       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1];
00151       curxy -= (nx-1);
00152       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1];
00153       for (int ix=1; ix< nx-1; ix++) {
00154   curxy++;
00155   y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx];
00156       }
00157     }
00158     else {
00159       int curxy = nx*myny_-1;
00160       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-importx2[nx-1];
00161       curxy -= (nx-1);
00162       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1]-importx2[0];
00163       for (int ix=1; ix< nx-1; ix++) {
00164   curxy++;
00165   y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-importx2[ix];
00166       }
00167     }
00168     for (int iy=1; iy< myny_-1; iy++) {
00169       int curxy = nx*(iy+1)-1;
00170       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy-1]-x[curxy+nx];
00171       curxy -= (nx-1);
00172       y[curxy] = 4.0*x[curxy]-x[curxy-nx]-x[curxy+1]-x[curxy+nx];
00173       for (int ix=1; ix< nx-1; ix++) {
00174   curxy++;
00175   y[curxy] = 4.0*x[curxy]-x[curxy-1]-x[curxy+1]-x[curxy-nx]-x[curxy+nx];
00176       }
00177     }
00178   }
00179   return(0);
00180 }
00181 //==============================================================================
00182 Epetra_CrsMatrix * Poisson2dOperator::GeneratePrecMatrix() const {
00183 
00184   // Generate a tridiagonal matrix as an Epetra_CrsMatrix
00185   // This method illustrates how to generate a matrix that is an approximation to the true
00186   // operator.  Given this matrix, we can use any of the Aztec or IFPACK preconditioners.
00187 
00188 
00189   // Create a Epetra_Matrix
00190   Epetra_CrsMatrix * A = new Epetra_CrsMatrix(Copy, *map_, 3);
00191 
00192   int NumMyElements = map_->NumMyElements();
00193   int NumGlobalElements = map_->NumGlobalElements();
00194 
00195   // Add  rows one-at-a-time
00196   double negOne = -1.0;
00197   double posTwo = 4.0;
00198   for (int i=0; i<NumMyElements; i++) {
00199     int GlobalRow = A->GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1;
00200 
00201     if (RowLess1!=-1) A->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
00202     if (RowPlus1!=NumGlobalElements) A->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
00203     A->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
00204   }
00205 
00206   // Finish up
00207   A->FillComplete();
00208 
00209   return(A);
00210 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines