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 // 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 #ifdef IFPACK_FLOPCOUNTERS
00122   ApplyInverseFlops_ += 2.0 * NumVectors_ * NumRows_ * NumRows_;
00123 #endif
00124   return(0);
00125 }
00126 
00127 //==============================================================================
00128 int& Ifpack_DenseContainer::ID(const int i)
00129 {
00130   return(ID_[i]);
00131 }
00132 
00133 //==============================================================================
00134 // FIXME: optimize performances of this guy...
00135 int Ifpack_DenseContainer::Extract(const Epetra_RowMatrix& Matrix_in)
00136 {
00137 
00138   for (int j = 0 ; j < NumRows_ ; ++j) {
00139     // be sure that the user has set all the ID's
00140     if (ID(j) == -1)
00141       IFPACK_CHK_ERR(-2);
00142     // be sure that all are local indices
00143     if (ID(j) > Matrix_in.NumMyRows())
00144       IFPACK_CHK_ERR(-2);
00145   }
00146 
00147   // allocate storage to extract matrix rows.
00148   int Length = Matrix_in.MaxNumEntries();
00149   vector<double> Values;
00150   Values.resize(Length);
00151   vector<int> Indices;
00152   Indices.resize(Length);
00153 
00154   for (int j = 0 ; j < NumRows_ ; ++j) {
00155 
00156     int LRID = ID(j);
00157 
00158     int NumEntries;
00159 
00160     int ierr = 
00161       Matrix_in.ExtractMyRowCopy(LRID, Length, NumEntries, 
00162                   &Values[0], &Indices[0]);
00163     IFPACK_CHK_ERR(ierr);
00164 
00165     for (int k = 0 ; k < NumEntries ; ++k) {
00166 
00167       int LCID = Indices[k];
00168 
00169       // skip off-processor elements
00170       if (LCID >= Matrix_in.NumMyRows()) 
00171     continue;
00172 
00173       // for local column IDs, look for each ID in the list
00174       // of columns hosted by this object
00175       // FIXME: use STL
00176       int jj = -1;
00177       for (int kk = 0 ; kk < NumRows_ ; ++kk)
00178     if (ID(kk) == LCID)
00179       jj = kk;
00180 
00181       if (jj != -1)
00182     SetMatrixElement(j,jj,Values[k]);
00183 
00184     }
00185   }
00186 
00187   return(0);
00188 }
00189 
00190 //==============================================================================
00191 int Ifpack_DenseContainer::Compute(const Epetra_RowMatrix& Matrix_in)
00192 {
00193   IsComputed_ = false;
00194   if (IsInitialized() == false) {
00195     IFPACK_CHK_ERR(Initialize());
00196   }
00197 
00198   if (KeepNonFactoredMatrix_)
00199     NonFactoredMatrix_ = Matrix_;
00200 
00201   // extract local rows and columns
00202   IFPACK_CHK_ERR(Extract(Matrix_in));
00203 
00204   if (KeepNonFactoredMatrix_)
00205     NonFactoredMatrix_ = Matrix_;
00206 
00207   // factorize the matrix using LAPACK
00208   if (NumRows_ != 0)
00209     IFPACK_CHK_ERR(Solver_.Factor());
00210 
00211   Label_ = "Ifpack_DenseContainer";
00212 
00213   // not sure of count
00214 #ifdef IFPACK_FLOPCOUNTERS
00215   ComputeFlops_ += 4.0 * NumRows_ * NumRows_ * NumRows_ / 3;
00216 #endif
00217   IsComputed_ = true;
00218 
00219   return(0);
00220 }
00221 
00222 //==============================================================================
00223 int Ifpack_DenseContainer::Apply()
00224 {
00225   if (IsComputed() == false)
00226     IFPACK_CHK_ERR(-3);
00227 
00228   if (KeepNonFactoredMatrix_) {
00229     IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,NonFactoredMatrix_,LHS_,0.0));
00230   }
00231   else
00232     IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,Matrix_,LHS_,0.0));
00233 
00234 #ifdef IFPACK_FLOPCOUNTERS
00235   ApplyFlops_ += 2 * NumRows_ * NumRows_;
00236 #endif
00237   return(0);
00238 }
00239 
00240 //==============================================================================
00241 ostream& Ifpack_DenseContainer::Print(ostream & os) const
00242 {
00243     os << "================================================================================" << endl;
00244   os << "Ifpack_DenseContainer" << endl;
00245   os << "Number of rows          = " << NumRows() << endl;
00246   os << "Number of vectors       = " << NumVectors() << endl;
00247   os << "IsInitialized()         = " << IsInitialized() << endl;
00248   os << "IsComputed()            = " << IsComputed() << endl;
00249 #ifdef IFPACK_FLOPCOUNTERS
00250   os << "Flops in Compute()      = " << ComputeFlops() << endl; 
00251   os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl; 
00252 #endif
00253   os << "================================================================================" << endl;
00254   os << endl;
00255 
00256   return(os);
00257 }
 All Classes Files Functions Variables Enumerations Friends