aztecoo_solve.cpp

Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) 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 "Teuchos_GlobalMPISession.hpp"
00031 #include "Teuchos_CommandLineProcessor.hpp"
00032 #include "Teuchos_ParameterList.hpp"
00033 #include "Teuchos_ParameterXMLFileReader.hpp"
00034 #include "Teuchos_RefCountPtr.hpp"
00035 #include "Teuchos_Time.hpp"
00036 #include "Teuchos_Comm.hpp"
00037 
00038 #ifdef HAVE_MPI
00039 #include "Epetra_MpiComm.h"
00040 #else
00041 #include "Epetra_SerialComm.h"
00042 #endif
00043 
00044 #include "Epetra_LinearProblem.h"
00045 
00046 #include "AztecOO.h"
00047 
00048 #include "ParameterHelper.hpp"
00049 
00050 #include "build_problem.hpp"
00051 #include "build_solver.hpp"
00052 
00053 void process_command_line(int argc, char*argv[], std::string& xml_file);
00054 
00055 int main(int argc, char*argv[])
00056 {
00057   Teuchos::GlobalMPISession mpisess(&argc,&argv,&std::cout);
00058 
00059   Teuchos::Time timer("total");
00060   timer.start();
00061 
00062 #ifdef HAVE_MPI
00063   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00064 #else
00065   Epetra_SerialComm Comm;
00066 #endif
00067 
00068   //Just get one parameter from the command-line: the name of an xml file
00069   //to get parameters from.
00070 
00071   std::string xml_file("xml_params.xml");
00072   process_command_line(argc, argv, xml_file);
00073 
00074   //Read the contents of the xml file into a ParameterList. That parameter list
00075   //should specify a matrix-file and optionally which Belos solver to use, and
00076   //which Tifpack preconditioner to use, etc. If there are sublists of parameters
00077   //for Belos and Tifpack, those will be passed to the respective destinations
00078   //from within the build_problem and build_solver functions.
00079 
00080   std::cout << "Every proc reading parameters from xml_file: "
00081             << xml_file << std::endl;
00082   Teuchos::ParameterList test_params =
00083       Teuchos::ParameterXMLFileReader(xml_file).getParameters();
00084 
00085   //The build_problem function is declared in build_problem.hpp.
00086   //Note that build_problem calls build_precond and sets a preconditioner on the
00087   //linear-problem, if a preconditioner is specified.
00088 
00089   Teuchos::RCP<Epetra_LinearProblem> problem = build_problem(test_params, Comm);
00090 
00091   //The build_solver function is declared in build_solver.hpp:
00092 
00093   Teuchos::RCP<AztecOO> solver = build_solver(test_params, problem);
00094 
00095   int max_iters = solver->GetAllAztecOptions()[AZ_max_iter];
00096   double tol    = solver->GetAllAztecParams()[AZ_tol];
00097 
00098   Teuchos::Time prec_time("precond");
00099   prec_time.start();
00100   double cond_est = 0;
00101   int err_code = -1;
00102   if (solver->GetAllAztecOptions()[AZ_precond] != AZ_user_precond) {
00103     err_code = solver->ConstructPreconditioner(cond_est);
00104   }
00105   prec_time.stop();
00106   if (Comm.MyPID() == 0 && err_code == 0) {
00107     std::cout << "Time to compute preconditioner: " << prec_time.totalElapsedTime() << "s"<<std::endl;
00108   }
00109 
00110   int ret = solver->Iterate(max_iters, tol);
00111 
00112   if (err_code == 0) {
00113     solver->DestroyPreconditioner();
00114   }
00115 
00116   int actual_iters = (int)solver->GetAztecStatus()[AZ_its];
00117 
00118   if (Comm.MyPID() == 0) {
00119     std::cout << "Converged in " << actual_iters << " iterations." << std::endl;
00120   }
00121 
00122   if (problem->GetLHS()->NumVectors() > 1) {
00123     throw std::runtime_error("ERROR: MultiVector->NumVectors()>1.");
00124   }
00125 
00126   Epetra_MultiVector* x = problem->GetLHS();
00127   Epetra_MultiVector* b = problem->GetRHS();
00128   Epetra_Operator* A = problem->GetOperator();
00129   Teuchos::RCP<Epetra_MultiVector> r = Teuchos::rcp(new Epetra_MultiVector(*b));
00130   A->Apply(*x, *r);
00131   r->Update(1.0, *b, -1.0);
00132   double norm = 0;
00133   r->Norm2(&norm);
00134 
00135   if (Comm.MyPID() == 0) {
00136     std::cout << "2-Norm of residual vec: " << norm << std::endl;
00137   }
00138 
00139   //Clean up by destroying heap-allocated objects that are not
00140   //held by Teuchos::RCP:
00141   delete A;
00142   delete x;
00143   delete b;
00144   Epetra_Operator* prec = solver->GetPrecOperator();
00145   delete prec;
00146 
00147   //If the xml file specified a number of iterations to expect, then we will
00148   //use that as a test pass/fail criteria.
00149 
00150   if (test_params.isParameter("expectNumIters")) {
00151     int expected_iters = 0;
00152     helper::GetParameter(test_params, "expectNumIters", expected_iters);
00153     double eps = 2.e-6;
00154     if (ret == 0 && actual_iters <= expected_iters && norm < eps) {
00155       if (Comm.MyPID() == 0) {
00156         std::cout << "End Result: TEST PASSED" << std::endl;
00157       }
00158     }
00159     else {
00160       if (Comm.MyPID() == 0) {
00161         std::cout << "Actual iters("<<actual_iters
00162            <<") > expected number of iterations ("
00163               <<expected_iters<<"), or resid-norm(" << norm << ") >= "<<eps <<std::endl;
00164       }
00165     }
00166   }
00167 
00168   timer.stop();
00169   if (Comm.MyPID() == 0) {
00170     std::cout << "proc 0 total program time: " << timer.totalElapsedTime() << std::endl;
00171   }
00172 
00173   return 0;
00174 }
00175 
00176 void process_command_line(int argc, char*argv[], std::string& xml_file)
00177 {
00178   Teuchos::CommandLineProcessor cmdp(false,true);
00179   cmdp.setOption("xml_file", &xml_file, "XML Parameters file");
00180   if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
00181     throw std::runtime_error("Error parsing command-line.");
00182   }
00183 }
00184 
 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