Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_InterlacedOpUnitTest.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Stokhos Package
6 // Copyright (2009) 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 Eric T. Phipps (etphipp@sandia.gov).
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 #include <Teuchos_ConfigDefs.hpp>
45 #include <Teuchos_UnitTestHarness.hpp>
46 #include <Teuchos_TimeMonitor.hpp>
47 #include <Teuchos_RCP.hpp>
48 
50 
51 // Stokhos Stochastic Galerkin
52 #include "Stokhos_Epetra.hpp"
55 
56 #ifdef HAVE_MPI
57 #include "Epetra_MpiComm.h"
58 #else
59 #include "Epetra_SerialComm.h"
60 #endif
61 #include "Epetra_CrsGraph.h"
62 #include "Epetra_Map.h"
63 #include "EpetraExt_BlockUtility.h"
64 #include "EpetraExt_RowMatrixOut.h"
65 
66 TEUCHOS_UNIT_TEST(interlaced_op, test)
67 {
68 #ifdef HAVE_MPI
69  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
70 #else
71  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
72 #endif
73 
74  //int rank = comm->MyPID();
75  int numProc = comm->NumProc();
76 
77  int num_KL = 1;
78  int porder = 5;
79  bool full_expansion = false;
80 
81  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
82  Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
83  Teuchos::RCP<Stokhos::ParallelData> sg_parallel_data;
84  Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expansion;
85  {
86  if(full_expansion)
87  Cijk = basis->computeTripleProductTensor();
88  else
89  Cijk = basis->computeLinearTripleProductTensor();
90 
91  Teuchos::ParameterList parallelParams;
92  parallelParams.set("Number of Spatial Processors", numProc);
93  sg_parallel_data = Teuchos::rcp(new Stokhos::ParallelData(basis, Cijk, comm,
94  parallelParams));
95 
96  expansion = Teuchos::rcp(new Stokhos::AlgebraicOrthogPolyExpansion<int,double>(basis,
97  Cijk));
98  }
99  Teuchos::RCP<const EpetraExt::MultiComm> sg_comm = sg_parallel_data->getMultiComm();
100 
101  // determinstic PDE graph
102  Teuchos::RCP<Epetra_Map> determRowMap = Teuchos::rcp(new Epetra_Map(-1,10,0,*comm));
103  Teuchos::RCP<Epetra_CrsGraph> determGraph = Teuchos::rcp(new Epetra_CrsGraph(Copy,*determRowMap,1));
104  for(int row=0;row<determRowMap->NumMyElements();row++) {
105  int gid = determRowMap->GID(row);
106  determGraph->InsertGlobalIndices(gid,1,&gid);
107  }
108  for(int row=1;row<determRowMap->NumMyElements()-1;row++) {
109  int gid = determRowMap->GID(row);
110  int indices[2] = {gid-1,gid+1};
111  determGraph->InsertGlobalIndices(gid,2,indices);
112  }
113  determGraph->FillComplete();
114 
115  Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList);
116  params->set("Scale Operator by Inverse Basis Norms", false);
117  params->set("Include Mean", true);
118  params->set("Only Use Linear Terms", false);
119 
120  Teuchos::RCP<Stokhos::EpetraSparse3Tensor> epetraCijk =
121  Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(basis,Cijk,sg_comm));
122  Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> W_sg_blocks =
123  Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(basis, epetraCijk->getStochasticRowMap(), determRowMap, determRowMap, sg_comm));
124  for(int i=0; i<W_sg_blocks->size(); i++) {
125  Teuchos::RCP<Epetra_CrsMatrix> crsMat = Teuchos::rcp(new Epetra_CrsMatrix(Copy,*determGraph));
126  crsMat->PutScalar(1.0 + i);
127  W_sg_blocks->setCoeffPtr(i,crsMat); // allocate a bunch of matrices
128  }
129 
130  Teuchos::RCP<const Epetra_Map> sg_map =
131  Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
132  *determRowMap, *(epetraCijk->getStochasticRowMap()),
133  *(epetraCijk->getMultiComm())));
134 
135  // build an interlaced operator (object under test) and a benchmark
136  // fully assembled operator
138 
139  Stokhos::InterlacedOperator op(sg_comm,basis,epetraCijk,determGraph,params);
140  op.PutScalar(0.0);
141  op.setupOperator(W_sg_blocks);
142 
143  Stokhos::FullyAssembledOperator full_op(sg_comm,basis,epetraCijk,determGraph,sg_map,sg_map,params);
144  full_op.PutScalar(0.0);
145  full_op.setupOperator(W_sg_blocks);
146 
147  // here we test interlaced operator against the fully assembled operator
149  bool result = true;
150  for(int i=0;i<100;i++) {
151  // build vector for fully assembled operator (blockwise)
152  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> x_vec_blocks =
153  Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(basis,epetraCijk->getStochasticRowMap(),determRowMap,epetraCijk->getMultiComm()));
154  Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f_vec_blocks =
155  Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(basis,epetraCijk->getStochasticRowMap(),determRowMap,epetraCijk->getMultiComm()));
156  Teuchos::RCP<Epetra_Vector> x_vec_blocked = x_vec_blocks->getBlockVector();
157  Teuchos::RCP<Epetra_Vector> f_vec_blocked = f_vec_blocks->getBlockVector();
158  x_vec_blocked->Random(); // build an initial vector
159  f_vec_blocked->PutScalar(0.0);
160 
161  // build interlaced vectors
162  Teuchos::RCP<Epetra_Vector> x_vec_inter = Teuchos::rcp(new Epetra_Vector(op.OperatorDomainMap()));
163  Teuchos::RCP<Epetra_Vector> f_vec_inter = Teuchos::rcp(new Epetra_Vector(op.OperatorRangeMap()));
164  Teuchos::RCP<Epetra_Vector> f_vec_blk_inter = Teuchos::rcp(new Epetra_Vector(op.OperatorRangeMap()));
165  Stokhos::SGModelEvaluator_Interlaced::copyToInterlacedVector(*x_vec_blocks,*x_vec_inter); // copy random x to
166  f_vec_inter->PutScalar(0.0);
167 
168  full_op.Apply(*x_vec_blocked,*f_vec_blocked);
169  op.Apply(*x_vec_inter,*f_vec_inter);
170 
171  // copy blocked action to interlaced for comparison
172  Stokhos::SGModelEvaluator_Interlaced::copyToInterlacedVector(*f_vec_blocks,*f_vec_blk_inter);
173 
174  // compute norm
175  double error = 0.0;
176  double true_norm = 0.0;
177  f_vec_blk_inter->NormInf(&true_norm);
178  f_vec_blk_inter->Update(-1.0,*f_vec_inter,1.0);
179  f_vec_blk_inter->NormInf(&error);
180 
181  out << "rel error = " << error/true_norm << " ( " << true_norm << " ), ";
182  result &= (error/true_norm < 1e-14);
183  }
184  out << std::endl;
185 
186  TEST_ASSERT(result);
187 }
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Copy
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
TEUCHOS_UNIT_TEST(interlaced_op, test)
Orthogonal polynomial expansions limited to algebraic operations.
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Wrap Apply() to add a timer.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< int, double > > buildBasis(int num_KL, int porder)
An Epetra operator representing the block stochastic Galerkin operator generated by fully assembling ...
static void copyToInterlacedVector(const Stokhos::EpetraVectorOrthogPoly &x_sg, Epetra_Vector &x)