Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
50
51// Stokhos Stochastic Galerkin
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
70Teuchos::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
103Teuchos::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
139template <typename ScalarT,typename OrdinalT>
140bool ord_func(std::pair<ScalarT,OrdinalT> const & a, std::pair<ScalarT,OrdinalT> const & b)
141{ return a.first < b.first; }
142
143template <typename ScalarT>
144void 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
162template <typename ScalarT>
163void 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
181TEUCHOS_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
226TEUCHOS_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
371TEUCHOS_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
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
468TEUCHOS_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 {
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
539TEUCHOS_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}
Copy
Teuchos::RCP< Epetra_CrsGraph > buildTridiagonalGraph(int numUnks, const Epetra_Comm &Comm)
void apply_ordering(std::vector< ScalarT > &values, const std::vector< std::size_t > &order)
bool ord_func(std::pair< ScalarT, OrdinalT > const &a, std::pair< ScalarT, OrdinalT > const &b)
TEUCHOS_UNIT_TEST(tAdaptivityManager, test_interface)
Teuchos::RCP< Epetra_CrsMatrix > buildTridiagonalOp(const Epetra_CrsGraph &graph, double *stencil)
void generate_sorted_order(const std::vector< ScalarT > &values, std::vector< std::size_t > &order)
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< int, double > > buildBasis(int num_KL, int porder)
int LID(int GID) const
int GID(int LID) const
int NumGlobalElements() const
int NumMyElements() const
Teuchos::RCP< Epetra_CrsMatrix > buildMatrixFromGraph() const
Teuchos::RCP< const Stokhos::ProductBasis< int, double > > getColStochasticBasis(int determLid) const
int getColStochasticBasisSize(int determLid) const
int getGlobalColId(int determLid, int basisIndex) const
int getRowStochasticBasisSize(int determLid) const
void sumInOperator(Epetra_CrsMatrix &A, const Stokhos::Sparse3Tensor< int, double > &Cijk, int k, const Epetra_CrsMatrix &J_k) const
int getGlobalRowId(int determLid, int basisIndex) const
std::size_t numNonZeros() const
How many non zeros are in this graph.
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)
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)
void buildAdaptedColOffsets(const Epetra_CrsGraph &determGraph, const std::vector< int > &myRowGidOffsets, std::vector< int > &myColGidOffsets)