Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
test_belos_epetra_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>
45#include <Teuchos_UnitTestHarness.hpp>
46#include <Teuchos_RCP.hpp>
47#include <Teuchos_ParameterList.hpp>
48#include "Teuchos_XMLParameterListHelpers.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"
66#include "Thyra_EpetraThyraWrappers.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
79using Teuchos::RCP;
80using Teuchos::rcp;
81
82Teuchos::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
111TEUCHOS_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
121 Teuchos::RCP<Epetra_CrsMatrix> mat = buildMatrix(100,Comm);
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
134 RCP<const Thyra::LinearOpBase<double> >
135 tA = Thyra::epetraLinearOp( mat );
136 RCP<Thyra::VectorBase<double> >
137 tx0 = Thyra::create_Vector( x0, tA->domain() );
138 RCP<Thyra::VectorBase<double> >
139 tx1 = Thyra::create_Vector( x1, tA->domain() );
140 RCP<const Thyra::VectorBase<double> >
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.
150 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
151 linearSolverBuilder.createLinearSolveStrategy("");
152
153 // Create a linear solver based on the forward operator A
154 RCP<Thyra::LinearOpWithSolveBase<double> > lows =
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
163TEUCHOS_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
173 Teuchos::RCP<Epetra_CrsMatrix> mat = buildMatrix(100,Comm);
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
182 RCP<const Thyra::LinearOpBase<double> > tA;
183 RCP<const Thyra::VectorBase<double> > tb;
184 {
185 // build blocked linear Op
186 RCP<const Thyra::LinearOpBase<double> > tA_sub
187 = Thyra::epetraLinearOp( mat );
188 RCP<const Thyra::LinearOpBase<double> > zero
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
194 RCP<const Thyra::VectorBase<double> > tb_sub
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 }
202 RCP<Thyra::VectorBase<double> > tx0;
203 RCP<Thyra::VectorBase<double> > tx1;
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.
219 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
220 linearSolverBuilder.createLinearSolveStrategy("");
221
222 // Create a linear solver based on the forward operator A
223 RCP<Thyra::LinearOpWithSolveBase<double> > lows =
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}
Concrete subclass of Thyra::LinearSolverBuilderBase for creating LinearOpWithSolveFactoryBase objects...
RCP< Thyra::LinearOpWithSolveFactoryBase< Scalar > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
void setParameterList(RCP< ParameterList > const &paramList)
TEUCHOS_UNIT_TEST(belos_gcrodr, multiple_solves)
Teuchos::RCP< Epetra_CrsMatrix > buildMatrix(int nx, Epetra_Comm &comm)