Testing Program for the Amesos/Thyra Linear Solver Adapters.
[Assorted Amesos to Thyra Linear Solver Adapter Example Code]

Here we show a testing program that takes the Amesos/%Thyra adatpers through the paces. More...
Below is a complete listing of the testing function with line numbers:

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
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
00045 
00046 bool Thyra::test_single_amesos_thyra_solver(
00047   const std::string                       matrixFile
00048   ,Teuchos::ParameterList                 *amesosLOWSFPL
00049   ,const bool                             testTranspose
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   using Teuchos::OSTab;
00060 
00061   bool result, success = true;
00062 
00063   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00064     out = Teuchos::rcp(out_arg,false);
00065 
00066 #ifndef __sun
00067 
00068   if(out.get()) {
00069     *out
00070       << "\n***"
00071       << "\n*** Testing Thyra::AmesosLinearOpWithSolveFactory (and Thyra::AmesosLinearOpWithSolve)"
00072       << "\n***\n"
00073       << "\nEchoing input options:"
00074       << "\n  matrixFile             = " << matrixFile
00075       << "\n  testTranspose          = " << testTranspose
00076       << "\n  numRandomVectors       = " << numRandomVectors
00077       << "\n  maxFwdError            = " << maxFwdError
00078       << "\n  maxError               = " << maxError
00079       << "\n  maxResid               = " << maxResid
00080       << "\n  showAllTests           = " << showAllTests
00081       << "\n  dumpAll                = " << dumpAll
00082       << std::endl;
00083     if(amesosLOWSFPL) {
00084       OSTab tab(out);
00085       *out
00086         << "amesosLOWSFPL:\n";
00087       amesosLOWSFPL->print(*OSTab(out).getOStream(),0,true);
00088     }
00089   }
00090   
00091   if(out.get()) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
00092   
00093   Epetra_SerialComm comm;
00094   Teuchos::RefCountPtr<Epetra_CrsMatrix> epetra_A;
00095   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00096 
00097   Teuchos::RefCountPtr<LinearOpBase<double> >
00098     A = Teuchos::rcp(new EpetraLinearOp(epetra_A));
00099 
00100   if(out.get() && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
00101 
00102   if(out.get()) *out << "\nB) Creating a AmesosLinearOpWithSolveFactory object lowsFactory ...\n";
00103 
00104   Teuchos::RefCountPtr<LinearOpWithSolveFactoryBase<double> >
00105     lowsFactory = Teuchos::rcp(new AmesosLinearOpWithSolveFactory());
00106 
00107   lowsFactory->setParameterList(Teuchos::rcp(amesosLOWSFPL,false));
00108 
00109   if(out.get()) {
00110     *out << "\nlowsFactory.description() = " << lowsFactory->description() << std::endl;
00111     *out << "\nlowsFactory.getValidParameters() =\n";
00112     lowsFactory->getValidParameters()->print(*OSTab(out).getOStream(),0,true,false);
00113   }
00114 
00115   if(out.get()) *out << "\nC) Creating a AmesosLinearOpWithSolve object nsA ...\n";
00116 
00117   Teuchos::RefCountPtr<LinearOpWithSolveBase<double> >
00118     nsA = lowsFactory->createOp();
00119 
00120   initializeOp<double>( *lowsFactory, A, &*nsA );
00121 
00122   if(out.get()) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
00123 
00124   Thyra::seed_randomize<double>(0);
00125 
00126   LinearOpTester<double> linearOpTester;
00127   linearOpTester.check_adjoint(testTranspose);
00128   linearOpTester.num_random_vectors(numRandomVectors);
00129   linearOpTester.set_all_error_tol(maxFwdError);
00130   linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
00131   linearOpTester.show_all_tests(showAllTests);
00132   linearOpTester.dump_all(dumpAll);
00133   result = linearOpTester.check(*nsA,out.get());
00134   if(!result) success = false;
00135 
00136   if(out.get()) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00137     
00138   LinearOpWithSolveTester<double> linearOpWithSolveTester;
00139   linearOpWithSolveTester.turn_off_all_tests();
00140   linearOpWithSolveTester.check_forward_default(true);
00141   linearOpWithSolveTester.forward_default_residual_error_tol(1.1*maxResid);
00142   linearOpWithSolveTester.forward_default_residual_warning_tol(2.0*maxResid);
00143   linearOpWithSolveTester.check_forward_residual(true);
00144   linearOpWithSolveTester.forward_residual_solve_tol(maxResid);
00145   linearOpWithSolveTester.forward_residual_slack_error_tol(1e-1*maxResid);
00146   linearOpWithSolveTester.forward_residual_slack_warning_tol(maxResid);
00147   linearOpWithSolveTester.check_forward_solution_error(true);
00148   linearOpWithSolveTester.forward_solution_error_solve_tol(maxError);
00149   linearOpWithSolveTester.forward_solution_error_slack_error_tol(1e-1*maxError);
00150   linearOpWithSolveTester.forward_solution_error_slack_warning_tol(maxError);
00151   if(testTranspose) {
00152     linearOpWithSolveTester.check_adjoint_default(true);
00153     linearOpWithSolveTester.adjoint_default_residual_error_tol(1.1*maxResid);
00154     linearOpWithSolveTester.adjoint_default_residual_warning_tol(2.0*maxResid);
00155     linearOpWithSolveTester.check_adjoint_residual(true);
00156     linearOpWithSolveTester.adjoint_residual_solve_tol(maxResid);
00157     linearOpWithSolveTester.adjoint_residual_slack_error_tol(1e-1*maxResid);
00158     linearOpWithSolveTester.adjoint_residual_slack_warning_tol(maxResid);
00159     linearOpWithSolveTester.check_adjoint_solution_error(true);
00160     linearOpWithSolveTester.adjoint_solution_error_solve_tol(maxError);
00161     linearOpWithSolveTester.adjoint_solution_error_slack_error_tol(1e-1*maxError);
00162     linearOpWithSolveTester.adjoint_solution_error_slack_warning_tol(maxError);
00163   }
00164   linearOpWithSolveTester.num_random_vectors(numRandomVectors);
00165   linearOpWithSolveTester.show_all_tests(showAllTests);
00166   linearOpWithSolveTester.dump_all(dumpAll);
00167   result = linearOpWithSolveTester.check(*nsA,out.get());
00168   if(!result) success = false;
00169 
00170   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";
00171 
00172   uninitializeOp<double>(*lowsFactory,&*nsA); // Optional call but a good idea if changing the operator
00173   epetra_A->Scale(2.5);
00174   initializeOp<double>(*lowsFactory,A,&*nsA);
00175   
00176   if(out.get()) *out << "\nG) Testing the LinearOpBase interface of nsA ...\n";
00177 
00178   Thyra::seed_randomize<double>(0);
00179 
00180   result = linearOpTester.check(*nsA,out.get());
00181   if(!result) success = false;
00182 
00183   if(out.get()) *out << "\nH) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00184     
00185   result = linearOpWithSolveTester.check(*nsA,out.get());
00186   if(!result) success = false;
00187 
00188   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";
00189 
00190   Teuchos::RefCountPtr<Epetra_CrsMatrix>
00191     epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
00192   epetra_A2->Scale(2.5);
00193   Teuchos::RefCountPtr<LinearOpBase<double> >
00194     A2 = Teuchos::rcp(new EpetraLinearOp(epetra_A2));
00195   initializeOp<double>(*lowsFactory,A2,&*nsA);
00196   
00197   if(out.get()) *out << "\nJ) Testing the LinearOpBase interface of nsA ...\n";
00198 
00199   Thyra::seed_randomize<double>(0);
00200 
00201   result = linearOpTester.check(*nsA,out.get());
00202   if(!result) success = false;
00203 
00204   if(out.get()) *out << "\nK) Testing the LinearOpWithSolveBase interface of nsA ...\n";
00205     
00206   result = linearOpWithSolveTester.check(*nsA,out.get());
00207   if(!result) success = false;
00208 
00209   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";
00210 
00211   Teuchos::RefCountPtr<const LinearOpBase<double> >
00212     A3 = scale<double>(2.5,Thyra::transpose<double>(A));
00213   Teuchos::RefCountPtr<LinearOpWithSolveBase<double> >
00214     nsA2 = linearOpWithSolve(*lowsFactory,A3);
00215   
00216   if(out.get()) *out << "\nM) Testing the LinearOpBase interface of nsA2 ...\n";
00217 
00218   Thyra::seed_randomize<double>(0);
00219 
00220   result = linearOpTester.check(*nsA2,out.get());
00221   if(!result) success = false;
00222 
00223   if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
00224     
00225   result = linearOpWithSolveTester.check(*nsA2,out.get());
00226   if(!result) success = false;
00227   
00228   if(out.get()) *out << "\nO) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
00229 
00230   result = linearOpTester.compare(
00231     *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
00232     ,out.get()
00233     );
00234   if(!result) success = false;
00235 
00236   if(out.get()) *out << "\nP) Running example use cases ...\n";
00237 
00238   nonExternallyPreconditionedLinearSolveUseCases(
00239     *A,*lowsFactory,true,*out
00240     );
00241 
00242   if(out.get()) *out << "\nQ) Creating a DefaultInverseLinearOp object from nsA and testing the LinearOpBase interface ...\n";
00243 
00244   Teuchos::RefCountPtr<const LinearOpBase<double> >
00245     invA = inverse<double>(nsA);
00246 
00247   result = linearOpTester.check(*invA,out.get());
00248   if(!result) success = false;
00249 
00250 #else // __sun
00251   
00252   if(out.get()) *out << "\nTest failed since is was not even compiled since __sun was defined!\n";
00253   success = false;
00254 
00255 #endif // __sun
00256 
00257   return success;
00258 
00259 }

Generated on Thu Sep 18 12:30:59 2008 for Amesos/Thyra Linear Solver Adapter Software by doxygen 1.3.9.1