Thyra Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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"
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 
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 
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
#define TEST_ASSERT(v1)
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ptr< T > ptr() const
Implements the Epetra_Operator interface with a Thyra LinearOperator.
TypeTo as(const TypeFrom &t)
#define TEST_EQUALITY(v1, v2)