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

Generated on Tue Jul 13 09:27:13 2010 for IFPACK by  doxygen 1.4.7