Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
45#include "Thyra_VectorStdOps.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"
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
70
71
72namespace Thyra {
73
74
75using Teuchos::null;
76using Teuchos::RCP;
77using Teuchos::rcp;
78using Teuchos::rcp_dynamic_cast;
79using Teuchos::inOutArg;
80using Teuchos::as;
81
82
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";
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 =
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(Comm, Issue1029)
Ptr< T > ptr() const
Implements the Epetra_Operator interface with a Thyra LinearOperator.
TypeTo as(const TypeFrom &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEST_ASSERT(castedDep1->getValuesAndValidators().size()==2)
TEST_EQUALITY(rcp_dynamic_cast< const EnhancedNumberValidator< double > >(castedDep1->getValuesAndValidators().find("val1") ->second, true) ->getMax(), double1Vali->getMax())