Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_AdaptivityToolsUnitTest.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 
49 #include "Stokhos_InterlacedTestSupport.hpp"
50 
51 // Stokhos Stochastic Galerkin
52 #include "Stokhos_BasisInteractionGraph.hpp"
53 #include "Stokhos_EpetraSparse3Tensor.hpp"
54 #include "Stokhos_ParallelData.hpp"
55 #include "Stokhos_AdaptivityManager.hpp"
56 #include "Stokhos_AdaptivityUtils.hpp"
57 
58 #ifdef HAVE_MPI
59  #include "Epetra_MpiComm.h"
60  #include "EpetraExt_MultiMpiComm.h"
61  #include "mpi.h"
62 #else
63  #include "Epetra_SerialComm.h"
64  #include "EpetraExt_MultiSerialComm.h"
65 #endif
66 #include "Epetra_CrsGraph.h"
67 
68 #include "EpetraExt_RowMatrixOut.h"
69 
70 Teuchos::RCP<Epetra_CrsGraph> buildTridiagonalGraph(int numUnks,const Epetra_Comm & Comm)
71 {
72  Epetra_Map map(-1,numUnks,0,Comm);
73  Teuchos::RCP<Epetra_CrsGraph> graph
74  = Teuchos::rcp(new Epetra_CrsGraph(Copy,map,0));
75 
76  // build tridiagonal graph
77  int colCnt = 3;
78  int * colPtr = 0;
79  int colIndices[3];
80  int colOffset[] = {-1, 0, 1};
81  for(int myRow=0;myRow<numUnks;myRow++) {
82  int row = map.GID(myRow);
83  for(int i=0;i<3;i++)
84  colIndices[i] = colOffset[i]+row;
85  colCnt = 3;
86  colPtr = colIndices;
87 
88  if(row==0) {
89  colCnt = 2;
90  colPtr = colIndices+1;
91  }
92  else if(row==map.NumGlobalElements()-1)
93  colCnt = 2;
94 
95  TEUCHOS_ASSERT(graph->InsertGlobalIndices(row,colCnt,colPtr)==0);
96  }
97 
98  graph->FillComplete();
99 
100  return graph;
101 }
102 
103 Teuchos::RCP<Epetra_CrsMatrix> buildTridiagonalOp(const Epetra_CrsGraph & graph,double * stencil)
104 {
105  Teuchos::RCP<Epetra_CrsMatrix> result = Teuchos::rcp(new Epetra_CrsMatrix(Copy,graph));
106  const Epetra_Map & map = result->RowMap();
107 
108  result->PutScalar(0.0);
109 
110  // build tridiagonal graph
111  int colCnt = 3;
112  int * colPtr = 0;
113  int colIndices[3];
114  int colOffset[] = {-1, 0, 1};
115  double * stencilPtr = 0;
116  for(int myRow=0;myRow<map.NumMyElements();myRow++) {
117  int row = map.GID(myRow);
118  for(int i=0;i<3;i++)
119  colIndices[i] = colOffset[i]+row;
120  colCnt = 3;
121  colPtr = colIndices;
122  stencilPtr = stencil;
123 
124  if(row==0) {
125  colCnt = 2;
126  colPtr = colIndices+1;
127  stencilPtr = stencil+1;
128  }
129  else if(row==map.NumGlobalElements()-1)
130  colCnt = 2;
131 
132  TEUCHOS_ASSERT(result->SumIntoGlobalValues(row,colCnt,stencilPtr,colPtr)==0);
133  }
134 
135  return result;
136 }
137 
138 // some testing infrastructure
139 template <typename ScalarT,typename OrdinalT>
140 bool ord_func(std::pair<ScalarT,OrdinalT> const & a, std::pair<ScalarT,OrdinalT> const & b)
141 { return a.first < b.first; }
142 
143 template <typename ScalarT>
144 void generate_sorted_order(const std::vector<ScalarT> & values,std::vector<std::size_t> & order)
145 {
146  typedef std::pair<ScalarT,std::size_t> Pair;
147 
148  // build vector to sort
149  std::vector<Pair> pairValues(values.size());
150  for(std::size_t i=0;i<values.size();i++)
151  pairValues[i] = std::make_pair(values[i],i);
152 
153  // build sorted ordering
154  std::sort(pairValues.begin(),pairValues.end(),ord_func<ScalarT,std::size_t>);
155 
156  // write out new order
157  order.resize(pairValues.size());
158  for(std::size_t i=0;i<pairValues.size();i++)
159  order[i] = pairValues[i].second;
160 }
161 
162 template <typename ScalarT>
163 void apply_ordering(std::vector<ScalarT> & values, const std::vector<std::size_t> & order)
164 {
165  typedef std::pair<std::size_t,ScalarT> Pair;
166 
167  // build vector to sort
168  std::vector<Pair> pairValues(values.size());
169  for(std::size_t i=0;i<order.size();i++)
170  pairValues[i] = std::make_pair(order[i],values[i]);
171 
172  // build sorted ordering
173  std::sort(pairValues.begin(),pairValues.end(),ord_func<std::size_t,ScalarT>);
174 
175  // write out new values
176  for(std::size_t i=0;i<pairValues.size();i++)
177  values[i] = pairValues[i].second;
178 }
179 
180 // Test construction of Linear Algebra (LA) data structures
181 TEUCHOS_UNIT_TEST(tAdaptivityManager, test_interface)
182 {
183  #ifdef HAVE_MPI
184  Teuchos::RCP<const Epetra_MpiComm> mpiComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
185  Teuchos::RCP<const EpetraExt::MultiComm> multiComm = Teuchos::rcp(new EpetraExt::MultiMpiComm(*mpiComm,-1));
186  Teuchos::RCP<const Epetra_Comm> comm = mpiComm;
187  #else
188  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
189  Teuchos::RCP<const EpetraExt::MultiComm> multiComm = Teuchos::rcp(new EpetraExt::MultiSerialComm(-1));
190  #endif
191 
192  int num_KL = 3;
193  int porder = 3;
194 
195  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
196  Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk = basis->computeTripleProductTensor();
197 
198  std::vector<int> order(3);
199  order[0] = 2; order[1] = 3; order[2] = 1;
200  std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > sa_BasisPerDRow(3);
201  sa_BasisPerDRow[0] = buildBasis(num_KL,1);
202  sa_BasisPerDRow[1] = buildBasis(num_KL,1);
203  sa_BasisPerDRow[2] = buildBasis(num_KL,order);
204 
205  Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList);
206 
207  Teuchos::RCP<Epetra_CrsGraph> determGraph = buildTridiagonalGraph(3,*comm);
208  Stokhos::AdaptivityManager adaptMngr(basis,sa_BasisPerDRow,*determGraph,false,-1);
209 
210  TEST_EQUALITY(adaptMngr.getGlobalRowId(0,0),0);
211  TEST_EQUALITY(adaptMngr.getGlobalRowId(1,0),sa_BasisPerDRow[0]->size());
212 
213  TEST_EQUALITY(adaptMngr.getGlobalColId(0,0),0);
214  TEST_EQUALITY(adaptMngr.getGlobalColId(1,0),sa_BasisPerDRow[0]->size());
215 
216  TEST_EQUALITY(adaptMngr.getRowStochasticBasisSize(0),sa_BasisPerDRow[0]->size());
217  TEST_EQUALITY(adaptMngr.getRowStochasticBasisSize(1),sa_BasisPerDRow[1]->size());
218  TEST_EQUALITY(adaptMngr.getRowStochasticBasisSize(2),sa_BasisPerDRow[2]->size());
219 
220  TEST_EQUALITY(adaptMngr.getColStochasticBasisSize(0),sa_BasisPerDRow[0]->size());
221  TEST_EQUALITY(adaptMngr.getColStochasticBasisSize(1),sa_BasisPerDRow[1]->size());
222  TEST_EQUALITY(adaptMngr.getColStochasticBasisSize(2),sa_BasisPerDRow[2]->size());
223 }
224 
225 // Test construction of Linear Algebra (LA) data structures
226 TEUCHOS_UNIT_TEST(tAdaptivityManager, sum_in_op_eq_order)
227 {
228  #ifdef HAVE_MPI
229  Teuchos::RCP<const Epetra_MpiComm> mpiComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
230  Teuchos::RCP<const EpetraExt::MultiComm> multiComm = Teuchos::rcp(new EpetraExt::MultiMpiComm(*mpiComm,-1));
231  Teuchos::RCP<const Epetra_Comm> comm = mpiComm;
232  #else
233  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
234  Teuchos::RCP<const EpetraExt::MultiComm> multiComm = Teuchos::rcp(new EpetraExt::MultiSerialComm(-1));
235  #endif
236 
237  int num_KL = 2;
238  int porder = 2;
239  int numDetermUnks = 3;
240 
241  double stencil[] = {-1.0,2.0,-1.0};
242  Teuchos::RCP<Epetra_CrsGraph> determGraph = buildTridiagonalGraph(numDetermUnks,*comm);
243  Teuchos::RCP<Epetra_CrsMatrix> determOp = buildTridiagonalOp(*determGraph,stencil);
244 
245  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
246  Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk = basis->computeTripleProductTensor();
247  Teuchos::RCP<Stokhos::EpetraSparse3Tensor> epetraCijk =
248  Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(basis,Cijk,multiComm));
249 
250  std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > sa_BasisPerDRow(numDetermUnks);
251  sa_BasisPerDRow[0] = buildBasis(num_KL,1);
252  sa_BasisPerDRow[1] = buildBasis(num_KL,1);
253  sa_BasisPerDRow[2] = buildBasis(num_KL,1);
254 
255  {
256  Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList);
257  params->set("Scale Operator by Inverse Basis Norms",true);
258  Stokhos::AdaptivityManager adaptMngr(basis,sa_BasisPerDRow,*determGraph,false,-1,true);
259 
260  Teuchos::RCP<Epetra_CrsMatrix> sg_A = adaptMngr.buildMatrixFromGraph();
261 
262  TEST_EQUALITY(sg_A->NumGlobalRows(),sg_A->NumGlobalCols()); // check for square
263  TEST_EQUALITY(sg_A->NumGlobalRows(),3*numDetermUnks); // check sizing
264 
265  adaptMngr.sumInOperator(*sg_A,*Cijk,0,*determOp);
266 
267  int determDof = 1;
268  // choose second determrow
269  for(int stochBasis=0;stochBasis<sa_BasisPerDRow[determDof]->size();stochBasis++) {
270  int numEntries = 0;
271  std::vector<std::size_t> order;
272  std::vector<int> stor_indices(9), indices;
273  std::vector<double> stor_values(9), values;
274 
275  sg_A->ExtractGlobalRowCopy(adaptMngr.getGlobalRowId(determDof,stochBasis),9,numEntries,&stor_values[0],&stor_indices[0]);
276 
277  TEST_EQUALITY(numEntries,9);
278 
279  // perform copy so things are correct size
280  for(int i=0;i<numEntries;i++) {
281  indices.push_back(stor_indices[i]);
282  values.push_back(stor_values[i]);
283  }
284 
285  // order indices and values
286  generate_sorted_order(indices,order);
287  apply_ordering(indices,order);
288  apply_ordering(values,order);
289 
290  int rowTerm = basis->index(sa_BasisPerDRow[determDof]->term(stochBasis));
291  int colTerm0 = basis->index(adaptMngr.getColStochasticBasis(determDof)->term(0));
292  int colTerm1 = basis->index(adaptMngr.getColStochasticBasis(determDof)->term(1));
293  int colTerm2 = basis->index(adaptMngr.getColStochasticBasis(determDof)->term(2));
294  double normValue = basis->norm_squared(rowTerm);
295 
296  // test middle row values
297  TEST_EQUALITY(values[0],stencil[0]*Cijk->getValue(rowTerm,colTerm0,0)/normValue);
298  TEST_EQUALITY(values[1],stencil[0]*Cijk->getValue(rowTerm,colTerm1,0)/normValue);
299  TEST_EQUALITY(values[2],stencil[0]*Cijk->getValue(rowTerm,colTerm2,0)/normValue);
300 
301  TEST_EQUALITY(values[3],stencil[1]*Cijk->getValue(rowTerm,colTerm0,0)/normValue);
302  TEST_EQUALITY(values[4],stencil[1]*Cijk->getValue(rowTerm,colTerm1,0)/normValue);
303  TEST_EQUALITY(values[5],stencil[1]*Cijk->getValue(rowTerm,colTerm2,0)/normValue);
304 
305  TEST_EQUALITY(values[6],stencil[2]*Cijk->getValue(rowTerm,colTerm0,0)/normValue);
306  TEST_EQUALITY(values[7],stencil[2]*Cijk->getValue(rowTerm,colTerm1,0)/normValue);
307  TEST_EQUALITY(values[8],stencil[2]*Cijk->getValue(rowTerm,colTerm2,0)/normValue);
308  }
309 
310  TEST_ASSERT(sg_A->Filled());
311  }
312 
313  {
314  Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList);
315  params->set("Scale Operator by Inverse Basis Norms",false);
316  Stokhos::AdaptivityManager adaptMngr(basis,sa_BasisPerDRow,*determGraph,false,-1,false);
317 
318  Teuchos::RCP<Epetra_CrsMatrix> sg_A = adaptMngr.buildMatrixFromGraph();
319 
320  TEST_EQUALITY(sg_A->NumGlobalRows(),sg_A->NumGlobalCols()); // check for square
321  TEST_EQUALITY(sg_A->NumGlobalRows(),3*numDetermUnks); // check sizing
322 
323  adaptMngr.sumInOperator(*sg_A,*Cijk,0,*determOp);
324 
325  int determDof = 1;
326  // choose second determrow
327  for(int stochBasis=0;stochBasis<sa_BasisPerDRow[determDof]->size();stochBasis++) {
328  int numEntries = 0;
329  std::vector<std::size_t> order;
330  std::vector<int> stor_indices(9), indices;
331  std::vector<double> stor_values(9), values;
332 
333  sg_A->ExtractGlobalRowCopy(adaptMngr.getGlobalRowId(determDof,stochBasis),9,numEntries,&stor_values[0],&stor_indices[0]);
334 
335  TEST_EQUALITY(numEntries,9);
336 
337  // perform copy so things are correct size
338  for(int i=0;i<numEntries;i++) {
339  indices.push_back(stor_indices[i]);
340  values.push_back(stor_values[i]);
341  }
342 
343  // order indices and values
344  generate_sorted_order(indices,order);
345  apply_ordering(indices,order);
346  apply_ordering(values,order);
347 
348  int rowTerm = basis->index(sa_BasisPerDRow[determDof]->term(stochBasis));
349  int colTerm0 = basis->index(adaptMngr.getColStochasticBasis(determDof)->term(0));
350  int colTerm1 = basis->index(adaptMngr.getColStochasticBasis(determDof)->term(1));
351  int colTerm2 = basis->index(adaptMngr.getColStochasticBasis(determDof)->term(2));
352 
353  // test middle row values
354  TEST_EQUALITY(values[0],stencil[0]*Cijk->getValue(rowTerm,colTerm0,0));
355  TEST_EQUALITY(values[1],stencil[0]*Cijk->getValue(rowTerm,colTerm1,0));
356  TEST_EQUALITY(values[2],stencil[0]*Cijk->getValue(rowTerm,colTerm2,0));
357 
358  TEST_EQUALITY(values[3],stencil[1]*Cijk->getValue(rowTerm,colTerm0,0));
359  TEST_EQUALITY(values[4],stencil[1]*Cijk->getValue(rowTerm,colTerm1,0));
360  TEST_EQUALITY(values[5],stencil[1]*Cijk->getValue(rowTerm,colTerm2,0));
361 
362  TEST_EQUALITY(values[6],stencil[2]*Cijk->getValue(rowTerm,colTerm0,0));
363  TEST_EQUALITY(values[7],stencil[2]*Cijk->getValue(rowTerm,colTerm1,0));
364  TEST_EQUALITY(values[8],stencil[2]*Cijk->getValue(rowTerm,colTerm2,0));
365  }
366 
367  TEST_ASSERT(sg_A->Filled());
368  }
369 }
370 
371 TEUCHOS_UNIT_TEST(tAdaptivityManager, sum_in_op_var_order)
372 {
373  #ifdef HAVE_MPI
374  Teuchos::RCP<const Epetra_MpiComm> mpiComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
375  Teuchos::RCP<const EpetraExt::MultiComm> multiComm = Teuchos::rcp(new EpetraExt::MultiMpiComm(*mpiComm,-1));
376  Teuchos::RCP<const Epetra_Comm> comm = mpiComm;
377  #else
378  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
379  Teuchos::RCP<const EpetraExt::MultiComm> multiComm = Teuchos::rcp(new EpetraExt::MultiSerialComm(-1));
380  #endif
381 
382  int num_KL = 4;
383  int porder = 3;
384  int numDetermUnks = 3;
385 
386  double stencil[] = {-1.0,2.0,-1.0};
387  Teuchos::RCP<Epetra_CrsGraph> determGraph = buildTridiagonalGraph(numDetermUnks,*comm);
388  Teuchos::RCP<Epetra_CrsMatrix> determOp = buildTridiagonalOp(*determGraph,stencil);
389 
390  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
391  Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk = basis->computeTripleProductTensor();
392  Teuchos::RCP<Stokhos::EpetraSparse3Tensor> epetraCijk =
393  Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(basis,Cijk,multiComm));
394 
395  Stokhos::BasisInteractionGraph big(*basis);
396 
397  std::vector<int> vorder(4);
398  vorder[0] = 2; vorder[1] = 3; vorder[2] = 2; vorder[3] = 0;
399  std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > sa_BasisPerDRow(numDetermUnks);
400  sa_BasisPerDRow[0] = buildBasis(num_KL,1);
401  sa_BasisPerDRow[1] = buildBasis(num_KL,2);
402  sa_BasisPerDRow[2] = buildBasis(num_KL,vorder);
403 
404  {
405  Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(new Teuchos::ParameterList);
406  params->set("Scale Operator by Inverse Basis Norms",false);
407  Stokhos::AdaptivityManager adaptMngr(basis,sa_BasisPerDRow,*determGraph,false,-1,false);
408 
409  Teuchos::RCP<Epetra_CrsMatrix> sg_A = adaptMngr.buildMatrixFromGraph();
410 
411  TEST_EQUALITY(sg_A->NumGlobalRows(),sg_A->NumGlobalCols()); // check for square
412 
413  adaptMngr.sumInOperator(*sg_A,*Cijk,0,*determOp);
414 
415  out << "Summed into" << std::endl;
416 
417  int determDof = 1;
418  // choose second determrow
419  for(int stochBasis=0;stochBasis<sa_BasisPerDRow[determDof]->size();stochBasis++) {
420  int numEntries = 0;
421  std::vector<std::size_t> order;
422  std::vector<int> stor_indices(400), indices;
423  std::vector<double> stor_values(400), values;
424 
425  out << "grabbing row " << stochBasis << " values" << std::endl;
426  TEST_ASSERT(sg_A->ExtractGlobalRowCopy(adaptMngr.getGlobalRowId(determDof,stochBasis),400,numEntries,&stor_values[0],&stor_indices[0])==0);
427  out << "num entries " << numEntries << std::endl;
428 
429  TEST_ASSERT(numEntries<400);
430 
431  // perform copy so things are correct size
432  for(int i=0;i<numEntries;i++) {
433  indices.push_back(stor_indices[i]);
434  values.push_back(stor_values[i]);
435  }
436 
437  // order indices and values
438  out << "sort row" << std::endl;
439  generate_sorted_order(indices,order);
440  apply_ordering(indices,order);
441  apply_ordering(values,order);
442 
443  out << "grabbing row index, and row norm" << std::endl;
444  int rowTerm = basis->index(sa_BasisPerDRow[determDof]->term(stochBasis));
445 
446  out << "checking matrix" << std::endl;
447  // test middle row values
448  int offset = 0;
449  for(int stochColBasisIndex = 0;stochColBasisIndex<3;stochColBasisIndex++) {
450  for(int stochCol=0;stochCol<adaptMngr.getColStochasticBasisSize(stochColBasisIndex);stochCol++) {
451  int colTerm = basis->index(adaptMngr.getColStochasticBasis(stochColBasisIndex)->term(stochCol));
452 
453  if(big(rowTerm,colTerm)) {
454  TEST_EQUALITY(indices[offset],adaptMngr.getGlobalColId(stochColBasisIndex,stochCol));
455  TEST_EQUALITY(values[offset],stencil[stochColBasisIndex]*Cijk->getValue(rowTerm,colTerm,0));
456  offset++;
457  }
458  }
459  out << "offset = " << offset << std::endl;
460  }
461  }
462 
463  TEST_ASSERT(sg_A->Filled());
464  }
465 }
466 
467 // Test construction of Linear Algebra (LA) data structures
468 TEUCHOS_UNIT_TEST(tBuildAdaptLA, test_uniform)
469 {
470  #ifdef HAVE_MPI
471  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
472  #else
473  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
474  #endif
475 
476  int num_KL = 3;
477  int porder = 3;
478 
479  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
480  Teuchos::RCP<Epetra_CrsGraph> determGraph = buildTridiagonalGraph(10,*comm);
481 
482  std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > sa_BasisPerDRow(determGraph->NumMyRows(),basis);
483 
484  // build adapted row map for uniform grid
485  Teuchos::RCP<Epetra_Map> rowMap;
486  std::vector<int> sa_RowGidOffsets;
487  {
488  int determMyRows = determGraph->RowMap().NumMyElements();
489  int determGlobalRows = determGraph->RowMap().NumGlobalElements();
490 
491  rowMap = Stokhos::adapt_utils::buildAdaptedRowMapAndOffsets(determGraph->Comm(),sa_BasisPerDRow,sa_RowGidOffsets);
492 
493  // tests some sizes
494  TEST_EQUALITY(rowMap->NumMyElements(),determMyRows*basis->size());
495  TEST_EQUALITY(rowMap->NumGlobalElements(),determGlobalRows*basis->size());
496  TEST_EQUALITY(int(sa_RowGidOffsets.size()),determMyRows);
497  TEST_ASSERT(rowMap->LinearMap());
498 
499  // test some content
500  {
501  bool result = true;
502  for(std::size_t i=0;i<sa_RowGidOffsets.size();i++)
503  result &= (sa_RowGidOffsets[i]==rowMap->GID(i*basis->size()));
504  TEST_ASSERT(result);
505  }
506  }
507 
508  // build adapted column map for uniform grid
509  std::vector<int> sa_ColGidOffsets;
510  {
511  Stokhos::adapt_utils::buildAdaptedColOffsets(*determGraph,
512  sa_RowGidOffsets,
513  sa_ColGidOffsets);
514 
515  const Epetra_BlockMap & d_rowMap = determGraph->RowMap();
516  const Epetra_BlockMap & d_colMap = determGraph->ColMap();
517 
518  int determMyCols = d_colMap.NumMyElements();
519 
520  TEST_EQUALITY(int(sa_ColGidOffsets.size()),determMyCols);
521 
522  bool result = true;
523  bool checkOne = false;
524  for(int localColId=0;localColId<determMyCols;localColId++) {
525  int localRowId = d_rowMap.LID(d_colMap.GID(localColId));
526 
527  if(localRowId==-1)
528  continue; // global comlumn not row on this processor
529 
530  checkOne = true;
531  result &= (sa_ColGidOffsets[localColId]==sa_RowGidOffsets[localRowId]);
532  }
533 
534  TEST_ASSERT(checkOne);
535  TEST_ASSERT(result);
536  }
537 }
538 
539 TEUCHOS_UNIT_TEST(tBuildAdaptLA, test_graph)
540 {
541  #ifdef HAVE_MPI
542  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
543  #else
544  Teuchos::RCP<const Epetra_Comm> comm = Teuchos::rcp(new Epetra_SerialComm);
545  #endif
546 
547  int numProc = comm->NumProc();
548  int rank = comm->MyPID();
549 
550  out << "NumProc = " << numProc << ", Rank = " << rank << std::endl;
551 
552  int numDetermRows = 3;
553  int num_KL = 3;
554  int porder = 3;
555 
556  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = buildBasis(num_KL,porder);
557  Teuchos::RCP<Epetra_CrsGraph> determGraph = buildTridiagonalGraph(numDetermRows,*comm);
558 
559  std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > sa_BasisPerDRow(numDetermRows);
560  for(int i=0;i<numDetermRows;i+=3) {
561  sa_BasisPerDRow[i] = buildBasis(num_KL,1);
562  if(i+1<numDetermRows)
563  sa_BasisPerDRow[i+1] = buildBasis(num_KL,2);
564  if(i+2<numDetermRows)
565  sa_BasisPerDRow[i+2] = buildBasis(num_KL,3);
566  }
567 
568  for(int i=0;i<numDetermRows;i++) {
569  out << "Row " << i << ":\n";
570  Stokhos::BasisInteractionGraph(*sa_BasisPerDRow[i]).printGraph(out);
571  out << std::endl;
572  }
573 
574  for(int i=1;i<numDetermRows;i++) {
575  out << "Pair row " << i-1 << ", col " << i << ":\n";
576  Stokhos::BasisInteractionGraph(*basis,*sa_BasisPerDRow[i-1],*sa_BasisPerDRow[i]).printGraph(out);
577  out << std::endl;
578 
579  out << "Pair row " << i << ", col " << i-1 << ":\n";
580  Stokhos::BasisInteractionGraph(*basis,*sa_BasisPerDRow[i],*sa_BasisPerDRow[i-1]).printGraph(out);
581  out << std::endl;
582 
583  }
584 
585  Teuchos::RCP<Epetra_CrsGraph> graph = Stokhos::adapt_utils::buildAdaptedGraph(*determGraph,basis,sa_BasisPerDRow);
586 
587  TEST_ASSERT(graph!=Teuchos::null);
588 
589  // compute expected number of non zeros
590  std::size_t nnz = 0;
591  for(int i=0;i<numDetermRows;i++) {
592  int gid = determGraph->GRID(i);
593  int indices[3];
594  int numRowEntries = 0;
595 
596  determGraph->ExtractGlobalRowCopy(gid,3,numRowEntries,indices);
597 
598  for(int c=0;c<numRowEntries;c++)
599  nnz += Stokhos::BasisInteractionGraph(*basis,*sa_BasisPerDRow[i],*sa_BasisPerDRow[indices[c]]).numNonZeros();
600  }
601 
602  TEST_EQUALITY(graph->NumMyNonzeros(),int(nnz));
603 }
Stokhos_InterlacedTestSupport.hpp
Epetra_BlockMap::NumGlobalElements
int NumGlobalElements() const
Epetra_Map
Stokhos::AdaptivityManager::getColStochasticBasis
Teuchos::RCP< const Stokhos::ProductBasis< int, double > > getColStochasticBasis(int determLid) const
Definition: Stokhos_AdaptivityManager.hpp:174
TEST_ASSERT
#define TEST_ASSERT(v1)
Epetra_CrsGraph::FillComplete
int FillComplete()
Epetra_CrsGraph::Comm
const Epetra_Comm & Comm() const
Epetra_CrsMatrix::ExtractGlobalRowCopy
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Epetra_CrsMatrix::NumGlobalRows
int NumGlobalRows() const
Epetra_CrsMatrix::SumIntoGlobalValues
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Teuchos_ConfigDefs.hpp
Epetra_CrsGraph::NumMyRows
int NumMyRows() const
Epetra_CrsGraph.h
Epetra_CrsGraph::ColMap
const Epetra_BlockMap & ColMap() const
Stokhos::BasisInteractionGraph::numNonZeros
std::size_t numNonZeros() const
How many non zeros are in this graph.
Definition: Stokhos_BasisInteractionGraph.cpp:190
ord_func
bool ord_func(std::pair< ScalarT, OrdinalT > const &a, std::pair< ScalarT, OrdinalT > const &b)
Definition: Stokhos_AdaptivityToolsUnitTest.cpp:140
Epetra_MpiComm
Epetra_CrsGraph::InsertGlobalIndices
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
Epetra_SerialComm.h
Epetra_Comm::MyPID
virtual int MyPID() const =0
Epetra_CrsMatrix::PutScalar
int PutScalar(double ScalarConstant)
Epetra_CrsMatrix::RowMap
const Epetra_Map & RowMap() const
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Stokhos::adapt_utils::buildAdaptedGraph
Teuchos::RCP< Epetra_CrsGraph > buildAdaptedGraph(const Epetra_CrsGraph &determGraph, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &masterBasis, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, bool onlyUseLinear=false, int kExpOrder=-1)
Definition: Stokhos_AdaptivityUtils.cpp:176
Stokhos_ParallelData.hpp
Stokhos::AdaptivityManager::getGlobalColId
int getGlobalColId(int determLid, int basisIndex) const
Definition: Stokhos_AdaptivityManager.hpp:97
Stokhos::ProductBasis::term
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const =0
Get orders of each coordinate polynomial given an index i.
Epetra_Comm
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Epetra_BlockMap::GID
int GID(int LID) const
Teuchos_UnitTestHarness.hpp
Epetra_BlockMap
Stokhos::adapt_utils::buildAdaptedRowMapAndOffsets
Teuchos::RCP< Epetra_Map > buildAdaptedRowMapAndOffsets(const Epetra_Comm &Comm, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, std::vector< int > &myRowGidOffsets)
Definition: Stokhos_AdaptivityUtils.cpp:51
Stokhos::AdaptivityManager::getColStochasticBasisSize
int getColStochasticBasisSize(int determLid) const
Definition: Stokhos_AdaptivityManager.hpp:154
Epetra_CrsGraph::RowMap
const Epetra_BlockMap & RowMap() const
generate_sorted_order
void generate_sorted_order(const std::vector< ScalarT > &values, std::vector< std::size_t > &order)
Definition: Stokhos_AdaptivityToolsUnitTest.cpp:144
Epetra_BlockMap::LID
int LID(int GID) const
Epetra_CrsMatrix::NumGlobalCols
int NumGlobalCols() const
buildTridiagonalGraph
Teuchos::RCP< Epetra_CrsGraph > buildTridiagonalGraph(int numUnks, const Epetra_Comm &Comm)
Definition: Stokhos_AdaptivityToolsUnitTest.cpp:70
Stokhos::adapt_utils::buildAdaptedColOffsets
void buildAdaptedColOffsets(const Epetra_CrsGraph &determGraph, const std::vector< int > &myRowGidOffsets, std::vector< int > &myColGidOffsets)
Definition: Stokhos_AdaptivityUtils.cpp:87
Epetra_SerialComm
Epetra_MpiComm.h
Stokhos::AdaptivityManager::getRowStochasticBasisSize
int getRowStochasticBasisSize(int determLid) const
Definition: Stokhos_AdaptivityManager.hpp:148
buildBasis
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< int, double > > buildBasis(int num_KL, int porder)
Definition: Stokhos_InterlacedTestSupport.cpp:111
Stokhos::AdaptivityManager::buildMatrixFromGraph
Teuchos::RCP< Epetra_CrsMatrix > buildMatrixFromGraph() const
Definition: Stokhos_AdaptivityManager.cpp:128
Epetra_CrsGraph::GRID
int GRID(int LRID_in) const
apply_ordering
void apply_ordering(std::vector< ScalarT > &values, const std::vector< std::size_t > &order)
Definition: Stokhos_AdaptivityToolsUnitTest.cpp:163
Stokhos_BasisInteractionGraph.hpp
Epetra_Comm::NumProc
virtual int NumProc() const =0
Epetra_CrsMatrix::Filled
bool Filled() const
Epetra_CrsMatrix
buildTridiagonalOp
Teuchos::RCP< Epetra_CrsMatrix > buildTridiagonalOp(const Epetra_CrsGraph &graph, double *stencil)
Definition: Stokhos_AdaptivityToolsUnitTest.cpp:103
Teuchos_TimeMonitor.hpp
Epetra_CrsGraph::ExtractGlobalRowCopy
int ExtractGlobalRowCopy(int_type Row, int LenOfIndices, int &NumIndices, int_type *Indices) const
Stokhos_AdaptivityManager.hpp
TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Stokhos_EpetraSparse3Tensor.hpp
TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(tAdaptivityManager, test_interface)
Definition: Stokhos_AdaptivityToolsUnitTest.cpp:181
Epetra_CrsGraph::NumMyNonzeros
int NumMyNonzeros() const
Stokhos::AdaptivityManager::sumInOperator
void sumInOperator(Epetra_CrsMatrix &A, const Stokhos::Sparse3Tensor< int, double > &Cijk, int k, const Epetra_CrsMatrix &J_k) const
Definition: Stokhos_AdaptivityManager.cpp:181
Teuchos_RCP.hpp
Stokhos_AdaptivityUtils.hpp
Stokhos::AdaptivityManager::getGlobalRowId
int getGlobalRowId(int determLid, int basisIndex) const
Definition: Stokhos_AdaptivityManager.hpp:111
Epetra_CrsGraph
Generated by   doxygen 1.8.5