test_belos_gcrodr.cpp

Go to the documentation of this file.
00001 #include <Teuchos_ConfigDefs.hpp>
00002 #include <Teuchos_UnitTestHarness.hpp>
00003 #include <Teuchos_RCP.hpp>
00004 #include <Teuchos_ParameterList.hpp>
00005 #include "Teuchos_XMLParameterListHelpers.hpp"
00006 
00007 #include <string>
00008 #include <iostream>
00009 
00010 #ifdef HAVE_MPI
00011    #include "Epetra_MpiComm.h"
00012 #else
00013    #include "Epetra_SerialComm.h"
00014 #endif
00015 
00016 #include "Epetra_Map.h"
00017 #include "Epetra_CrsMatrix.h"
00018 #include "Epetra_Vector.h"
00019 // #include "EpetraExt_RowMatrixOut.h"
00020 
00021 #include "Thyra_LinearOpBase.hpp"
00022 #include "Thyra_VectorBase.hpp"
00023 #include "Thyra_EpetraThyraWrappers.hpp"
00024 #include "Thyra_EpetraLinearOp.hpp"
00025 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00026 #include "Thyra_LinearOpWithSolveHelpers.hpp"
00027 #include "Thyra_DefaultZeroLinearOp.hpp"
00028 #include "Thyra_DefaultProductVector.hpp"
00029 #include "Thyra_DefaultProductVectorSpace.hpp"
00030 #include "Thyra_MultiVectorStdOps.hpp"
00031 #include "Thyra_VectorStdOps.hpp"
00032 #include "Thyra_DefaultBlockedLinearOp.hpp"
00033 
00034 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00035 
00036 using Teuchos::RCP;
00037 using Teuchos::rcp;
00038 
00039 Teuchos::RCP<Epetra_CrsMatrix> buildMatrix(int nx, Epetra_Comm & comm)
00040 {
00041    Epetra_Map map(nx*comm.NumProc(),0,comm);
00042    Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(new Epetra_CrsMatrix(Copy,map,3));
00043 
00044    int offsets[3] = {-1, 0, 1 };
00045    double values[3] = { -1, 2, -1};
00046    int maxGid = map.MaxAllGID();
00047    for(int lid=0;lid<nx;lid++) {
00048       int gid = mat->GRID(lid);
00049       int numEntries = 3, offset = 0;
00050       int indices[3] = { gid+offsets[0],
00051                          gid+offsets[1],
00052                          gid+offsets[2] };
00053       if(gid==0) { // left end point
00054          numEntries = 2;
00055          offset = 1;
00056       }            // right end point
00057       else if(gid==maxGid)
00058          numEntries = 2;
00059 
00060       // insert rows
00061       mat->InsertGlobalValues(gid,numEntries,values+offset,indices+offset);
00062    }
00063 
00064    mat->FillComplete();
00065    return mat;
00066 }
00067 
00068 TEUCHOS_UNIT_TEST(belos_gcrodr, multiple_solves)
00069 {
00070    // build global (or serial communicator)
00071    #ifdef HAVE_MPI
00072       Epetra_MpiComm Comm(MPI_COMM_WORLD);
00073    #else
00074       Epetra_SerialComm Comm;
00075    #endif
00076 
00077    // build and allocate linear system
00078    Teuchos::RCP<Epetra_CrsMatrix> mat = buildMatrix(100,Comm);
00079    Teuchos::RCP<Epetra_Vector> x0 = rcp(new Epetra_Vector(mat->OperatorDomainMap()));
00080    Teuchos::RCP<Epetra_Vector> x1 = rcp(new Epetra_Vector(mat->OperatorDomainMap()));
00081    Teuchos::RCP<Epetra_Vector> b = rcp(new Epetra_Vector(mat->OperatorRangeMap()));
00082 
00083    x0->Random();
00084    x1->Random();
00085    b->PutScalar(0.0);
00086 
00087    // sanity check
00088    // EpetraExt::RowMatrixToMatrixMarketFile("mat_output.mm",*mat);
00089 
00090    // build Thyra wrappers
00091    RCP<const Thyra::LinearOpBase<double> >
00092       tA = Thyra::epetraLinearOp( mat );
00093    RCP<Thyra::VectorBase<double> >
00094       tx0 = Thyra::create_Vector( x0, tA->domain() );
00095    RCP<Thyra::VectorBase<double> >
00096       tx1 = Thyra::create_Vector( x1, tA->domain() );
00097    RCP<const Thyra::VectorBase<double> >
00098       tb = Thyra::create_Vector( b, tA->range() );
00099 
00100    // now comes Stratimikos
00101    RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("BelosGCRODRTest.xml");
00102    Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00103    linearSolverBuilder.setParameterList(paramList);
00104  
00105    // Create a linear solver factory given information read from the
00106    // parameter list.
00107    RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
00108          linearSolverBuilder.createLinearSolveStrategy("");
00109 
00110    // Create a linear solver based on the forward operator A
00111    RCP<Thyra::LinearOpWithSolveBase<double> > lows =
00112          Thyra::linearOpWithSolve(*lowsFactory, tA);
00113 
00114    // Solve the linear system 
00115    Thyra::SolveStatus<double> status; 
00116    status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx0.ptr());
00117    status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx1.ptr());
00118 }
00119 
00120 TEUCHOS_UNIT_TEST(belos_gcrodr, 2x2_multiple_solves)
00121 {
00122    // build global (or serial communicator)
00123    #ifdef HAVE_MPI
00124       Epetra_MpiComm Comm(MPI_COMM_WORLD);
00125    #else
00126       Epetra_SerialComm Comm;
00127    #endif
00128 
00129    // build and allocate linear system
00130    Teuchos::RCP<Epetra_CrsMatrix> mat = buildMatrix(100,Comm);
00131    Teuchos::RCP<Epetra_Vector> b = rcp(new Epetra_Vector(mat->OperatorRangeMap()));
00132 
00133    b->PutScalar(0.0);
00134 
00135    // sanity check
00136    // EpetraExt::RowMatrixToMatrixMarketFile("mat_output.mm",*mat);
00137 
00138    // build Thyra wrappers
00139    RCP<const Thyra::LinearOpBase<double> > tA;
00140    RCP<const Thyra::VectorBase<double> > tb;
00141    {
00142       // build blocked linear Op
00143       RCP<const Thyra::LinearOpBase<double> > tA_sub 
00144             = Thyra::epetraLinearOp( mat );
00145       RCP<const Thyra::LinearOpBase<double> > zero 
00146             = Thyra::zero(tA_sub->range(),tA_sub->domain());
00147       
00148       tA = Thyra::block2x2(tA_sub,zero,zero,tA_sub);
00149 
00150       // build blocked vector
00151       RCP<const Thyra::VectorBase<double> > tb_sub 
00152             = Thyra::create_Vector( b, tA_sub->range() );
00153 
00154       RCP<Thyra::VectorBase<double> > tb_m = Thyra::createMember(tA->range());
00155       Thyra::randomize(-1.0,1.0,tb_m.ptr());
00156 
00157       tb = tb_m;
00158    }
00159    RCP<Thyra::VectorBase<double> > tx0;
00160    RCP<Thyra::VectorBase<double> > tx1;
00161    {
00162       tx0 = Thyra::createMember(tA->domain());
00163       tx1 = Thyra::createMember(tA->domain());
00164      
00165       Thyra::randomize(-1.0,1.0,tx0.ptr());
00166       Thyra::randomize(-1.0,1.0,tx1.ptr());
00167    }
00168 
00169    // now comes Stratimikos
00170    RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("BelosGCRODRTest.xml");
00171    Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00172    linearSolverBuilder.setParameterList(paramList);
00173  
00174    // Create a linear solver factory given information read from the
00175    // parameter list.
00176    RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
00177          linearSolverBuilder.createLinearSolveStrategy("");
00178 
00179    // Create a linear solver based on the forward operator A
00180    RCP<Thyra::LinearOpWithSolveBase<double> > lows =
00181          Thyra::linearOpWithSolve(*lowsFactory, tA);
00182 
00183    // Solve the linear system 
00184    Thyra::SolveStatus<double> status; 
00185    status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx0.ptr());
00186    status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx1.ptr());
00187 }
 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