test/ILU/cxx_main.cpp

Go to the documentation of this file.
00001 
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //            Trilinos: An Object-Oriented Solver Framework
00006 //                 Copyright (2001) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #include "Ifpack_ConfigDefs.h"
00031 
00032 #include "Epetra_ConfigDefs.h"
00033 #ifdef HAVE_MPI
00034 #include "mpi.h"
00035 #include "Epetra_MpiComm.h"
00036 #else
00037 #include "Epetra_SerialComm.h"
00038 #endif
00039 #include "Epetra_Comm.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_Time.h"
00042 #include "Epetra_BlockMap.h"
00043 #include "Epetra_MultiVector.h"
00044 #include "Epetra_Vector.h"
00045 #include "Epetra_Export.h"
00046 #include "AztecOO.h"
00047 #include "Galeri_Maps.h"
00048 #include "Galeri_CrsMatrices.h"
00049 #include "Ifpack_CrsRiluk.h"
00050 #include "Ifpack.h"
00051 #include "Teuchos_RefCountPtr.hpp"
00052 
00053 // function for fancy output
00054 
00055 string toString(const int& x) {
00056   char s[100];
00057   sprintf(s, "%d", x);
00058   return string(s);
00059 }
00060 
00061 string toString(const double& x) {
00062   char s[100];
00063   sprintf(s, "%g", x);
00064   return string(s);
00065 }
00066 
00067 // main driver
00068 
00069 int main(int argc, char *argv[]) {
00070 
00071 #ifdef HAVE_MPI
00072   MPI_Init(&argc,&argv);
00073   Epetra_MpiComm Comm (MPI_COMM_WORLD);
00074 #else
00075   Epetra_SerialComm Comm;
00076 #endif
00077 
00078   int MyPID = Comm.MyPID();
00079   bool verbose = false; 
00080   if (MyPID==0) verbose = true;
00081 
00082   Teuchos::ParameterList GaleriList;
00083   int nx = 30; 
00084 
00085   GaleriList.set("nx", nx);
00086   GaleriList.set("ny", nx * Comm.NumProc());
00087   GaleriList.set("mx", 1);
00088   GaleriList.set("my", Comm.NumProc());
00089   Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
00090   Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
00091   Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
00092   Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
00093   LHS->PutScalar(0.0); RHS->Random();
00094 
00095   // ============================ //
00096   // Construct ILU preconditioner //
00097   // ---------------------------- //
00098 
00099   // I wanna test funky values to be sure that they have the same
00100   // influence on the algorithms, both old and new
00101   int    LevelFill = 2;
00102   double DropTol = 0.3333;
00103   double Athresh = 0.0123;
00104   double Rthresh = 0.9876;
00105   double Relax   = 0.1;
00106   int    Overlap = 2;
00107   
00108   Teuchos::RefCountPtr<Ifpack_IlukGraph> Graph;
00109   Teuchos::RefCountPtr<Ifpack_CrsRiluk> RILU;
00110 
00111   Graph = Teuchos::rcp( new Ifpack_IlukGraph(A->Graph(), LevelFill, Overlap) );
00112   int ierr;
00113   ierr = Graph->ConstructFilledGraph();
00114   IFPACK_CHK_ERR(ierr);
00115 
00116   RILU = Teuchos::rcp( new Ifpack_CrsRiluk(*Graph) );
00117   RILU->SetAbsoluteThreshold(Athresh);
00118   RILU->SetRelativeThreshold(Rthresh);
00119   RILU->SetRelaxValue(Relax);
00120   int initerr = RILU->InitValues(*A);
00121   if (initerr!=0) cout << Comm << "*ERR* InitValues = " << initerr;
00122 
00123   RILU->Factor();
00124 
00125   // Define label for printing out during the solve phase
00126   string label = "Ifpack_CrsRiluk Preconditioner: LevelFill = " + toString(LevelFill) +
00127                                                  " Overlap = " + toString(Overlap) +
00128                                                  " Athresh = " + toString(Athresh) +
00129                                                  " Rthresh = " + toString(Rthresh);
00130   // Here we create an AztecOO object
00131   LHS->PutScalar(0.0);
00132 
00133   int Niters = 1200;
00134 
00135   AztecOO solver;
00136   solver.SetUserMatrix(&*A);
00137   solver.SetLHS(&*LHS);
00138   solver.SetRHS(&*RHS);
00139   solver.SetAztecOption(AZ_solver,AZ_gmres);
00140   solver.SetPrecOperator(&*RILU);
00141   solver.SetAztecOption(AZ_output, 16); 
00142   solver.Iterate(Niters, 5.0e-5);
00143 
00144   int OldIters = solver.NumIters();
00145 
00146   // now rebuild the same preconditioner using RILU, we expect the same
00147   // number of iterations
00148 
00149   Ifpack Factory;
00150   Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create("ILU", &*A, Overlap) );
00151 
00152   Teuchos::ParameterList List;
00153   List.get("fact: level-of-fill", LevelFill);
00154   List.get("fact: drop tolerance", DropTol);
00155   List.get("fact: absolute threshold", Athresh);
00156   List.get("fact: relative threshold", Rthresh);
00157   List.get("fact: relax value", Relax);
00158 
00159   IFPACK_CHK_ERR(Prec->SetParameters(List));
00160   IFPACK_CHK_ERR(Prec->Compute());
00161 
00162   // Here we create an AztecOO object
00163   LHS->PutScalar(0.0);
00164 
00165   solver.SetUserMatrix(&*A);
00166   solver.SetLHS(&*LHS);
00167   solver.SetRHS(&*RHS);
00168   solver.SetAztecOption(AZ_solver,AZ_gmres);
00169   solver.SetPrecOperator(&*Prec);
00170   solver.SetAztecOption(AZ_output, 16); 
00171   solver.Iterate(Niters, 5.0e-5);
00172 
00173   int NewIters = solver.NumIters();
00174 
00175   if (OldIters != NewIters)
00176     IFPACK_CHK_ERR(-1);
00177 
00178 
00179 #ifdef HAVE_IFPACK_SUPERLU
00180   // Now test w/ SuperLU's ILU, if we've got it
00181   Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec2 = Teuchos::rcp( Factory.Create("SILU", &*A,0) ); 
00182   Teuchos::ParameterList SList;
00183   SList.set("fact: drop tolerance",1e-4);
00184   SList.set("fact: zero pivot threshold",.1);
00185   SList.set("fact: maximum fill factor",10.0);
00186   // NOTE: There is a bug in SuperLU 4.0 which will crash the code if the maximum fill factor is set too low.
00187   // This bug was reported to Sherry Li on 4/8/10.
00188   SList.set("fact: silu drop rule",9);
00189 
00190   IFPACK_CHK_ERR(Prec2->SetParameters(SList));
00191   IFPACK_CHK_ERR(Prec2->Compute());
00192 
00193   LHS->PutScalar(0.0);
00194   solver.SetPrecOperator(&*Prec2);
00195   solver.Iterate(Niters, 5.0e-5);
00196   Prec2->Print(cout);
00197 
00198 #endif
00199 
00200 
00201 #ifdef HAVE_MPI
00202   MPI_Finalize() ;
00203 #endif
00204 
00205   return(EXIT_SUCCESS);
00206 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:33 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3