Thyra Package Browser (Single Doxygen Collection)
Version of the Day
Loading...
Searching...
No Matches
adapters
epetra
test
EpetraOperatorWrapper
EpetraOperatorWrapper_UnitTests.cpp
Go to the documentation of this file.
1
/*
2
// @HEADER
3
// ***********************************************************************
4
//
5
// Thyra: Interfaces and Support for Abstract Numerical Algorithms
6
// Copyright (2004) 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 (bartlettra@ornl.gov)
39
//
40
// ***********************************************************************
41
// @HEADER
42
*/
43
44
#include "
Thyra_EpetraOperatorWrapper.hpp
"
45
#include "Thyra_VectorStdOps.hpp"
46
#include "
Thyra_EpetraThyraWrappers.hpp
"
47
#include "
Thyra_EpetraLinearOp.hpp
"
48
#include "Thyra_DefaultBlockedLinearOp.hpp"
49
#include "Thyra_ProductVectorBase.hpp"
50
#include "Thyra_SpmdVectorSpaceBase.hpp"
51
#include "Thyra_DetachedSpmdVectorView.hpp"
52
#include "Thyra_TestingTools.hpp"
53
#include "Thyra_LinearOpTester.hpp"
54
#include "Trilinos_Util_CrsMatrixGallery.h"
55
#include "
Teuchos_GlobalMPISession.hpp
"
56
#include "
Teuchos_VerboseObject.hpp
"
57
#include "
Teuchos_XMLParameterListHelpers.hpp
"
58
#include "
Teuchos_CommandLineProcessor.hpp
"
59
#include "
Teuchos_StandardCatchMacros.hpp
"
60
61
#ifdef HAVE_MPI
62
# include "Epetra_MpiComm.h"
63
#else
64
# include "Epetra_SerialComm.h"
65
#endif
66
#include "Epetra_Vector.h"
67
#include "Epetra_CrsMatrix.h"
68
69
#include "
Teuchos_UnitTestHarness.hpp
"
70
71
72
namespace
Thyra
{
73
74
75
using
Teuchos::null;
76
using
Teuchos::RCP
;
77
using
Teuchos::rcp
;
78
using
Teuchos::rcp_dynamic_cast;
79
using
Teuchos::inOutArg;
80
using
Teuchos::as
;
81
82
83
TEUCHOS_UNIT_TEST
(
EpetraOperatorWrapper
,
basic
)
84
{
85
86
#ifdef HAVE_MPI
87
Epetra_MpiComm
comm
(
MPI_COMM_WORLD
);
88
#else
89
Epetra_SerialComm
comm
;
90
#endif
91
92
out <<
"\nRunning on "
<<
comm
.NumProc() <<
" processors\n"
;
93
94
int
nx
= 39;
// essentially random values
95
int
ny
= 53;
96
97
out <<
"Using Trilinos_Util to create test matrices\n"
;
98
99
// create some big blocks to play with
100
Trilinos_Util::CrsMatrixGallery
FGallery
(
"recirc_2d"
,
comm
,
false
);
// CJ TODO FIXME: change for Epetra64
101
FGallery
.Set(
"nx"
,
nx
);
102
FGallery
.Set(
"ny"
,
ny
);
103
RCP<Epetra_CrsMatrix>
F
= rcp(
FGallery
.GetMatrix(),
false
);
104
105
Trilinos_Util::CrsMatrixGallery
CGallery
(
"laplace_2d"
,
comm
,
false
);
// CJ TODO FIXME: change for Epetra64
106
CGallery
.Set(
"nx"
,
nx
);
107
CGallery
.Set(
"ny"
,
ny
);
108
RCP<Epetra_CrsMatrix>
C
= rcp(
CGallery
.GetMatrix(),
false
);
109
110
Trilinos_Util::CrsMatrixGallery
BGallery
(
"diag"
,
comm
,
false
);
// CJ TODO FIXME: change for Epetra64
111
BGallery
.Set(
"nx"
,
nx
*
ny
);
112
BGallery
.Set(
"a"
,5.0);
113
RCP<Epetra_CrsMatrix>
B
= rcp(
BGallery
.GetMatrix(),
false
);
114
115
Trilinos_Util::CrsMatrixGallery
BtGallery
(
"diag"
,
comm
,
false
);
// CJ TODO FIXME: change for Epetra64
116
BtGallery
.Set(
"nx"
,
nx
*
ny
);
117
BtGallery
.Set(
"a"
,3.0);
118
RCP<Epetra_CrsMatrix>
Bt
= rcp(
BtGallery
.GetMatrix(),
false
);
119
120
// load'em up in a thyra operator
121
out <<
"Building block2x2 Thyra matrix ... wrapping in EpetraOperatorWrapper\n"
;
122
const
RCP<const LinearOpBase<double>
>
A
=
123
Thyra::block2x2<double>
(
124
Thyra::epetraLinearOp(
F
),
125
Thyra::epetraLinearOp(
Bt
),
126
Thyra::epetraLinearOp(
B
),
127
Thyra::epetraLinearOp(
C
),
128
"A"
129
);
130
131
const
RCP<Thyra::EpetraOperatorWrapper>
epetra_A
=
132
rcp(
new
Thyra::EpetraOperatorWrapper
(
A
));
133
134
// begin the tests!
135
const
Epetra_Map
&
rangeMap
=
epetra_A
->OperatorRangeMap();
136
const
Epetra_Map
&
domainMap
=
epetra_A
->OperatorDomainMap();
137
138
// check to see that the number of global elements is correct
139
TEST_EQUALITY(
rangeMap
.NumGlobalElements(), 2*
nx
*
ny
);
140
TEST_EQUALITY(
domainMap
.NumGlobalElements(), 2*
nx
*
ny
);
141
142
// largest global ID should be one less then the # of elements
143
TEST_EQUALITY(
rangeMap
.NumGlobalElements()-1,
rangeMap
.MaxAllGID());
144
TEST_EQUALITY(
domainMap
.NumGlobalElements()-1,
domainMap
.MaxAllGID());
145
146
// create a vector to test: copyThyraIntoEpetra
147
{
148
const
RCP<VectorBase<double>
>
tv
= Thyra::createMember(
A
->domain());
149
Thyra::randomize
(-100.0, 100.0,
tv
.
ptr
());
150
const
RCP<const VectorBase<double>
>
tv_0
=
151
Thyra::productVectorBase<double>
(
tv
)->getVectorBlock(0);
152
const
RCP<const VectorBase<double>
>
tv_1
=
153
Thyra::productVectorBase<double>
(
tv
)->getVectorBlock(1);
154
const
Thyra::ConstDetachedSpmdVectorView<double>
vv_0
(
tv_0
);
155
const
Thyra::ConstDetachedSpmdVectorView<double>
vv_1
(
tv_1
);
156
157
int
off_0
=
vv_0
.globalOffset();
158
int
off_1
=
vv_1
.globalOffset();
159
160
// create its Epetra counter part
161
Epetra_Vector
ev
(
epetra_A
->OperatorDomainMap());
162
epetra_A
->copyThyraIntoEpetra(*
tv
,
ev
);
163
164
// compare handle_tv to ev!
165
TEST_EQUALITY(
tv
->space()->dim(),
as<Ordinal>
(
ev
.GlobalLength()));
166
const
int
numMyElements
=
domainMap
.NumMyElements();
167
double
tval
= 0.0;
168
for
(
int
i
=0;
i
<
numMyElements
;
i
++) {
169
int
gid
=
domainMap
.GID(
i
);
170
if
(
gid
<
nx
*
ny
)
171
tval
=
vv_0
[
gid
-
off_0
];
172
else
173
tval
=
vv_1
[
gid
-
off_1
-
nx
*
ny
];
174
TEST_EQUALITY(
ev
[
i
],
tval
);
175
}
176
}
177
178
// create a vector to test: copyEpetraIntoThyra
179
{
180
// create an Epetra vector
181
Epetra_Vector
ev
(
epetra_A
->OperatorDomainMap());
182
ev
.Random();
183
184
// create its thyra counterpart
185
const
RCP<VectorBase<double>
>
tv
= Thyra::createMember(
A
->domain());
186
const
RCP<const VectorBase<double>
>
tv_0
=
187
Thyra::productVectorBase<double>
(
tv
)->getVectorBlock(0);
188
const
RCP<const VectorBase<double>
>
tv_1
=
189
Thyra::productVectorBase<double>
(
tv
)->getVectorBlock(1);
190
const
Thyra::ConstDetachedSpmdVectorView<double>
vv_0
(
tv_0
);
191
const
Thyra::ConstDetachedSpmdVectorView<double>
vv_1
(
tv_1
);
192
193
int
off_0
=
rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double>
>(
194
tv_0
->space())->localOffset();
195
int
off_1
=
rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double>
>(
196
tv_1
->space())->localOffset();
197
198
epetra_A
->copyEpetraIntoThyra(
ev
,
tv
.
ptr
());
199
200
// compare tv to ev!
201
TEST_EQUALITY(
tv
->space()->dim(),
as<Ordinal>
(
ev
.GlobalLength()));
202
int
numMyElements
=
domainMap
.NumMyElements();
203
double
tval
= 0.0;
204
for
(
int
i
=0;
i
<
numMyElements
;
i
++) {
205
int
gid
=
domainMap
.GID(
i
);
206
if
(
gid
<
nx
*
ny
)
207
tval
=
vv_0
[
gid
-
off_0
];
208
else
209
tval
=
vv_1
[
gid
-
off_1
-
nx
*
ny
];
210
TEST_EQUALITY(
ev
[
i
],
tval
);
211
}
212
}
213
214
// Test using Thyra::LinearOpTester
215
const
RCP<const LinearOpBase<double>
>
thyraEpetraOp
= epetraLinearOp(
epetra_A
);
216
LinearOpTester<double>
linearOpTester
;
217
linearOpTester
.show_all_tests(
true
);
218
const
bool
checkResult
=
linearOpTester
.check(*
thyraEpetraOp
, inOutArg(out));
219
TEST_ASSERT(
checkResult
);
220
221
}
222
223
224
}
// namespace Thyra
TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(Comm, Issue1029)
Teuchos_CommandLineProcessor.hpp
Teuchos_GlobalMPISession.hpp
Teuchos_StandardCatchMacros.hpp
Teuchos_UnitTestHarness.hpp
Teuchos_VerboseObject.hpp
Teuchos_XMLParameterListHelpers.hpp
Thyra_EpetraLinearOp.hpp
Thyra_EpetraOperatorWrapper.hpp
Thyra_EpetraThyraWrappers.hpp
Teuchos::RCP
Teuchos::RCP::ptr
Ptr< T > ptr() const
Thyra::EpetraOperatorWrapper
Implements the Epetra_Operator interface with a Thyra LinearOperator.
Definition
Thyra_EpetraOperatorWrapper.hpp:64
A
B
C
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra
Definition
Thyra_DiagonalEpetraLinearOpWithSolveFactory.cpp:53
Generated by
1.10.0