Ifpack_DenseContainer.cpp

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