Ifpack_LocalFilter.cpp

00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner 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 
00030 #include "Ifpack_ConfigDefs.h"
00031 
00032 #include "Epetra_MultiVector.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_RowMatrix.h"
00035 #include "Epetra_Map.h"
00036 #include "Epetra_BlockMap.h"
00037 #include "Ifpack_LocalFilter.h"
00038 
00039 //==============================================================================
00040 Ifpack_LocalFilter::Ifpack_LocalFilter(const Teuchos::RefCountPtr<const Epetra_RowMatrix>& Matrix) :
00041   Matrix_(Matrix),
00042   NumRows_(0),
00043   NumNonzeros_(0),
00044   MaxNumEntries_(0),
00045   MaxNumEntriesA_(0)
00046 {
00047   sprintf(Label_,"%s","Ifpack_LocalFilter");
00048 
00049 #ifdef HAVE_MPI
00050   SerialComm_ = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_SELF) );
00051 #else
00052   SerialComm_ = Teuchos::rcp( new Epetra_SerialComm );
00053 #endif
00054 
00055   // localized matrix has all the local rows of Matrix
00056   NumRows_ = Matrix->NumMyRows();
00057 
00058   // build a linear map, based on the serial communicator
00059   Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,*SerialComm_) );
00060 
00061   // NumEntries_ will contain the actual number of nonzeros
00062   // for each localized row (that is, without external nodes,
00063   // and always with the diagonal entry)
00064   NumEntries_.resize(NumRows_);
00065 
00066   // want to store the diagonal vector. FIXME: am I really useful?
00067   Diagonal_ = Teuchos::rcp( new Epetra_Vector(*Map_) );
00068   if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
00069   
00070   // store this for future access to ExtractMyRowCopy().
00071   // This is the # of nonzeros in the non-local matrix
00072   MaxNumEntriesA_ = Matrix->MaxNumEntries();
00073   // tentative value for MaxNumEntries. This is the number of
00074   // nonzeros in the local matrix
00075   MaxNumEntries_ = Matrix->MaxNumEntries();
00076 
00077   // ExtractMyRowCopy() will use these vectors
00078   Indices_.resize(MaxNumEntries_);
00079   Values_.resize(MaxNumEntries_);
00080 
00081   // now compute:
00082   // - the number of nonzero per row
00083   // - the total number of nonzeros
00084   // - the diagonal entries
00085 
00086   // compute nonzeros (total and per-row), and store the
00087   // diagonal entries (already modified)
00088   int ActualMaxNumEntries = 0;
00089 
00090   for (int i = 0 ; i < NumRows_ ; ++i) {
00091     
00092     NumEntries_[i] = 0;
00093     int Nnz, NewNnz = 0;
00094     IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntries_,Nnz,&Values_[0],&Indices_[0]));
00095 
00096     for (int j = 0 ; j < Nnz ; ++j) {
00097       if (Indices_[j] < NumRows_ ) ++NewNnz;
00098 
00099       if (Indices_[j] == i)
00100     (*Diagonal_)[i] = Values_[j];
00101     }
00102 
00103     if (NewNnz > ActualMaxNumEntries)
00104       ActualMaxNumEntries = NewNnz;
00105 
00106     NumNonzeros_ += NewNnz;
00107     NumEntries_[i] = NewNnz;
00108 
00109   }
00110  
00111   MaxNumEntries_ = ActualMaxNumEntries;
00112 }
00113 
00114 //==============================================================================
00115 int Ifpack_LocalFilter::
00116 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00117          double *Values, int * Indices) const
00118 {
00119   if ((MyRow < 0) || (MyRow >= NumRows_)) {
00120     IFPACK_CHK_ERR(-1); // range not valid
00121   }
00122 
00123   if (Length < NumEntries_[MyRow])
00124     return(-1);
00125 
00126   // always extract using the object Values_ and Indices_.
00127   // This is because I need more space than that given by
00128   // the user (for the external nodes)
00129   int Nnz;
00130   int ierr = Matrix_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz,
00131                        &Values_[0],&Indices_[0]);
00132 
00133   IFPACK_CHK_ERR(ierr);
00134 
00135   // populate the user's vectors, add diagonal if not found
00136   NumEntries = 0;
00137 
00138   for (int j = 0 ; j < Nnz ; ++j) {
00139     // only local indices
00140     if (Indices_[j] < NumRows_ ) {
00141       Indices[NumEntries] = Indices_[j];
00142       Values[NumEntries] = Values_[j];
00143       ++NumEntries;
00144     }
00145   }
00146     
00147   return(0);
00148 
00149 }
00150 
00151 //==============================================================================
00152 int Ifpack_LocalFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00153 {
00154   if (!Diagonal.Map().SameAs(*Map_))
00155     IFPACK_CHK_ERR(-1);
00156   Diagonal = *Diagonal_;
00157   return(0);
00158 }
00159 
00160 //==============================================================================
00161 int Ifpack_LocalFilter::Apply(const Epetra_MultiVector& X,
00162       Epetra_MultiVector& Y) const 
00163 {
00164 
00165   // skip expensive checks, I suppose input data are ok
00166 
00167   Y.PutScalar(0.0);
00168   int NumVectors = Y.NumVectors();
00169 
00170   double** X_ptr;
00171   double** Y_ptr;
00172   X.ExtractView(&X_ptr);
00173   Y.ExtractView(&Y_ptr);
00174 
00175   for (int i = 0 ; i < NumRows_ ; ++i) {
00176     
00177     int Nnz;
00178     int ierr = Matrix_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,&Values_[0],
00179                                          &Indices_[0]);
00180     IFPACK_CHK_ERR(ierr);
00181 
00182     for (int j = 0 ; j < Nnz ; ++j) {
00183       if (Indices_[j] < NumRows_ ) {
00184     for (int k = 0 ; k < NumVectors ; ++k)
00185       Y_ptr[k][i] += Values_[j] * X_ptr[k][Indices_[j]];
00186       }
00187     }
00188   }
00189 
00190   return(0);
00191 }
00192 
00193 //==============================================================================
00194 int Ifpack_LocalFilter::ApplyInverse(const Epetra_MultiVector& X,
00195          Epetra_MultiVector& Y) const
00196 {
00197   IFPACK_CHK_ERR(-1); // not implemented
00198 }
00199 
00200 //==============================================================================
00201 const Epetra_BlockMap& Ifpack_LocalFilter::Map() const
00202 {
00203   return(*Map_);
00204 }

Generated on Tue Oct 20 12:48:54 2009 for IFPACK by doxygen 1.4.7