Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MatrixFreeOperatorUnitTest.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44// Test the exponential of a linear Hermite expansion using an analytic
45// closed form solution
46
47#include "Teuchos_UnitTestHarness.hpp"
48#include "Teuchos_TestingHelpers.hpp"
49#include "Teuchos_UnitTestRepository.hpp"
50#include "Teuchos_GlobalMPISession.hpp"
51
52#include "Stokhos_Epetra.hpp"
53#include "EpetraExt_BlockUtility.h"
55
56#ifdef HAVE_MPI
57#include "Epetra_MpiComm.h"
58#else
59#include "Epetra_SerialComm.h"
60#endif
61
63
65 Teuchos::RCP<const Epetra_Map> sg_x_map;
66 Teuchos::RCP<const Epetra_Map> sg_f_map;
67 Teuchos::RCP<Stokhos::MatrixFreeOperator> mat_free_op;
68 Teuchos::RCP<Stokhos::FullyAssembledOperator> assembled_op;
69 double tol;
70
71 // Can't be a constructor because MPI will not be initialized
72 void setup() {
73
75
76 // Test tolerance
77 tol = 1.0e-12;
78
79 // Basis of dimension 3, order 5
80 const int d = 2;
81 const int p = 3;
82 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
83 for (int i=0; i<d; i++) {
84 bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<int,double>(p));
85 }
86 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
87 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
88
89 // Triple product tensor
90 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
91 basis->computeTripleProductTensor();
92
93 // Create a communicator for Epetra objects
94 Teuchos::RCP<const Epetra_Comm> globalComm;
95#ifdef HAVE_MPI
96 globalComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
97#else
98 globalComm = Teuchos::rcp(new Epetra_SerialComm);
99#endif
100
101 // Create stochastic parallel distribution
102 int num_spatial_procs = -1;
103 int num_procs = globalComm->NumProc();
104 if (num_procs > 1)
105 num_spatial_procs = num_procs / 2;
106 Teuchos::ParameterList parallelParams;
107 parallelParams.set("Number of Spatial Processors", num_spatial_procs);
108 Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data =
109 Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, globalComm,
110 parallelParams));
111 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
112 sg_parallel_data->getMultiComm();
113 Teuchos::RCP<const Epetra_Comm> app_comm =
114 sg_parallel_data->getSpatialComm();
115 Teuchos::RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
116 sg_parallel_data->getEpetraCijk();
117
118 // Deterministic domain map
119 const int num_x = 5;
120 Teuchos::RCP<Epetra_Map> x_map =
121 Teuchos::rcp(new Epetra_Map(num_x, 0, *app_comm));
122
123 // Deterministic column map
124 Teuchos::RCP<Epetra_Map> x_overlap_map =
125 Teuchos::rcp(new Epetra_LocalMap(num_x, 0, *app_comm));
126
127 // Deterministic range map
128 const int num_f = 3;
129 Teuchos::RCP<Epetra_Map> f_map =
130 Teuchos::rcp(new Epetra_Map(num_f, 0, *app_comm));
131
132 // Product domain & range maps
133 Teuchos::RCP<const Epetra_BlockMap> stoch_row_map =
134 epetraCijk->getStochasticRowMap();
135 sg_x_map = Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
136 *x_map, *stoch_row_map, *sg_comm));
137 sg_f_map = Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
138 *f_map, *stoch_row_map, *sg_comm));
139
140 // Deterministic matrix graph
141 const int num_indices = num_x;
142 Teuchos::RCP<Epetra_CrsGraph> graph =
143 Teuchos::rcp(new Epetra_CrsGraph(Copy, *f_map, num_indices));
144 int indices[num_indices];
145 for (int j=0; j<num_indices; j++)
146 indices[j] = x_overlap_map->GID(j);
147 for (int i=0; i<f_map->NumMyElements(); i++)
148 graph->InsertGlobalIndices(f_map->GID(i), num_indices, indices);
149 graph->FillComplete(*x_map, *f_map);
150
151 // Create matrix expansion
152 Teuchos::RCP<Epetra_BlockMap> sg_overlap_map =
153 Teuchos::rcp(new Epetra_LocalMap(
154 basis->size(), 0,
155 *(sg_parallel_data->getStochasticComm())));
156 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > mat_sg =
157 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
158 basis, sg_overlap_map, x_map, f_map, sg_f_map, sg_comm));
159 for (int block=0; block<basis->size(); block++) {
160 Teuchos::RCP<Epetra_CrsMatrix> mat =
161 Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
162 TEUCHOS_TEST_FOR_EXCEPTION(!mat->IndicesAreLocal(), std::logic_error,
163 "Indices are not local!");
164 double values[num_indices];
165 for (int i=0; i<f_map->NumMyElements(); i++) {
166 for (int j=0; j<num_indices; j++) {
167 indices[j] = x_overlap_map->GID(j);
168 values[j] = 0.1*(i+1)*(j+1)*(block+1);
169 }
170 mat->ReplaceMyValues(i, num_indices, values, indices);
171 }
172 mat->FillComplete(*x_map, *f_map);
173 mat_sg->setCoeffPtr(block, mat);
174 }
175
176 // Matrix-free operator
177 Teuchos::RCP<Teuchos::ParameterList> op_params =
178 Teuchos::rcp(new Teuchos::ParameterList);
179 mat_free_op =
180 Teuchos::rcp(new Stokhos::MatrixFreeOperator(
181 sg_comm, basis, epetraCijk, x_map, f_map,
182 sg_x_map, sg_f_map, op_params));
183 mat_free_op->setupOperator(mat_sg);
184
185 // Fully assembled operator
186 assembled_op =
187 Teuchos::rcp(new Stokhos::FullyAssembledOperator(
188 sg_comm, basis, epetraCijk, graph, sg_x_map, sg_f_map,
189 op_params));
190 assembled_op->setupOperator(mat_sg);
191 }
192
193 };
194
196
197 TEUCHOS_UNIT_TEST( Stokhos_MatrixFreeOperator, ApplyUnitTest ) {
198 // Test Apply()
199 Epetra_Vector input(*setup.sg_x_map), result1(*setup.sg_f_map),
200 result2(*setup.sg_f_map), diff(*setup.sg_f_map);
201 input.Random();
202 setup.mat_free_op->Apply(input, result1);
203 setup.assembled_op->Apply(input, result2);
204 diff.Update(1.0, result1, -1.0, result2, 0.0);
205 double nrm;
206 diff.NormInf(&nrm);
207 success = std::fabs(nrm) < setup.tol;
208 out << "Apply infinity norm of difference: " << nrm << std::endl;
209 out << "Matrix-free result = " << std::endl << result1 << std::endl
210 << "Assebled result = " << std::endl << result2 << std::endl;
211 }
212
213 TEUCHOS_UNIT_TEST( Stokhos_MatrixFreeOperator, ApplyTransposeUnitTest ) {
214 // Test tranposed Apply()
215 Epetra_Vector input(*setup.sg_f_map), result1(*setup.sg_x_map),
216 result2(*setup.sg_x_map), diff(*setup.sg_x_map);
217 input.Random();
218 setup.mat_free_op->SetUseTranspose(true);
219 setup.assembled_op->SetUseTranspose(true);
220 setup.mat_free_op->Apply(input, result1);
221 setup.assembled_op->Apply(input, result2);
222 diff.Update(1.0, result1, -1.0, result2, 0.0);
223 double nrm;
224 diff.NormInf(&nrm);
225 success = std::fabs(nrm) < setup.tol;
226 out << "Apply-transpose infinity norm of difference: " << nrm << std::endl;
227 out << "Matrix-free result = " << std::endl << result1 << std::endl
228 << "Assebled result = " << std::endl << result2 << std::endl;
229 }
230
231}
232
233int main( int argc, char* argv[] ) {
234 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
236 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
237}
Copy
int main(int argc, char *argv[])
int NormInf(double *Result) const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
static void SetTracebackMode(int TracebackModeValue)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
Legendre polynomial basis.
An Epetra operator representing the block stochastic Galerkin operator.
TEUCHOS_UNIT_TEST(Stokhos_MatrixFreeOperator, ApplyUnitTest)
Teuchos::RCP< Stokhos::MatrixFreeOperator > mat_free_op
Teuchos::RCP< Stokhos::FullyAssembledOperator > assembled_op