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