test_single_amesos_thyra_solver.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) 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 #include "test_single_amesos_thyra_solver.hpp"
00030 
00031 #ifndef SUN_CXX
00032 
00033 #include "Thyra_AmesosLinearOpWithSolveFactory.hpp"
00034 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00035 #include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
00036 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00037 #include "Thyra_DefaultInverseLinearOp.hpp"
00038 #include "Thyra_EpetraLinearOp.hpp"
00039 #include "Thyra_LinearOpTester.hpp"
00040 #include "Thyra_LinearOpWithSolveTester.hpp"
00041 #include "EpetraExt_readEpetraLinearSystem.h"
00042 #include "Epetra_SerialComm.h"
00043 
00044 #endif // SUN_CXX
00045 
00046 bool Thyra::test_single_amesos_thyra_solver(
00047   const std::string                       matrixFile
00048   ,Teuchos::ParameterList                 *amesosLOWSFPL
00049   ,const bool                             testTranspose_in
00050   ,const int                              numRandomVectors
00051   ,const double                           maxFwdError
00052   ,const double                           maxError
00053   ,const double                           maxResid
00054   ,const bool                             showAllTests
00055   ,const bool                             dumpAll
00056   ,Teuchos::FancyOStream                  *out_arg
00057   )
00058 {
00059 
00060   using Teuchos::RCP;
00061   using Teuchos::OSTab;
00062 
00063   bool result, success = true;
00064 
00065   RCP<Teuchos::FancyOStream>
00066     out = Teuchos::rcp(out_arg,false);
00067 
00068 #ifndef SUN_CXX
00069 
00070   if(out.get()) {
00071     *out
00072       << "\n***"
00073       << "\n*** Testing Thyra::AmesosLinearOpWithSolveFactory (and Thyra::AmesosLinearOpWithSolve)"
00074       << "\n***\n"
00075       << "\nEchoing input options:"
00076       << "\n  matrixFile             = " << matrixFile
00077       << "\n  testTranspose          = " << testTranspose_in
00078       << "\n  numRandomVectors       = " << numRandomVectors
00079       << "\n  maxFwdError            = " << maxFwdError
00080       << "\n  maxError               = " << maxError
00081       << "\n  maxResid               = " << maxResid
00082       << "\n  showAllTests           = " << showAllTests
00083       << "\n  dumpAll                = " << dumpAll
00084       << std::endl;
00085     if(amesosLOWSFPL) {
00086       OSTab tab(out);
00087       *out
00088         << "amesosLOWSFPL:\n";
00089       amesosLOWSFPL->print(OSTab(out).o(),0,true);
00090     }
00091   }
00092   
00093   if(out.get()) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
00094   
00095   Epetra_SerialComm comm;
00096   RCP<Epetra_CrsMatrix> epetra_A;
00097   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00098 
00099   RCP<const LinearOpBase<double> >
00100     A = Thyra::epetraLinearOp(epetra_A);
00101 
00102   if(out.get() && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
00103 
00104   if(out.get()) *out << "\nB) Creating a AmesosLinearOpWithSolveFactory object lowsFactory ...\n";
00105 
00106   RCP<LinearOpWithSolveFactoryBase<double> >
00107     lowsFactory = Teuchos::rcp(new AmesosLinearOpWithSolveFactory());
00108 
00109   lowsFactory->setParameterList(Teuchos::rcp(amesosLOWSFPL,false));
00110 
00111   if(out.get()) {
00112     *out << "\nlowsFactory.description() = " << lowsFactory->description() << std::endl;
00113     *out << "\nlowsFactory.getValidParameters() =\n";
00114     lowsFactory->getValidParameters()->print(OSTab(out).o(),0,true,false);
00115   }
00116 
00117   if(out.get()) *out << "\nC) Creating a AmesosLinearOpWithSolve object nsA ...\n";
00118 
00119   RCP<LinearOpWithSolveBase<double> >
00120     nsA = lowsFactory->createOp();
00121 
00122   initializeOp<double>( *lowsFactory, A, &*nsA );
00123 
00124   bool testTranspose = testTranspose_in;
00125   if (testTranspose_in && !Thyra::solveSupports(*nsA, Thyra::CONJTRANS)) {
00126     if(out.get()) *out
00127       << "\nChanging testTranspose=false since "
00128       << nsA->description() << " does not support adjoint solves!\n";
00129     testTranspose = false;
00130   }
00131   
00132   if(out.get()) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
00133 
00134   Thyra::seed_randomize<double>(0);
00135 
00136   LinearOpTester<double> linearOpTester;
00137   linearOpTester.check_adjoint(testTranspose);
00138   linearOpTester.num_random_vectors(numRandomVectors);
00139   linearOpTester.set_all_error_tol(maxFwdError);
00140   linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
00141   linearOpTester.show_all_tests(showAllTests);
00142   linearOpTester.dump_all(dumpAll);
00143   result = linearOpTester.check(*nsA,out.get());
00144   if(!result) success = false;
00145 
00146   if(out.get()) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00147     
00148   LinearOpWithSolveTester<double> linearOpWithSolveTester;
00149   linearOpWithSolveTester.turn_off_all_tests();
00150   linearOpWithSolveTester.check_forward_default(true);
00151   linearOpWithSolveTester.forward_default_residual_error_tol(1.1*maxResid);
00152   linearOpWithSolveTester.forward_default_residual_warning_tol(2.0*maxResid);
00153   linearOpWithSolveTester.check_forward_residual(true);
00154   linearOpWithSolveTester.forward_residual_solve_tol(maxResid);
00155   linearOpWithSolveTester.forward_residual_slack_error_tol(1e-1*maxResid);
00156   linearOpWithSolveTester.forward_residual_slack_warning_tol(maxResid);
00157   linearOpWithSolveTester.check_adjoint_default(testTranspose);
00158   linearOpWithSolveTester.adjoint_default_residual_error_tol(1.1*maxResid);
00159   linearOpWithSolveTester.adjoint_default_residual_warning_tol(2.0*maxResid);
00160   linearOpWithSolveTester.check_adjoint_residual(testTranspose);
00161   linearOpWithSolveTester.adjoint_residual_solve_tol(maxResid);
00162   linearOpWithSolveTester.adjoint_residual_slack_error_tol(1e-1*maxResid);
00163   linearOpWithSolveTester.adjoint_residual_slack_warning_tol(maxResid);
00164   linearOpWithSolveTester.num_random_vectors(numRandomVectors);
00165   linearOpWithSolveTester.show_all_tests(showAllTests);
00166   linearOpWithSolveTester.dump_all(dumpAll);
00167 
00168   LinearOpWithSolveTester<double> adjLinearOpWithSolveTester(linearOpWithSolveTester);
00169   adjLinearOpWithSolveTester.check_forward_default(testTranspose);
00170   adjLinearOpWithSolveTester.check_forward_residual(testTranspose);
00171 
00172   result = linearOpWithSolveTester.check(*nsA,out.get());
00173   if(!result) success = false;
00174 
00175   if(out.get()) *out << "\nF) Uninitialize the matrix object nsA, scale the epetra_A object by 2.5, and then refactor nsA with epetra_A ...\n";
00176 
00177   uninitializeOp<double>(*lowsFactory,&*nsA); // Optional call but a good idea if changing the operator
00178   epetra_A->Scale(2.5);
00179   initializeOp<double>(*lowsFactory,A,&*nsA);
00180   
00181   if(out.get()) *out << "\nG) Testing the LinearOpBase interface of nsA ...\n";
00182 
00183   Thyra::seed_randomize<double>(0);
00184 
00185   result = linearOpTester.check(*nsA,out.get());
00186   if(!result) success = false;
00187 
00188   if(out.get()) *out << "\nH) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00189     
00190   result = linearOpWithSolveTester.check(*nsA,out.get());
00191   if(!result) success = false;
00192 
00193   if(out.get()) *out << "\nI) Uninitialize the matrix object nsA, create a scaled (by 2.5) copy  epetra_A2 of epetra_A, and then refactor nsA with epetra_A2 ...\n";
00194 
00195   RCP<Epetra_CrsMatrix>
00196     epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
00197   epetra_A2->Scale(2.5);
00198   RCP<const LinearOpBase<double> >
00199     A2 = Thyra::epetraLinearOp(epetra_A2);
00200   initializeOp<double>(*lowsFactory,A2,&*nsA);
00201   
00202   if(out.get()) *out << "\nJ) Testing the LinearOpBase interface of nsA ...\n";
00203 
00204   Thyra::seed_randomize<double>(0);
00205 
00206   result = linearOpTester.check(*nsA,out.get());
00207   if(!result) success = false;
00208 
00209   if(out.get()) *out << "\nK) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00210     
00211   result = linearOpWithSolveTester.check(*nsA,out.get());
00212   if(!result) success = false;
00213 
00214   if(out.get()) *out << "\nL) Create an implicitly scaled (by 2.5) and transposed matrix A3 = scale(2.5,transpose(A)) and initialize nsA2 ...\n";
00215 
00216   RCP<const LinearOpBase<double> >
00217     A3 = scale<double>(2.5,Thyra::transpose<double>(A));
00218   RCP<LinearOpWithSolveBase<double> >
00219     nsA2 = linearOpWithSolve(*lowsFactory,A3);
00220   
00221   if(out.get()) *out << "\nM) Testing the LinearOpBase interface of nsA2 ...\n";
00222 
00223   Thyra::seed_randomize<double>(0);
00224 
00225   result = linearOpTester.check(*nsA2,out.get());
00226   if(!result) success = false;
00227 
00228   if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
00229     
00230   result = adjLinearOpWithSolveTester.check(*nsA2,out.get());
00231   if(!result) success = false;
00232   
00233   if(out.get()) *out << "\nO) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
00234 
00235   result = linearOpTester.compare(
00236     *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
00237     ,out.get()
00238     );
00239   if(!result) success = false;
00240 
00241   if(out.get()) *out << "\nP) Running example use cases ...\n";
00242 
00243   nonExternallyPreconditionedLinearSolveUseCases(
00244     *A,*lowsFactory,true,*out
00245     );
00246 
00247   if(out.get()) *out << "\nQ) Creating a DefaultInverseLinearOp object from nsA and testing the LinearOpBase interface ...\n";
00248 
00249   RCP<const LinearOpBase<double> >
00250     invA = inverse<double>(nsA.getConst());
00251 
00252   result = linearOpTester.check(*invA,out.get());
00253   if(!result) success = false;
00254 
00255 #else // SUN_CXX
00256   
00257   if(out.get()) *out << "\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
00258   success = false;
00259 
00260 #endif // SUN_CXX
00261 
00262   return success;
00263 
00264 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:19:51 2011 for Stratimikos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3