Ifpack_DenseContainer.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_DenseContainer.h"
00003 #include "Epetra_RowMatrix.h"
00004 
00005 //==============================================================================
00006 int Ifpack_DenseContainer::NumRows() const
00007 {
00008   return(NumRows_);
00009 }
00010 
00011 //==============================================================================
00012 int Ifpack_DenseContainer::Initialize()
00013 {
00014   
00015   IsInitialized_ = false;
00016 
00017   IFPACK_CHK_ERR(LHS_.Reshape(NumRows_,NumVectors_));
00018   IFPACK_CHK_ERR(RHS_.Reshape(NumRows_,NumVectors_));
00019   IFPACK_CHK_ERR(ID_.Reshape(NumRows_,NumVectors_));
00020   IFPACK_CHK_ERR(Matrix_.Reshape(NumRows_,NumRows_));
00021 
00022   // zero out matrix elements
00023   for (int i = 0 ; i < NumRows_ ; ++i)
00024     for (int j = 0 ; j < NumRows_ ; ++j)
00025       Matrix_(i,j) = 0.0;
00026 
00027   // zero out vector elements
00028   for (int i = 0 ; i < NumRows_ ; ++i)
00029     for (int j = 0 ; j < NumVectors_ ; ++j) {
00030       LHS_(i,j) = 0.0;
00031       RHS_(i,j) = 0.0;
00032     }
00033 
00034   // Set to -1 ID_'s
00035   for (int i = 0 ; i < NumRows_ ; ++i)
00036     ID_(i) = -1;  
00037 
00038   if (NumRows_ != 0) {
00039     IFPACK_CHK_ERR(Solver_.SetMatrix(Matrix_));
00040     IFPACK_CHK_ERR(Solver_.SetVectors(LHS_,RHS_));
00041   }
00042 
00043   IsInitialized_ = true;
00044   return(0);
00045   
00046 }
00047 
00048 //==============================================================================
00049 double& Ifpack_DenseContainer::LHS(const int i, const int Vector)
00050 {
00051   return(LHS_.A()[Vector * NumRows_ + i]);
00052 }
00053   
00054 //==============================================================================
00055 double& Ifpack_DenseContainer::RHS(const int i, const int Vector)
00056 {
00057   return(RHS_.A()[Vector * NumRows_ + i]);
00058 }
00059 
00060 //==============================================================================
00061 int Ifpack_DenseContainer::
00062 SetMatrixElement(const int row, const int col, const double value)
00063 {
00064   if (IsInitialized() == false)
00065     IFPACK_CHK_ERR(Initialize());
00066 
00067   if ((row < 0) || (row >= NumRows())) {
00068     IFPACK_CHK_ERR(-2); // not in range
00069   }
00070 
00071   if ((col < 0) || (col >= NumRows())) {
00072     IFPACK_CHK_ERR(-2); // not in range
00073   }
00074 
00075   Matrix_(row, col) = value;
00076 
00077   return(0);
00078 
00079 }
00080 
00081 //==============================================================================
00082 int Ifpack_DenseContainer::ApplyInverse()
00083 {
00084 
00085   if (!IsComputed()) {
00086     IFPACK_CHK_ERR(-1);
00087   }
00088   
00089   if (NumRows_ != 0)
00090     IFPACK_CHK_ERR(Solver_.Solve());
00091 
00092   ApplyInverseFlops_ += 2.0 * NumVectors_ * NumRows_ * NumRows_;
00093   return(0);
00094 }
00095 
00096 //==============================================================================
00097 int& Ifpack_DenseContainer::ID(const int i)
00098 {
00099   return(ID_[i]);
00100 }
00101 
00102 //==============================================================================
00103 // FIXME: optimize performances of this guy...
00104 int Ifpack_DenseContainer::Extract(const Epetra_RowMatrix& Matrix)
00105 {
00106 
00107   for (int j = 0 ; j < NumRows_ ; ++j) {
00108     // be sure that the user has set all the ID's
00109     if (ID(j) == -1)
00110       IFPACK_CHK_ERR(-2);
00111     // be sure that all are local indices
00112     if (ID(j) > Matrix.NumMyRows())
00113       IFPACK_CHK_ERR(-2);
00114   }
00115 
00116   // allocate storage to extract matrix rows.
00117   int Length = Matrix.MaxNumEntries();
00118   vector<double> Values;
00119   Values.resize(Length);
00120   vector<int> Indices;
00121   Indices.resize(Length);
00122 
00123   for (int j = 0 ; j < NumRows_ ; ++j) {
00124 
00125     int LRID = ID(j);
00126 
00127     int NumEntries;
00128 
00129     int ierr = 
00130       Matrix.ExtractMyRowCopy(LRID, Length, NumEntries, 
00131                   &Values[0], &Indices[0]);
00132     IFPACK_CHK_ERR(ierr);
00133 
00134     for (int k = 0 ; k < NumEntries ; ++k) {
00135 
00136       int LCID = Indices[k];
00137 
00138       // skip off-processor elements
00139       if (LCID >= Matrix.NumMyRows()) 
00140     continue;
00141 
00142       // for local column IDs, look for each ID in the list
00143       // of columns hosted by this object
00144       // FIXME: use STL
00145       int jj = -1;
00146       for (int kk = 0 ; kk < NumRows_ ; ++kk)
00147     if (ID(kk) == LCID)
00148       jj = kk;
00149 
00150       if (jj != -1)
00151     SetMatrixElement(j,jj,Values[k]);
00152 
00153     }
00154   }
00155 
00156   return(0);
00157 }
00158 
00159 //==============================================================================
00160 int Ifpack_DenseContainer::Compute(const Epetra_RowMatrix& Matrix)
00161 {
00162   IsComputed_ = false;
00163   if (IsInitialized() == false) {
00164     IFPACK_CHK_ERR(Initialize());
00165   }
00166 
00167   if (KeepNonFactoredMatrix_)
00168     NonFactoredMatrix_ = Matrix_;
00169 
00170   // extract local rows and columns
00171   IFPACK_CHK_ERR(Extract(Matrix));
00172 
00173   if (KeepNonFactoredMatrix_)
00174     NonFactoredMatrix_ = Matrix_;
00175 
00176   // factorize the matrix using LAPACK
00177   if (NumRows_ != 0)
00178     IFPACK_CHK_ERR(Solver_.Factor());
00179 
00180   Label_ = "Ifpack_DenseContainer";
00181 
00182   // not sure of count
00183   ComputeFlops_ += 4.0 * NumRows_ * NumRows_ * NumRows_ / 3;
00184   IsComputed_ = true;
00185 
00186   return(0);
00187 }
00188 
00189 //==============================================================================
00190 int Ifpack_DenseContainer::Apply()
00191 {
00192   if (IsComputed() == false)
00193     IFPACK_CHK_ERR(-3);
00194 
00195   if (KeepNonFactoredMatrix_) {
00196     IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,NonFactoredMatrix_,LHS_,0.0));
00197   }
00198   else
00199     IFPACK_CHK_ERR(RHS_.Multiply('N','N', 1.0,Matrix_,LHS_,0.0));
00200 
00201   ApplyFlops_ += 2 * NumRows_ * NumRows_;
00202   return(0);
00203 }
00204 
00205 //==============================================================================
00206 ostream& Ifpack_DenseContainer::Print(ostream & os) const
00207 {
00208     os << "================================================================================" << endl;
00209   os << "Ifpack_DenseContainer" << endl;
00210   os << "Number of rows          = " << NumRows() << endl;
00211   os << "Number of vectors       = " << NumVectors() << endl;
00212   os << "IsInitialized()         = " << IsInitialized() << endl;
00213   os << "IsComputed()            = " << IsComputed() << endl;
00214   os << "Flops in Compute()      = " << ComputeFlops() << endl; 
00215   os << "Flops in ApplyInverse() = " << ApplyInverseFlops() << endl; 
00216   os << "================================================================================" << endl;
00217   os << endl;
00218 
00219   return(os);
00220 }

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1