IFPACK Development
Ifpack_DenseContainer.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_DenseContainer.h"
00045 #include "Epetra_RowMatrix.h"
00046 
00047 //==============================================================================
00048 int Ifpack_DenseContainer::NumRows() const
00049 {
00050   return(NumRows_);
00051 }
00052 
00053 //==============================================================================
00054 int Ifpack_DenseContainer::Initialize()
00055 {
00056   
00057   IsInitialized_ = false;
00058 
00059   IFPACK_CHK_ERR(LHS_.Reshape(NumRows_,NumVectors_));
00060   IFPACK_CHK_ERR(RHS_.Reshape(NumRows_,NumVectors_));
00061   IFPACK_CHK_ERR(ID_.Reshape(NumRows_,NumVectors_));
00062   IFPACK_CHK_ERR(Matrix_.Reshape(NumRows_,NumRows_));
00063 
00064   // zero out matrix elements
00065   for (int i = 0 ; i < NumRows_ ; ++i)
00066     for (int j = 0 ; j < NumRows_ ; ++j)
00067       Matrix_(i,j) = 0.0;
00068 
00069   // zero out vector elements
00070   for (int i = 0 ; i < NumRows_ ; ++i)
00071     for (int j = 0 ; j < NumVectors_ ; ++j) {
00072       LHS_(i,j) = 0.0;
00073       RHS_(i,j) = 0.0;
00074     }
00075 
00076   // Set to -1 ID_'s
00077   for (int i = 0 ; i < NumRows_ ; ++i)
00078     ID_(i) = -1;  
00079 
00080   if (NumRows_ != 0) {
00081     IFPACK_CHK_ERR(Solver_.SetMatrix(Matrix_));
00082     IFPACK_CHK_ERR(Solver_.SetVectors(LHS_,RHS_));
00083   }
00084 
00085   IsInitialized_ = true;
00086   return(0);
00087   
00088 }
00089 
00090 //==============================================================================
00091 double& Ifpack_DenseContainer::LHS(const int i, const int Vector)
00092 {
00093   return(LHS_.A()[Vector * NumRows_ + i]);
00094 }
00095   
00096 //==============================================================================
00097 double& Ifpack_DenseContainer::RHS(const int i, const int Vector)
00098 {
00099   return(RHS_.A()[Vector * NumRows_ + i]);
00100 }
00101 
00102 //==============================================================================
00103 int Ifpack_DenseContainer::
00104 SetMatrixElement(const int row, const int col, const double value)
00105 {
00106   if (IsInitialized() == false)
00107     IFPACK_CHK_ERR(Initialize());
00108 
00109   if ((row < 0) || (row >= NumRows())) {
00110     IFPACK_CHK_ERR(-2); // not in range
00111   }
00112 
00113   if ((col < 0) || (col >= NumRows())) {
00114     IFPACK_CHK_ERR(-2); // not in range
00115   }
00116 
00117   Matrix_(row, col) = value;
00118 
00119   return(0);
00120 
00121 }
00122 
00123 //==============================================================================
00124 int Ifpack_DenseContainer::ApplyInverse()
00125 {
00126 
00127   if (!IsComputed()) {
00128     IFPACK_CHK_ERR(-1);
00129   }
00130   
00131   if (NumRows_ != 0)
00132     IFPACK_CHK_ERR(Solver_.Solve());
00133 
00134 #ifdef IFPACK_FLOPCOUNTERS
00135   ApplyInverseFlops_ += 2.0 * NumVectors_ * NumRows_ * NumRows_;
00136 #endif
00137   return(0);
00138 }
00139 
00140 //==============================================================================
00141 int& Ifpack_DenseContainer::ID(const int i)
00142 {
00143   return(ID_[i]);
00144 }
00145 
00146 //==============================================================================
00147 // FIXME: optimize performances of this guy...
00148 int Ifpack_DenseContainer::Extract(const Epetra_RowMatrix& Matrix_in)
00149 {
00150 
00151   for (int j = 0 ; j < NumRows_ ; ++j) {
00152     // be sure that the user has set all the ID's
00153     if (ID(j) == -1)
00154       IFPACK_CHK_ERR(-2);
00155     // be sure that all are local indices
00156     if (ID(j) > Matrix_in.NumMyRows())
00157       IFPACK_CHK_ERR(-2);
00158   }
00159 
00160   // allocate storage to extract matrix rows.
00161   int Length = Matrix_in.MaxNumEntries();
00162   std::vector<double> Values;
00163   Values.resize(Length);
00164   std::vector<int> Indices;
00165   Indices.resize(Length);
00166 
00167   for (int j = 0 ; j < NumRows_ ; ++j) {
00168 
00169     int LRID = ID(j);
00170 
00171     int NumEntries;
00172 
00173     int ierr = 
00174       Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries, 
00175                   &Values[0], &Indices[0]);
00176     IFPACK_CHK_ERR(ierr);
00177 
00178     for (int k = 0 ; k < NumEntries ; ++k) {
00179 
00180       int LCID = Indices[k];
00181 
00182       // skip off-processor elements
00183       if (LCID >= Matrix_in.NumMyRows()) 
00184     continue;
00185 
00186       // for local column IDs, look for each ID in the list
00187       // of columns hosted by this object
00188       // FIXME: use STL
00189       int jj = -1;
00190       for (int kk = 0 ; kk < NumRows_ ; ++kk)
00191     if (ID(kk) == LCID)
00192       jj = kk;
00193 
00194       if (jj != -1)
00195     SetMatrixElement(j,jj,Values[k]);
00196 
00197     }
00198   }
00199 
00200   return(0);
00201 }
00202 
00203 //==============================================================================
00204 int Ifpack_DenseContainer::Compute(const Epetra_RowMatrix& Matrix_in)
00205 {
00206   IsComputed_ = false;
00207   if (IsInitialized() == false) {
00208     IFPACK_CHK_ERR(Initialize());
00209   }
00210 
00211   if (KeepNonFactoredMatrix_)
00212     NonFactoredMatrix_ = Matrix_;
00213 
00214   // extract local rows and columns
00215   IFPACK_CHK_ERR(Extract(Matrix_in));
00216 
00217   if (KeepNonFactoredMatrix_)
00218     NonFactoredMatrix_ = Matrix_;
00219 
00220   // factorize the matrix using LAPACK
00221   if (NumRows_ != 0)
00222     IFPACK_CHK_ERR(Solver_.Factor());
00223 
00224   Label_ = "Ifpack_DenseContainer";
00225 
00226   // not sure of count
00227 #ifdef IFPACK_FLOPCOUNTERS
00228   ComputeFlops_ += 4.0 * NumRows_ * NumRows_ * NumRows_ / 3;
00229 #endif
00230   IsComputed_ = true;
00231 
00232   return(0);
00233 }
00234 
00235 //==============================================================================
00236 int Ifpack_DenseContainer::Apply()
00237 {
00238   if (IsComputed() == false)
00239     IFPACK_CHK_ERR(-3);
00240 
00241   if (KeepNonFactoredMatrix_) {
00242     IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,NonFactoredMatrix_,LHS_,0.0));
00243   }
00244   else
00245     IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,Matrix_,LHS_,0.0));
00246 
00247 #ifdef IFPACK_FLOPCOUNTERS
00248   ApplyFlops_ += 2 * NumRows_ * NumRows_;
00249 #endif
00250   return(0);
00251 }
00252 
00253 //==============================================================================
00254 ostream& Ifpack_DenseContainer::Print(ostream & os) const
00255 {
00256     os << "================================================================================" << endl;
00257   os << "Ifpack_DenseContainer" << endl;
00258   os << "Number of rows          = " << NumRows() << endl;
00259   os << "Number of vectors       = " << NumVectors() << endl;
00260   os << "IsInitialized()         = " << IsInitialized() << endl;
00261   os << "IsComputed()            = " << IsComputed() << endl;
00262 #ifdef IFPACK_FLOPCOUNTERS
00263   os << "Flops in Compute()      = " << ComputeFlops() << endl; 
00264   os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl; 
00265 #endif
00266   os << "================================================================================" << endl;
00267   os << endl;
00268 
00269   return(os);
00270 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends