Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test_belos_gcrodr.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Stratimikos: Thyra-based strategies for linear solvers
6 // Copyright (2006) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 #include <Teuchos_ConfigDefs.hpp>
46 #include <Teuchos_RCP.hpp>
49 
50 #include <string>
51 #include <iostream>
52 
53 #ifdef HAVE_MPI
54  #include "Epetra_MpiComm.h"
55 #else
56  #include "Epetra_SerialComm.h"
57 #endif
58 
59 #include "Epetra_Map.h"
60 #include "Epetra_CrsMatrix.h"
61 #include "Epetra_Vector.h"
62 // #include "EpetraExt_RowMatrixOut.h"
63 
64 #include "Thyra_LinearOpBase.hpp"
65 #include "Thyra_VectorBase.hpp"
67 #include "Thyra_EpetraLinearOp.hpp"
68 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
69 #include "Thyra_LinearOpWithSolveHelpers.hpp"
70 #include "Thyra_DefaultZeroLinearOp.hpp"
71 #include "Thyra_DefaultProductVector.hpp"
72 #include "Thyra_DefaultProductVectorSpace.hpp"
73 #include "Thyra_MultiVectorStdOps.hpp"
74 #include "Thyra_VectorStdOps.hpp"
75 #include "Thyra_DefaultBlockedLinearOp.hpp"
76 
78 
79 using Teuchos::RCP;
80 using Teuchos::rcp;
81 
82 Teuchos::RCP<Epetra_CrsMatrix> buildMatrix(int nx, Epetra_Comm & comm)
83 {
84  Epetra_Map map(nx*comm.NumProc(),0,comm);
85  Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(new Epetra_CrsMatrix(Copy,map,3));
86 
87  int offsets[3] = {-1, 0, 1 };
88  double values[3] = { -1, 2, -1};
89  int maxGid = map.MaxAllGID();
90  for(int lid=0;lid<nx;lid++) {
91  int gid = mat->GRID(lid);
92  int numEntries = 3, offset = 0;
93  int indices[3] = { gid+offsets[0],
94  gid+offsets[1],
95  gid+offsets[2] };
96  if(gid==0) { // left end point
97  numEntries = 2;
98  offset = 1;
99  } // right end point
100  else if(gid==maxGid)
101  numEntries = 2;
102 
103  // insert rows
104  mat->InsertGlobalValues(gid,numEntries,values+offset,indices+offset);
105  }
106 
107  mat->FillComplete();
108  return mat;
109 }
110 
111 TEUCHOS_UNIT_TEST(belos_gcrodr, multiple_solves)
112 {
113  // build global (or serial communicator)
114  #ifdef HAVE_MPI
115  Epetra_MpiComm Comm(MPI_COMM_WORLD);
116  #else
117  Epetra_SerialComm Comm;
118  #endif
119 
120  // build and allocate linear system
122  Teuchos::RCP<Epetra_Vector> x0 = rcp(new Epetra_Vector(mat->OperatorDomainMap()));
123  Teuchos::RCP<Epetra_Vector> x1 = rcp(new Epetra_Vector(mat->OperatorDomainMap()));
124  Teuchos::RCP<Epetra_Vector> b = rcp(new Epetra_Vector(mat->OperatorRangeMap()));
125 
126  x0->Random();
127  x1->Random();
128  b->PutScalar(0.0);
129 
130  // sanity check
131  // EpetraExt::RowMatrixToMatrixMarketFile("mat_output.mm",*mat);
132 
133  // build Thyra wrappers
135  tA = Thyra::epetraLinearOp( mat );
137  tx0 = Thyra::create_Vector( x0, tA->domain() );
139  tx1 = Thyra::create_Vector( x1, tA->domain() );
141  tb = Thyra::create_Vector( b, tA->range() );
142 
143  // now comes Stratimikos
144  RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("BelosGCRODRTest.xml");
145  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
146  linearSolverBuilder.setParameterList(paramList);
147 
148  // Create a linear solver factory given information read from the
149  // parameter list.
151  linearSolverBuilder.createLinearSolveStrategy("");
152 
153  // Create a linear solver based on the forward operator A
155  Thyra::linearOpWithSolve(*lowsFactory, tA);
156 
157  // Solve the linear system
158  Thyra::SolveStatus<double> status;
159  status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx0.ptr());
160  status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx1.ptr());
161 }
162 
163 TEUCHOS_UNIT_TEST(belos_gcrodr, 2x2_multiple_solves)
164 {
165  // build global (or serial communicator)
166  #ifdef HAVE_MPI
167  Epetra_MpiComm Comm(MPI_COMM_WORLD);
168  #else
169  Epetra_SerialComm Comm;
170  #endif
171 
172  // build and allocate linear system
174  Teuchos::RCP<Epetra_Vector> b = rcp(new Epetra_Vector(mat->OperatorRangeMap()));
175 
176  b->PutScalar(0.0);
177 
178  // sanity check
179  // EpetraExt::RowMatrixToMatrixMarketFile("mat_output.mm",*mat);
180 
181  // build Thyra wrappers
184  {
185  // build blocked linear Op
187  = Thyra::epetraLinearOp( mat );
189  = Thyra::zero(tA_sub->range(),tA_sub->domain());
190 
191  tA = Thyra::block2x2(tA_sub,zero,zero,tA_sub);
192 
193  // build blocked vector
195  = Thyra::create_Vector( b, tA_sub->range() );
196 
197  RCP<Thyra::VectorBase<double> > tb_m = Thyra::createMember(tA->range());
198  Thyra::randomize(-1.0,1.0,tb_m.ptr());
199 
200  tb = tb_m;
201  }
204  {
205  tx0 = Thyra::createMember(tA->domain());
206  tx1 = Thyra::createMember(tA->domain());
207 
208  Thyra::randomize(-1.0,1.0,tx0.ptr());
209  Thyra::randomize(-1.0,1.0,tx1.ptr());
210  }
211 
212  // now comes Stratimikos
213  RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("BelosGCRODRTest.xml");
214  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
215  linearSolverBuilder.setParameterList(paramList);
216 
217  // Create a linear solver factory given information read from the
218  // parameter list.
220  linearSolverBuilder.createLinearSolveStrategy("");
221 
222  // Create a linear solver based on the forward operator A
224  Thyra::linearOpWithSolve(*lowsFactory, tA);
225 
226  // Solve the linear system
227  Thyra::SolveStatus<double> status;
228  status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx0.ptr());
229  status = Thyra::solve<double>(*lows, Thyra::NOTRANS, *tb, tx1.ptr());
230 }
void setParameterList(RCP< ParameterList > const &paramList)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
RCP< Thyra::LinearOpWithSolveFactoryBase< double > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space=Teuchos::null)
TEUCHOS_UNIT_TEST(belos_gcrodr, multiple_solves)
Ptr< T > ptr() const
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...
Teuchos::RCP< Epetra_CrsMatrix > buildMatrix(int nx, Epetra_Comm &comm)