test/IC/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_CrsIct.h"
00050 #include "Ifpack.h"
00051 
00052 // function for fancy output
00053 
00054 string toString(const int& x) {
00055   char s[100];
00056   sprintf(s, "%d", x);
00057   return string(s);
00058 }
00059 
00060 string toString(const double& x) {
00061   char s[100];
00062   sprintf(s, "%g", x);
00063   return string(s);
00064 }
00065 
00066 // main driver
00067 
00068 int main(int argc, char *argv[]) {
00069 
00070 #ifdef HAVE_MPI
00071   MPI_Init(&argc,&argv);
00072   Epetra_MpiComm Comm (MPI_COMM_WORLD);
00073 #else
00074   Epetra_SerialComm Comm;
00075 #endif
00076 
00077   int MyPID = Comm.MyPID();
00078   bool verbose = false; 
00079   if (MyPID==0) verbose = true;
00080 
00081   // The problem is defined on a 2D grid, global size is nx * nx.
00082   int nx = 30;
00083   Teuchos::ParameterList GaleriList;
00084   GaleriList.set("nx", nx);
00085   GaleriList.set("ny", nx * Comm.NumProc());
00086   GaleriList.set("mx", 1);
00087   GaleriList.set("my", Comm.NumProc());
00088   Epetra_Map* Map = Galeri::CreateMap("Cartesian2D", Comm, GaleriList);
00089   Epetra_CrsMatrix* A = Galeri::CreateCrsMatrix("Laplace2D", Map, GaleriList);
00090   Epetra_MultiVector* LHS = new Epetra_MultiVector(*Map, 1);
00091   Epetra_MultiVector* RHS = new Epetra_MultiVector(*Map, 1);
00092   LHS->PutScalar(0.0); RHS->Random();
00093 
00094   // ============================ //
00095   // Construct ILU preconditioner //
00096   // ---------------------------- //
00097 
00098   // I wanna test funky values to be sure that they have the same
00099   // influence on the algorithms, both old and new
00100   int    LevelFill = 2;
00101   double DropTol = 0.3333;
00102   double Condest;
00103   
00104   Ifpack_CrsIct * ICT = NULL;
00105   ICT = new Ifpack_CrsIct(*A,DropTol,LevelFill);
00106   ICT->SetAbsoluteThreshold(0.00123);
00107   ICT->SetRelativeThreshold(0.9876);
00108   // Init values from A
00109   ICT->InitValues(*A);
00110   // compute the factors
00111   ICT->Factor();
00112   // and now estimate the condition number
00113   ICT->Condest(false,Condest);
00114   
00115   if( Comm.MyPID() == 0 ) {
00116     cout << "Condition number estimate (level-of-fill = "
00117    << LevelFill <<  ") = " << Condest << endl;
00118   }
00119 
00120   // Define label for printing out during the solve phase
00121   string label = "Ifpack_CrsIct Preconditioner: LevelFill = " + toString(LevelFill) + 
00122                                                  " Overlap = 0"; 
00123   ICT->SetLabel(label.c_str());
00124   
00125   // Here we create an AztecOO object
00126   LHS->PutScalar(0.0);
00127 
00128   int Niters = 1200;
00129 
00130   AztecOO solver;
00131   solver.SetUserMatrix(A);
00132   solver.SetLHS(LHS);
00133   solver.SetRHS(RHS);
00134   solver.SetAztecOption(AZ_solver,AZ_cg);
00135   solver.SetPrecOperator(ICT);
00136   solver.SetAztecOption(AZ_output, 16); 
00137   solver.Iterate(Niters, 5.0e-5);
00138 
00139   int OldIters = solver.NumIters();
00140 
00141   if (ICT!=0) delete ICT;
00142                
00143   // now rebuild the same preconditioner using ICT, we expect the same
00144   // number of iterations
00145 
00146   Ifpack Factory;
00147   Ifpack_Preconditioner* Prec = Factory.Create("IC", A);
00148 
00149   Teuchos::ParameterList List;
00150   List.get("fact: level-of-fill", 2);
00151   List.get("fact: drop tolerance", 0.3333);
00152   List.get("fact: absolute threshold", 0.00123);
00153   List.get("fact: relative threshold", 0.9876);
00154   List.get("fact: relaxation value", 0.0);
00155 
00156   IFPACK_CHK_ERR(Prec->SetParameters(List));
00157   IFPACK_CHK_ERR(Prec->Compute());
00158 
00159   // Here we create an AztecOO object
00160   LHS->PutScalar(0.0);
00161 
00162   solver.SetUserMatrix(A);
00163   solver.SetLHS(LHS);
00164   solver.SetRHS(RHS);
00165   solver.SetAztecOption(AZ_solver,AZ_cg);
00166   solver.SetPrecOperator(Prec);
00167   solver.SetAztecOption(AZ_output, 16); 
00168   solver.Iterate(Niters, 5.0e-5);
00169 
00170   int NewIters = solver.NumIters();
00171 
00172   if (OldIters != NewIters)
00173     IFPACK_CHK_ERR(-1);
00174 
00175   delete Prec;
00176   delete A;
00177   delete LHS;
00178   delete RHS;
00179   delete Map;
00180 
00181 #ifdef HAVE_MPI
00182   MPI_Finalize() ;
00183 #endif
00184 
00185   return(EXIT_SUCCESS);
00186 }

Generated on Thu Sep 18 12:37:21 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1