IFPACK Development
Ifpack_SingletonFilter.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 // 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 
00043 #include "Ifpack_ConfigDefs.h"
00044 #include "Ifpack_SingletonFilter.h"
00045 #include "Epetra_ConfigDefs.h"
00046 #include "Epetra_RowMatrix.h"
00047 #include "Epetra_Comm.h"
00048 #include "Epetra_Map.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_Vector.h"
00051 #include "Epetra_Import.h"
00052 
00053 //==============================================================================
00054 Ifpack_SingletonFilter::Ifpack_SingletonFilter(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix) :
00055   A_(Matrix),
00056   NumSingletons_(0),
00057   NumRows_(0),
00058   NumRowsA_(0),
00059   MaxNumEntries_(0),
00060   MaxNumEntriesA_(0),
00061   NumNonzeros_(0)
00062 {
00063   // use this filter only on serial matrices
00064   if (A_->Comm().NumProc() != 1) {
00065     cerr << "Ifpack_SingletonFilter can be used with Comm().NumProc() == 1" << endl;
00066     cerr << "only. This class is a tool for Ifpack_AdditiveSchwarz," << endl;
00067     cerr << "and it is not meant to be used otherwise." << endl;
00068     exit(EXIT_FAILURE);
00069   }
00070   
00071   if ((A_->NumMyRows() != A_->NumGlobalRows64()) ||
00072      (A_->NumMyRows() != A_->NumMyCols()))
00073     IFPACK_CHK_ERRV(-1);
00074   
00075   NumRowsA_ = (A_->NumMyRows());
00076   MaxNumEntriesA_ = A_->MaxNumEntries();
00077 
00078   Indices_.resize(MaxNumEntriesA_);
00079   Values_.resize(MaxNumEntriesA_);
00080   Reorder_.resize(A_->NumMyRows());
00081 
00082   for (int i = 0 ; i < NumRowsA_ ; ++i)
00083     Reorder_[i] = -1;
00084 
00085   // first check how may singletons I do have
00086   for (int i = 0 ; i < NumRowsA_ ; ++i) {
00087     int Nnz;
00088     IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00089                     &Values_[0], &Indices_[0]));
00090     if (Nnz != 1) {
00091       Reorder_[i] = NumRows_++;
00092     }
00093     else {
00094       NumSingletons_++;
00095     }
00096   }
00097 
00098   InvReorder_.resize(NumRows_);
00099   for (int i = 0 ; i < NumRowsA_ ; ++i) {
00100     if (Reorder_[i] < 0)
00101       continue;
00102     InvReorder_[Reorder_[i]] = i;
00103   }
00104  
00105   NumEntries_.resize(NumRows_);
00106   SingletonIndex_.resize(NumSingletons_);
00107 
00108   // now compute the nonzeros per row
00109   int count = 0;
00110   for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
00111 
00112     int Nnz;
00113     IFPACK_CHK_ERRV(A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00114                       &Values_[0], &Indices_[0]));
00115     int ii = Reorder_[i];
00116     if (ii >= 0) {
00117       assert (Nnz != 1);
00118 
00119       NumEntries_[ii] = Nnz;
00120       NumNonzeros_ += Nnz;
00121       if (Nnz > MaxNumEntries_)
00122     MaxNumEntries_ = Nnz;
00123     }
00124     else {
00125       SingletonIndex_[count] = i;
00126       count++;
00127     }
00128   }
00129 
00130 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
00131   Map_ = Teuchos::rcp( new Epetra_Map(NumRows_,0,Comm()) );
00132 #endif
00133 
00134   // and finish up with the diagonal entry
00135   Diagonal_ = Teuchos::rcp( new Epetra_Vector(*Map_) );
00136 
00137   Epetra_Vector Diagonal(A_->Map());
00138   A_->ExtractDiagonalCopy(Diagonal);
00139   for (int i = 0 ; i < NumRows_ ; ++i) {
00140     int ii = InvReorder_[i];
00141     (*Diagonal_)[i] = Diagonal[ii];
00142   }
00143   
00144 }
00145 
00146 //==============================================================================
00147 int Ifpack_SingletonFilter::
00148 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00149          double *Values, int * Indices) const
00150 {
00151   int Nnz;
00152 
00153   if (Length < NumEntries_[MyRow])
00154     IFPACK_CHK_ERR(-1);
00155 
00156   int Row = InvReorder_[MyRow];
00157   IFPACK_CHK_ERR(A_->ExtractMyRowCopy(Row,MaxNumEntriesA_,Nnz,
00158                      &Values_[0],&Indices_[0]));
00159   NumEntries = 0;
00160   for (int i = 0 ; i < Nnz ; ++i) {
00161     int ii = Reorder_[Indices_[i]];
00162     if ( ii >= 0) {
00163       Indices[NumEntries] = ii;
00164       Values[NumEntries] = Values_[i];
00165       NumEntries++;
00166     }
00167   }
00168   return(0);   
00169 }
00170 
00171 //==============================================================================
00172 int Ifpack_SingletonFilter::
00173 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00174 {
00175   Diagonal = *Diagonal_;
00176   return(0);
00177 }
00178 
00179 //==============================================================================
00180 int Ifpack_SingletonFilter::
00181 Multiply(bool TransA, const Epetra_MultiVector& X, 
00182      Epetra_MultiVector& Y) const
00183 {
00184 
00185   int NumVectors = X.NumVectors();
00186   if (NumVectors != Y.NumVectors())
00187     IFPACK_CHK_ERR(-1);
00188 
00189   Y.PutScalar(0.0);
00190 
00191   std::vector<int> Indices(MaxNumEntries_);
00192   std::vector<double> Values(MaxNumEntries_);
00193 
00194   // cycle over all rows of the original matrix
00195   for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
00196 
00197     if (Reorder_[i] < 0)
00198       continue; // skip singleton rows
00199     
00200     int Nnz;
00201     A_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,
00202              &Values[0], &Indices[0]);
00203     if (!TransA) {
00204       // no transpose first
00205       for (int j = 0 ; j < NumVectors ; ++j) {
00206     for (int k = 0 ; k < Nnz ; ++k) {
00207       if (Reorder_[i] >= 0)
00208         Y[j][i] += Values[k] * X[j][Reorder_[Indices[k]]];
00209     }
00210       }
00211     }
00212     else {
00213       // transpose here
00214       for (int j = 0 ; j < NumVectors ; ++j) {
00215     for (int k = 0 ; k < Nnz ; ++k) {
00216       if (Reorder_[i] >= 0)
00217         Y[j][Reorder_[Indices[k]]] += Values[k] * X[j][i];
00218     }
00219       }
00220     }
00221   }
00222 
00223   return(0);
00224 }
00225 
00226 //==============================================================================
00227 int Ifpack_SingletonFilter::
00228 Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00229       const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00230 {
00231   return(-1);
00232 }
00233 
00234 //==============================================================================
00235 int Ifpack_SingletonFilter::
00236 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00237 {
00238   IFPACK_CHK_ERR(Multiply(false,X,Y));
00239   return(0);
00240 }
00241 
00242 //==============================================================================
00243 int Ifpack_SingletonFilter::
00244 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00245 {
00246   return(-1); // NOT IMPLEMENTED AT THIS STAGE
00247 }
00248 
00249 //==============================================================================
00250 int Ifpack_SingletonFilter::
00251 SolveSingletons(const Epetra_MultiVector& RHS, 
00252         Epetra_MultiVector& LHS)
00253 {
00254   for (int i = 0 ; i < NumSingletons_ ; ++i) {
00255     int ii = SingletonIndex_[i];
00256     // get the diagonal value for the singleton
00257     int Nnz;
00258     A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
00259             &Values_[0], &Indices_[0]);
00260     for (int j = 0 ; j < Nnz ; ++j) {
00261       if (Indices_[j] == ii) {
00262     for (int k = 0 ; k < LHS.NumVectors() ; ++k)
00263       LHS[k][ii] = RHS[k][ii] / Values_[j];
00264       }
00265     }
00266   }
00267 
00268   return(0);
00269 }
00270 
00271 //==============================================================================
00272 int Ifpack_SingletonFilter::
00273 CreateReducedRHS(const Epetra_MultiVector& LHS,
00274          const Epetra_MultiVector& RHS, 
00275          Epetra_MultiVector& ReducedRHS)
00276 {
00277   int NumVectors = LHS.NumVectors();
00278 
00279   for (int i = 0 ; i < NumRows_ ; ++i)
00280     for (int k = 0 ; k < NumVectors ; ++k)
00281       ReducedRHS[k][i] = RHS[k][InvReorder_[i]];
00282 
00283   for (int i = 0 ; i < NumRows_ ; ++i) {
00284     int ii = InvReorder_[i];
00285     int Nnz;
00286     IFPACK_CHK_ERR(A_->ExtractMyRowCopy(ii,MaxNumEntriesA_,Nnz,
00287                     &Values_[0], &Indices_[0]));
00288 
00289     for (int j = 0 ; j < Nnz ; ++j) {
00290       if (Reorder_[Indices_[j]] == -1) {
00291     for (int k = 0 ; k < NumVectors ; ++k)
00292       ReducedRHS[k][i] -= Values_[j] * LHS[k][Indices_[j]];
00293       }
00294     }
00295   }
00296   return(0);
00297 }
00298 
00299 //==============================================================================
00300 int Ifpack_SingletonFilter::
00301 UpdateLHS(const Epetra_MultiVector& ReducedLHS,
00302       Epetra_MultiVector& LHS)
00303 {
00304   for (int i = 0 ; i < NumRows_ ; ++i)
00305     for (int k = 0 ; k < LHS.NumVectors() ; ++k)
00306       LHS[k][InvReorder_[i]] = ReducedLHS[k][i];
00307 
00308   return(0);
00309 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends