Ifpack_CrsIct.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_CrsIct.h"
00031 #include "Epetra_Comm.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_CrsGraph.h"
00034 #include "Epetra_CrsMatrix.h"
00035 #include "Epetra_Vector.h"
00036 #include "Epetra_MultiVector.h"
00037 #include "Epetra_Util.h"
00038 #include "icrout_cholesky_mex.h"
00039 
00040 #ifdef HAVE_IFPACK_TEUCHOS
00041 #include <Teuchos_ParameterList.hpp>
00042 #include <ifp_parameters.h>
00043 #endif
00044 
00045 //==============================================================================
00046 Ifpack_CrsIct::Ifpack_CrsIct(const Epetra_CrsMatrix & A, double Droptol, int Lfil) 
00047   : A_(A),
00048     Comm_(A.Comm()),
00049     Allocated_(false),
00050     ValuesInitialized_(false),
00051     Factored_(false),
00052     Condest_(-1.0),
00053     Athresh_(0.0),
00054     Rthresh_(1.0),
00055     Droptol_(Droptol),
00056     Lfil_(Lfil),
00057     OverlapX_(0),
00058     OverlapY_(0),
00059     LevelOverlap_(0),
00060     OverlapMode_(Zero),
00061     Aict_(0),
00062     Lict_(0),
00063     Ldiag_(0)
00064 {
00065   Allocate();
00066 }
00067 
00068 //==============================================================================
00069 Ifpack_CrsIct::Ifpack_CrsIct(const Ifpack_CrsIct & FactoredMatrix) 
00070   : A_(FactoredMatrix.A_),
00071     Comm_(FactoredMatrix.Comm_),
00072     Allocated_(FactoredMatrix.Allocated_),
00073     ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
00074     Factored_(FactoredMatrix.Factored_),
00075     Condest_(FactoredMatrix.Condest_),
00076     Athresh_(FactoredMatrix.Athresh_),
00077     Rthresh_(FactoredMatrix.Rthresh_),
00078     Droptol_(FactoredMatrix.Droptol_),
00079     Lfil_(FactoredMatrix.Lfil_),
00080     OverlapX_(0),
00081     OverlapY_(0),
00082     LevelOverlap_(FactoredMatrix.LevelOverlap_),
00083     OverlapMode_(FactoredMatrix.OverlapMode_),
00084     Aict_(0),
00085     Lict_(0),
00086     Ldiag_(0)
00087   
00088 {
00089   U_ = new Epetra_CrsMatrix(FactoredMatrix.U());
00090   D_ = new Epetra_Vector(FactoredMatrix.D());
00091   
00092 }
00093 
00094 //==============================================================================
00095 int Ifpack_CrsIct::Allocate() {
00096 
00097   // Allocate Epetra_CrsMatrix using ILUK graphs
00098   if (LevelOverlap_==0) {
00099     U_ = new Epetra_CrsMatrix(Copy, A_.RowMatrixRowMap(), A_.RowMatrixRowMap(), 0);
00100     D_ = new Epetra_Vector(A_.RowMatrixRowMap());
00101   }
00102   else {
00103     EPETRA_CHK_ERR(-1); // LevelOverlap > 0 not implemented yet
00104     //    U_ = new Epetra_CrsMatrix(Copy, OverlapRowMap());
00105     //    D_ = new Epetra_Vector(OverlapRowMap());
00106   }
00107     
00108   
00109   
00110     SetAllocated(true);
00111     return(0);
00112 }
00113 //==============================================================================
00114 Ifpack_CrsIct::~Ifpack_CrsIct(){
00115 
00116 
00117   delete U_;
00118   delete D_; // Diagonal is stored separately.  We store the inverse.
00119 
00120   if (OverlapX_!=0) delete OverlapX_;
00121   if (OverlapY_!=0) delete OverlapY_;
00122 
00123   if (Lict_!=0) {
00124     Matrix * Lict = (Matrix *) Lict_;
00125     free(Lict->ptr);
00126     free(Lict->col);
00127     free(Lict->val);
00128     delete Lict;
00129   }
00130   if (Aict_!=0) {
00131     Matrix * Aict = (Matrix *) Aict_;
00132     delete Aict;
00133   }
00134   if (Ldiag_!=0) free(Ldiag_);
00135 
00136   ValuesInitialized_ = false;
00137   Factored_ = false;
00138   Allocated_ = false;
00139 }
00140 
00141 #ifdef HAVE_IFPACK_TEUCHOS
00142 //==========================================================================
00143 int Ifpack_CrsIct::SetParameters(const Teuchos::ParameterList& parameterlist,
00144                  bool cerr_warning_if_unused)
00145 {
00146   Ifpack::param_struct params;
00147   params.int_params[Ifpack::level_fill-FIRST_INT_PARAM] = Lfil_;
00148   params.double_params[Ifpack::absolute_threshold] = Athresh_;
00149   params.double_params[Ifpack::relative_threshold] = Rthresh_;
00150   params.double_params[Ifpack::drop_tolerance] = Droptol_;
00151   params.overlap_mode = OverlapMode_;
00152 
00153   Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
00154 
00155   Lfil_ = params.int_params[Ifpack::level_fill-FIRST_INT_PARAM];
00156   Athresh_ = params.double_params[Ifpack::absolute_threshold];
00157   Rthresh_ = params.double_params[Ifpack::relative_threshold];
00158   Droptol_ = params.double_params[Ifpack::drop_tolerance];
00159   OverlapMode_ = params.overlap_mode;
00160 
00161   return(0);
00162 }
00163 #endif
00164 
00165 //==========================================================================
00166 int Ifpack_CrsIct::InitValues(const Epetra_CrsMatrix & A) {
00167 
00168   int ierr = 0;
00169   int i, j;
00170   int * InI=0, * UI = 0;
00171   double * InV=0, * UV = 0;
00172   int NumIn, NumL, NumU;
00173   bool DiagFound;
00174   int NumNonzeroDiags = 0;
00175 
00176   Epetra_CrsMatrix * OverlapA = (Epetra_CrsMatrix *) &A_;
00177 
00178   if (LevelOverlap_>0) {
00179     EPETRA_CHK_ERR(-1); // Not implemented yet
00180     //OverlapA = new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph());
00181     //EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));
00182     //EPETRA_CHK_ERR(OverlapA->FillComplete());
00183   }
00184   // Get Maximun Row length
00185   int MaxNumEntries = OverlapA->MaxNumEntries();
00186 
00187   InI = new int[MaxNumEntries]; // Allocate temp space
00188   UI = new int[MaxNumEntries];
00189   InV = new double[MaxNumEntries];
00190   UV = new double[MaxNumEntries];
00191 
00192   double *DV;
00193   ierr = D_->ExtractView(&DV); // Get view of diagonal
00194     
00195 
00196   // First we copy the user's matrix into diagonal vector and U, regardless of fill level
00197 
00198   int NumRows = OverlapA->NumMyRows();
00199 
00200   for (i=0; i< NumRows; i++) {
00201 
00202     OverlapA->ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI); // Get Values and Indices
00203     
00204     // Split into L and U (we don't assume that indices are ordered).
00205     
00206     NumL = 0; 
00207     NumU = 0; 
00208     DiagFound = false;
00209     
00210     for (j=0; j< NumIn; j++) {
00211       int k = InI[j];
00212 
00213       if (k==i) {
00214     DiagFound = true;
00215     DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
00216       }
00217 
00218       else if (k < 0) return(-1); // Out of range
00219       else if (i<k && k<NumRows) {
00220     UI[NumU] = k;
00221     UV[NumU] = InV[j];
00222     NumU++;
00223       }
00224     }
00225     
00226     // Check in things for this row of L and U
00227 
00228     if (DiagFound) NumNonzeroDiags++;
00229     if (NumU) U_->InsertMyValues(i, NumU, UV, UI);
00230     
00231   }
00232 
00233   delete [] UI;
00234   delete [] UV;
00235   delete [] InI;
00236   delete [] InV;
00237 
00238   if (LevelOverlap_>0 && U().DistributedGlobal()) delete OverlapA;
00239 
00240 
00241   U_->FillComplete(A_.OperatorDomainMap(), A_.OperatorRangeMap());
00242   SetValuesInitialized(true);
00243   SetFactored(false);
00244 
00245   int ierr1 = 0;
00246   if (NumNonzeroDiags<U_->NumMyRows()) ierr1 = 1;
00247   A_.Comm().MaxAll(&ierr1, &ierr, 1);
00248   EPETRA_CHK_ERR(ierr);
00249   return(0);
00250 }
00251 
00252 //==========================================================================
00253 int Ifpack_CrsIct::Factor() {
00254 
00255   // if (!Allocated()) return(-1); // This test is not needed at this time.  All constructors allocate.
00256   if (!ValuesInitialized_) EPETRA_CHK_ERR(-2); // Must have values initialized.
00257   if (Factored()) EPETRA_CHK_ERR(-3); // Can't have already computed factors.
00258 
00259   SetValuesInitialized(false);
00260 
00261   int i;
00262 
00263   int m, n, nz, Nrhs, ldrhs, ldlhs;
00264   int * ptr=0, * ind;
00265   double * val, * rhs, * lhs;
00266 
00267   int ierr = Epetra_Util_ExtractHbData(U_, 0, 0, m, n, nz, ptr, ind,
00268                 val, Nrhs, rhs, ldrhs, lhs, ldlhs);
00269   if (ierr<0) EPETRA_CHK_ERR(ierr);
00270 
00271   Matrix * Aict;
00272   if (Aict_==0) {
00273     Aict = new Matrix;
00274     Aict_ = (void *) Aict;
00275   }
00276   else Aict = (Matrix *) Aict_;
00277   Matrix * Lict;
00278   if (Lict_==0) {
00279     Lict = new Matrix;
00280     Lict_ = (void *) Lict;
00281   }
00282   else Lict = (Matrix *) Lict_;
00283   Aict->val = val;
00284   Aict->col = ind;
00285   Aict->ptr = ptr;
00286   double *DV;
00287   EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
00288     
00289   crout_ict(m, Aict, DV, Droptol_, Lfil_, Lict, &Ldiag_);
00290 
00291   // Get rid of unnecessary data
00292   delete [] ptr;
00293   delete U_;
00294   delete D_;
00295 
00296   // Create Epetra View of L from crout_ict
00297 
00298   if (LevelOverlap_==0) {
00299     U_ = new Epetra_CrsMatrix(View, A_.RowMatrixRowMap(), A_.RowMatrixRowMap(),0);
00300     D_ = new Epetra_Vector(View, A_.RowMatrixRowMap(), Ldiag_);
00301   }
00302   else {
00303     EPETRA_CHK_ERR(-1); // LevelOverlap > 0 not implemented yet
00304     //    U_ = new Epetra_CrsMatrix(Copy, OverlapRowMap());
00305     //    D_ = new Epetra_Vector(OverlapRowMap());
00306   }
00307 
00308   ptr = Lict->ptr;
00309   ind = Lict->col;
00310   val = Lict->val;
00311     
00312   for (i=0; i< m; i++) {
00313     int NumEntries = ptr[i+1]-ptr[i];
00314     int * Indices = ind+ptr[i];
00315     double * Values = val+ptr[i];
00316     U_->InsertMyValues(i, NumEntries, Values, Indices);
00317   }
00318 
00319   U_->FillComplete(A_.OperatorDomainMap(), A_.OperatorRangeMap());
00320   
00321   D_->Reciprocal(*D_); // Put reciprocal of diagonal in this vector
00322   // Add up flops
00323  
00324   double current_flops = 2 * nz; // Just an estimate
00325   double total_flops = 0;
00326     
00327   A_.Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs
00328 
00329   // Now count the rest
00330   total_flops += (double) U_->NumGlobalNonzeros(); // Accounts for multiplier above
00331   total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
00332 
00333   UpdateFlops(total_flops); // Update flop count
00334 
00335   SetFactored(true);
00336 
00337   return(0);
00338 
00339 }
00340 
00341 //=============================================================================
00342 int Ifpack_CrsIct::Solve(bool Trans, const Epetra_MultiVector& X, 
00343                 Epetra_MultiVector& Y) const {
00344 //
00345 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00346 //
00347 
00348   if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
00349 
00350   bool Upper = true;
00351   bool UnitDiagonal = true;
00352 
00353   Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
00354   Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
00355 
00356 
00357   U_->Solve(Upper, true, UnitDiagonal, *X1, *Y1);
00358   Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00359   U_->Solve(Upper, false, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
00360   return(0);
00361 }
00362 //=============================================================================
00363 int Ifpack_CrsIct::Multiply(bool Trans, const Epetra_MultiVector& X, 
00364                 Epetra_MultiVector& Y) const {
00365 //
00366 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00367 //
00368 
00369   if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
00370 
00371   //bool Upper = true;
00372   //bool Lower = false;
00373   //bool UnitDiagonal = true;
00374 
00375   Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
00376   Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
00377 
00378   U_->Multiply(false, *X1, *Y1);
00379   Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00380   Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00381   Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00382   U_->Multiply(true, Y1temp, *Y1);
00383   Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00384   return(0);
00385 }
00386 //=============================================================================
00387 int Ifpack_CrsIct::Condest(bool Trans, double & ConditionNumberEstimate) const {
00388 
00389   if (Condest_>=0.0) {
00390     ConditionNumberEstimate = Condest_;
00391     return(0);
00392   }
00393   // Create a vector with all values equal to one
00394   Epetra_Vector Ones(A_.RowMap());
00395   Epetra_Vector OnesResult(Ones);
00396   Ones.PutScalar(1.0);
00397 
00398   EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
00399   EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
00400   EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
00401   Condest_ = ConditionNumberEstimate; // Save value for possible later calls
00402   return(0);
00403 }
00404 //=============================================================================
00405 // Non-member functions
00406 
00407 ostream& operator << (ostream& os, const Ifpack_CrsIct& A)
00408 {
00409   // int LevelOverlap = A.LevelOverlap();
00410   Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
00411   Epetra_Vector & D = (Epetra_Vector &) A.D();
00412 
00413   //os.width(14);
00414   //os <<  "     Level of Overlap = "; os << LevelOverlap;
00415   //os << endl;
00416 
00417   os.width(14);
00418   os <<  "     Inverse of Diagonal = ";
00419   os << endl;
00420   os << D; // Let Epetra_Vector handle the rest.
00421   os << endl;
00422 
00423   os.width(14);
00424   os <<  "     Upper Triangle = ";
00425   os << endl;
00426   os << U; // Let Epetra_CrsMatrix handle the rest.
00427   os << endl;
00428  
00429   // Reset os flags
00430 
00431   return os;
00432 }

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