Thyra Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraLinearOp_UnitTests.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
4 //
5 // Copyright 2004 NTESS and the Thyra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Thyra_EpetraLinearOp.hpp"
11 #include "Thyra_ScaledLinearOpBase.hpp"
12 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
13 #include "Thyra_RowStatLinearOpBase.hpp"
14 #include "Thyra_LinearOpTester.hpp"
15 #include "Thyra_DefaultBlockedLinearOp.hpp"
16 #include "Thyra_VectorBase.hpp"
17 #include "Thyra_VectorStdOps.hpp"
18 #include "Thyra_MultiVectorBase.hpp"
19 #include "Thyra_MultiVectorStdOps.hpp"
20 #include "Thyra_DefaultProductVectorSpace.hpp"
21 #include "Thyra_DefaultProductVector.hpp"
22 #include "Thyra_DefaultDiagonalLinearOp.hpp"
23 #include "Thyra_DefaultMultipliedLinearOp.hpp"
24 #include "Thyra_DefaultZeroLinearOp.hpp"
25 #include "Thyra_LinearOpTester.hpp"
26 #include "Thyra_UnitTestHelpers.hpp"
27 #include "Thyra_DefaultSpmdVectorSpace.hpp"
28 
29 #include "EpetraLinearOpTestHelpers.hpp"
30 
31 #include "Thyra_DefaultBlockedLinearOp.hpp"
32 
33 #include "Teuchos_UnitTestHarness.hpp"
34 #include "Teuchos_DefaultComm.hpp"
35 
36 namespace {
37 
38 
39 
40 } // namespace
41 
42 
43 namespace Thyra {
44 
45 
46 //
47 // Unit Tests
48 //
49 
50 
51 TEUCHOS_UNIT_TEST( EpetraLinearOp, ScaledLinearOpBase )
52 {
53  using Teuchos::null;
54  using Teuchos::inOutArg;
55  using Teuchos::updateSuccess;
56  using Teuchos::rcp_dynamic_cast;
57  typedef ScalarTraits<double> ST;
58 
59  // Set up an EpetraLinearOp
60 
61  const RCP<const Epetra_Comm> comm = getEpetraComm();
62  const int numLocalRows = g_localDim;
63  const int numRows = numLocalRows * comm->NumProc();
64  const int numCols = numLocalRows / 2;
65  const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
66  const RCP<LinearOpBase<double> > epetraOp = nonconstEpetraLinearOp(epetraCrsM);
67 
68  const RCP<ScaledLinearOpBase<double> > scaledOp =
69  rcp_dynamic_cast<ScaledLinearOpBase<double> >(epetraOp, true);
70 
71  // Get the original mat-vec
72 
73  const double two = 2.0;
74 
75  const RCP<VectorBase<double> > rhs_vec =
76  createMember<double>(epetraOp->domain());
77  assign<double>(rhs_vec.ptr(), two);
78 
79  const RCP<VectorBase<double> > lhs_orig_vec =
80  createMember<double>(epetraOp->range());
81 
82  apply<double>(*epetraOp, Thyra::NOTRANS, *rhs_vec, lhs_orig_vec.ptr());
83 
84  if (g_dumpAll) {
85  out << "epetraOp = " << *epetraOp;
86  out << "rhs_vec = " << *rhs_vec;
87  out << "lhs_orig_vec = " << *lhs_orig_vec;
88  }
89 
90  // Scale the op from the left (row scaling)
91 
92  const double three = 3.0;
93 
94  const RCP<VectorBase<double> > row_scaling =
95  createMember<double>(epetraOp->range());
96  assign<double>(row_scaling.ptr(), three);
97 
98  scaledOp->scaleLeft(*row_scaling);
99 
100  if (g_dumpAll) {
101  out << "row_scaling = " << *row_scaling;
102  out << "epetraOp left scaled = " << *epetraOp;
103  }
104 
105  // Test that resulting left scaling
106 
107  const RCP<VectorBase<double> > lhs_left_scaled_vec =
108  createMember<double>(epetraOp->range());
109 
110  apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_left_scaled_vec.ptr());
111 
112  if (g_dumpAll) {
113  out << "lhs_left_scaled_vec = " << *lhs_left_scaled_vec;
114  }
115 
116  TEST_FLOATING_EQUALITY(
117  three * sum<double>(*lhs_orig_vec),
118  sum<double>(*lhs_left_scaled_vec),
119  as<double>(10.0 * ST::eps())
120  );
121 
122  // Left scale the matrix back
123 
124  const RCP<VectorBase<double> > inv_row_scaling =
125  createMember<double>(epetraOp->range());
126  reciprocal<double>(*row_scaling, inv_row_scaling.ptr());
127 
128  scaledOp->scaleLeft(*inv_row_scaling);
129 
130  if (g_dumpAll) {
131  out << "inv_row_scaling = " << *row_scaling;
132  out << "epetraOp left scaled back to orig = " << *epetraOp;
133  }
134 
135  const RCP<VectorBase<double> > lhs_orig2_vec =
136  createMember<double>(epetraOp->range());
137 
138  apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_orig2_vec.ptr());
139 
140  if (g_dumpAll) {
141  out << "lhs_orig2_vec = " << *lhs_orig2_vec;
142  }
143 
144  TEST_FLOATING_EQUALITY(
145  sum<double>(*lhs_orig_vec),
146  sum<double>(*lhs_orig2_vec),
147  as<double>(10.0 * ST::eps())
148  );
149  // NOTE: Above, it would ask for exact binary match except if one uses
150  // threading it will not match exactly!
151 
152  // Scale the op from the right (col scaling)
153 
154  const double four = 4.0;
155 
156  const RCP<VectorBase<double> > col_scaling =
157  createMember<double>(epetraOp->domain());
158  assign<double>(col_scaling.ptr(), four);
159 
160  scaledOp->scaleRight(*col_scaling);
161 
162  if (g_dumpAll) {
163  out << "col_scaling = " << *col_scaling;
164  out << "epetraOp right scaled = " << *epetraOp;
165  }
166 
167  // Test that resulting right scaling
168 
169  const RCP<VectorBase<double> > lhs_right_scaled_vec =
170  createMember<double>(epetraOp->range());
171 
172  apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_right_scaled_vec.ptr());
173 
174  if (g_dumpAll) {
175  out << "lhs_right_scaled_vec = " << *lhs_right_scaled_vec;
176  }
177 
178  TEST_FLOATING_EQUALITY(
179  four * sum<double>(*lhs_orig_vec),
180  sum<double>(*lhs_right_scaled_vec),
181  as<double>(10.0 * ST::eps())
182  );
183 
184  // Right scale the matrix back
185 
186  const RCP<VectorBase<double> > inv_col_scaling =
187  createMember<double>(epetraOp->domain());
188  reciprocal<double>(*col_scaling, inv_col_scaling.ptr());
189 
190  scaledOp->scaleRight(*inv_col_scaling);
191 
192  if (g_dumpAll) {
193  out << "inv_col_scaling = " << *col_scaling;
194  out << "epetraOp right scaled back to orig = " << *epetraOp;
195  }
196 
197  const RCP<VectorBase<double> > lhs_orig3_vec =
198  createMember<double>(epetraOp->range());
199 
200  apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_orig3_vec.ptr());
201 
202  if (g_dumpAll) {
203  out << "lhs_orig3_vec = " << *lhs_orig3_vec;
204  }
205 
206  TEST_FLOATING_EQUALITY(
207  sum<double>(*lhs_orig_vec),
208  sum<double>(*lhs_orig3_vec),
209  as<double>(10.0 * ST::eps())
210  );
211  // NOTE: Above, it would ask for exact binary match except if one uses
212  // threading it will not match exactly!
213 
214 }
215 
216 
217 TEUCHOS_UNIT_TEST( EpetraLinearOp, RowStatLinearOpBase )
218 {
219  using Teuchos::null;
220  using Teuchos::inOutArg;
221  using Teuchos::updateSuccess;
222  using Teuchos::rcp_dynamic_cast;
223  typedef ScalarTraits<double> ST;
224 
225  // Set up the EpetraLinearOp
226 
227  const RCP<const Epetra_Comm> comm = getEpetraComm();
228  const int numLocalRows = g_localDim;
229  const int numRows = numLocalRows * comm->NumProc();
230  const int numCols = numLocalRows / 2;
231  const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
232  const double two = 2.0;
233  epetraCrsM->PutScalar(-two); // put in negative two just to be extra "tricky"
234  const RCP<LinearOpBase<double> > epetraOp = nonconstEpetraLinearOp(epetraCrsM);
235 
236  const RCP<RowStatLinearOpBase<double> > rowStatOp =
237  rcp_dynamic_cast<RowStatLinearOpBase<double> >(epetraOp, true);
238 
239  if (g_dumpAll) {
240  out << "epetraOp = " << *epetraOp;
241  }
242 
243  // Get the inverse row sums
244 
245  const RCP<VectorBase<double> > inv_row_sums =
246  createMember<double>(epetraOp->range());
247  const RCP<VectorBase<double> > row_sums =
248  createMember<double>(epetraOp->range());
249 
250  rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
251  inv_row_sums.ptr());
252  rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
253  row_sums.ptr());
254 
255  if (g_dumpAll) {
256  out << "inv_row_sums = " << *inv_row_sums;
257  out << "row_sums = " << *row_sums;
258  }
259 
260  TEST_FLOATING_EQUALITY(
261  sum<double>(*inv_row_sums),
262  as<double>((1.0 / (two * numCols)) * numRows),
263  as<double>(10.0 * ST::eps())
264  );
265 
266  TEST_FLOATING_EQUALITY(
267  sum<double>(*row_sums),
268  as<double>(two * numCols * numRows),
269  as<double>(10.0 * ST::eps())
270  );
271 }
272 
273 RCP<Epetra_CrsMatrix> getMyEpetraMatrix(int numRows, int numCols, double shift=0.0)
274 {
275  const RCP<const Epetra_Comm> comm = getEpetraComm();
276 
277  const Epetra_Map rowMap(numRows, 0, *comm);
278  const Epetra_Map domainMap(numCols, numCols, 0, *comm);
279 
280  const RCP<Epetra_CrsMatrix> epetraCrsM =
281  rcp(new Epetra_CrsMatrix(Copy, rowMap,domainMap,0));
282 
283  Array<double> rowEntries(numCols);
284  Array<int> columnIndices(numCols);
285  for (int j = 0; j < numCols; ++j)
286  columnIndices[j] = j;
287 
288  const int numLocalRows = rowMap.NumMyElements();
289  for (int i = 0; i < numLocalRows; ++i) {
290  for (int j = 0; j < numCols; ++j) {
291  rowEntries[j] = as<double>(i+1) + as<double>(j+1) / 10 + shift;
292  }
293 
294  epetraCrsM->InsertMyValues( i, numCols, &rowEntries[0], &columnIndices[0] );
295  }
296 
297  epetraCrsM->FillComplete(domainMap,rowMap);
298  return epetraCrsM;
299 }
300 
301 TEUCHOS_UNIT_TEST( EpetraLinearOp, Blocked_ScaledLinearOpBase)
302 {
303  using Teuchos::null;
304  using Teuchos::inOutArg;
305  using Teuchos::updateSuccess;
306  using Teuchos::rcp_dynamic_cast;
307  // typedef ScalarTraits<double> ST; // unused
308 
309  // Set up the EpetraLinearOp
310 
311  const RCP<const Epetra_Comm> comm = getEpetraComm();
312  const int numLocalRows = g_localDim;
313  const int numRows = numLocalRows * comm->NumProc();
314  const int numCols = numLocalRows ;
315 
316  out << "numRows = " << numRows << ", numCols = " << numCols << std::endl;
317 
318  const RCP<Epetra_CrsMatrix> epetraCrsM00_base = getMyEpetraMatrix(numRows, numRows);
319  const RCP<Epetra_CrsMatrix> epetraCrsM01_base = getMyEpetraMatrix(numRows, numCols);
320  const RCP<Epetra_CrsMatrix> epetraCrsM10_base = getMyEpetraMatrix(numCols, numRows);
321  const RCP<Epetra_CrsMatrix> epetraCrsM11_base = getMyEpetraMatrix(numCols, numCols);
322  epetraCrsM00_base->PutScalar(2.0);
323  epetraCrsM01_base->PutScalar(-8.0);
324  epetraCrsM10_base->PutScalar(-9.0);
325  epetraCrsM11_base->PutScalar(3.0);
326  const RCP<Epetra_CrsMatrix> epetraCrsM00 = getMyEpetraMatrix(numRows, numRows);
327  const RCP<Epetra_CrsMatrix> epetraCrsM01 = getMyEpetraMatrix(numRows, numCols);
328  const RCP<Epetra_CrsMatrix> epetraCrsM10 = getMyEpetraMatrix(numCols, numRows);
329  const RCP<Epetra_CrsMatrix> epetraCrsM11 = getMyEpetraMatrix(numCols, numCols);
330  epetraCrsM00->PutScalar(2.0);
331  epetraCrsM01->PutScalar(-8.0);
332  epetraCrsM10->PutScalar(-9.0);
333  epetraCrsM11->PutScalar(3.0);
334 
335  const RCP<const LinearOpBase<double> > op00_base = epetraLinearOp(epetraCrsM00_base);
336  const RCP<const LinearOpBase<double> > op01_base = epetraLinearOp(epetraCrsM01_base);
337  const RCP<const LinearOpBase<double> > op10_base = epetraLinearOp(epetraCrsM10_base);
338  const RCP<const LinearOpBase<double> > op11_base = epetraLinearOp(epetraCrsM11_base);
339 
340  const RCP<LinearOpBase<double> > op00 = nonconstEpetraLinearOp(epetraCrsM00);
341  const RCP<LinearOpBase<double> > op01 = nonconstEpetraLinearOp(epetraCrsM01);
342  const RCP<LinearOpBase<double> > op10 = nonconstEpetraLinearOp(epetraCrsM10);
343  const RCP<LinearOpBase<double> > op11 = nonconstEpetraLinearOp(epetraCrsM11);
344 
345  const RCP<const LinearOpBase<double> > blocked_base = block2x2(op00_base,op01_base,op10_base,op11_base);
346  const RCP<LinearOpBase<double> > blocked = nonconstBlock2x2(op00,op01,op10,op11);
347 
348  const RCP<VectorBase<double> > left_scale = createMember<double>(blocked_base->range());
349  const RCP<VectorBase<double> > right_scale = createMember<double>(blocked_base->domain());
350 
351  put_scalar(7.0,left_scale.ptr());
352  put_scalar(-4.0,right_scale.ptr());
353 
354  rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleLeft(*left_scale);
355 
356  {
357  LinearOpTester<double> tester;
358  tester.set_all_error_tol(1e-10);
359  tester.show_all_tests(true);
360  tester.dump_all(true);
361  tester.num_random_vectors(5);
362  const RCP<const LinearOpBase<double> > left_op = Thyra::diagonal(left_scale);
363  const RCP<const LinearOpBase<double> > ref_op = multiply(left_op,blocked_base);
364 
365  updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
366  }
367 
368  rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleRight(*right_scale);
369 
370  {
371  LinearOpTester<double> tester;
372  tester.set_all_error_tol(1e-10);
373  tester.show_all_tests(true);
374  tester.dump_all(true);
375  tester.num_random_vectors(5);
376  const RCP<const LinearOpBase<double> > left_op = Thyra::diagonal(left_scale);
377  const RCP<const LinearOpBase<double> > right_op = Thyra::diagonal(right_scale);
378  const RCP<const LinearOpBase<double> > ref_op = multiply(left_op,blocked_base,right_op);
379 
380  updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
381  }
382 }
383 
384 TEUCHOS_UNIT_TEST( EpetraLinearOp, Blocked_RowStatLinearOpBase )
385 {
386  using Teuchos::null;
387  using Teuchos::inOutArg;
388  using Teuchos::updateSuccess;
389  using Teuchos::rcp_dynamic_cast;
390  typedef ScalarTraits<double> ST;
391 
392  // Set up the EpetraLinearOp
393 
394  const RCP<const Epetra_Comm> comm = getEpetraComm();
395  const int numLocalRows = g_localDim;
396  const int numRows = numLocalRows * comm->NumProc();
397  const int numCols = numLocalRows / 2;
398 
399  out << "numRows = " << numRows << ", numCols = " << numCols << std::endl;
400 
401  const RCP<Epetra_CrsMatrix> epetraCrsM00 = getEpetraMatrix(numRows, numRows);
402  const RCP<Epetra_CrsMatrix> epetraCrsM01 = getEpetraMatrix(numRows, numCols);
403  const RCP<Epetra_CrsMatrix> epetraCrsM10 = getEpetraMatrix(numCols, numRows);
404  const RCP<Epetra_CrsMatrix> epetraCrsM11 = getEpetraMatrix(numCols, numCols);
405  epetraCrsM00->PutScalar(2.0);
406  epetraCrsM01->PutScalar(-8.0);
407  epetraCrsM10->PutScalar(-9.0);
408  epetraCrsM11->PutScalar(3.0);
409 
410  const RCP<const LinearOpBase<double> > op00 = epetraLinearOp(epetraCrsM00);
411  const RCP<const LinearOpBase<double> > op01 = epetraLinearOp(epetraCrsM01);
412  const RCP<const LinearOpBase<double> > op10 = epetraLinearOp(epetraCrsM10);
413  const RCP<const LinearOpBase<double> > op11 = epetraLinearOp(epetraCrsM11);
414 
415  const RCP<const LinearOpBase<double> > blocked = block2x2(op00,op01,op10,op11);
416 
417  const RCP<const RowStatLinearOpBase<double> > rowStatOp =
418  rcp_dynamic_cast<const RowStatLinearOpBase<double> >(blocked, true);
419 
420  if (g_dumpAll) {
421  out << "epetraOp = " << *blocked;
422  }
423 
424  // Get the inverse row sums
425 
426  const RCP<VectorBase<double> > inv_row_sums =
427  createMember<double>(blocked->range());
428  const RCP<VectorBase<double> > row_sums =
429  createMember<double>(blocked->range());
430 
431  rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
432  inv_row_sums.ptr());
433  rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
434  row_sums.ptr());
435 
436  if (g_dumpAll) {
437  out << "inv_row_sums = " << *inv_row_sums;
438  out << "row_sums = " << *row_sums;
439  }
440 
441  TEST_FLOATING_EQUALITY(
442  sum<double>(*inv_row_sums),
443  as<double>((1.0/(numRows*2.0+numCols*8.0))*numRows + (1.0/(numRows*9.0+numCols*3.0))*numCols),
444  as<double>(10.0 * ST::eps())
445  );
446  TEST_FLOATING_EQUALITY(
447  sum<double>(*row_sums),
448  as<double>((numRows*2.0+numCols*8.0)*numRows + (numRows*9.0+numCols*3.0)*numCols),
449  as<double>(10.0 * ST::eps())
450  );
451 }
452 
453 TEUCHOS_UNIT_TEST( EpetraLinearOp, Blocked_ScalingWithMultiVectors)
454 {
455  using Teuchos::null;
456  using Teuchos::inOutArg;
457  using Teuchos::updateSuccess;
458  using Teuchos::rcp_dynamic_cast;
459  typedef ScalarTraits<double> ST;
460 
461  // Set up the EpetraLinearOp
462 
463  const RCP<const Epetra_Comm> comm = getEpetraComm();
464  const RCP<const Teuchos::Comm<Ordinal> > tComm =
465  Teuchos::DefaultComm<Ordinal>::getComm();
466  const int numLocalRows = 4;
467  const int numRows = numLocalRows * comm->NumProc();
468  const int numCols = numLocalRows / 2;
469 
470  out << "numRows = " << numRows << ", numCols = " << numCols << std::endl;
471 
472  const RCP<Epetra_CrsMatrix> epetraCrsM00 = getMyEpetraMatrix(numRows, numRows);
473  const RCP<Epetra_CrsMatrix> epetraCrsM00_base = getMyEpetraMatrix(numRows, numRows);
474  epetraCrsM00->PutScalar(2.0);
475  epetraCrsM00_base->PutScalar(2.0);
476 
477  const RCP<LinearOpBase<double> > op00 = nonconstEpetraLinearOp(epetraCrsM00);
478  const RCP<const LinearOpBase<double> > op00_base = epetraLinearOp(epetraCrsM00_base);
479 
480  RCP<const Thyra::VectorSpaceBase<double> > vs_0 = op00->range();
481  RCP<const Thyra::VectorSpaceBase<double> > vs_1 = Thyra::locallyReplicatedDefaultSpmdVectorSpace<double>(tComm,numCols);
482 
483  RCP<Thyra::MultiVectorBase<double> > vec_01 = Thyra::createMembers(vs_0,numCols);
484  RCP<Thyra::MultiVectorBase<double> > vec_10t = Thyra::createMembers(op00->domain(),numCols); // tranposed
485  RCP<Thyra::MultiVectorBase<double> > vec_01_base = Thyra::createMembers(vs_0,numCols);
486  RCP<Thyra::MultiVectorBase<double> > vec_10t_base = Thyra::createMembers(op00->domain(),numCols); // tranposed
487  const RCP<LinearOpBase<double> > op10t = vec_10t;
488  const RCP<const LinearOpBase<double> > op10t_base = vec_10t_base;
489  assign(vec_01.ptr(),-8.0);
490  assign(vec_10t.ptr(),-9.0);
491  assign(vec_01_base.ptr(),-8.0);
492  assign(vec_10t_base.ptr(),-9.0);
493 
494  const RCP<LinearOpBase<double> > op01 = vec_01;
495  const RCP<LinearOpBase<double> > op10 = nonconstAdjoint(op10t);
496  const RCP<LinearOpBase<double> > op11 = nonconstZero(vec_01->domain(),vec_01->domain());
497 
498  const RCP<const LinearOpBase<double> > op01_base = vec_01_base;
499  const RCP<const LinearOpBase<double> > op10_base = adjoint(op10t_base);
500  const RCP<const LinearOpBase<double> > op11_base = zero(vec_01->domain(),vec_01->domain());
501 
502  out << "FIRST" << std::endl;
503  const RCP<LinearOpBase<double> > blocked = nonconstBlock2x2(op00,op01,op10,op11);
504  out << "SECOND" << Teuchos::describe(*blocked,Teuchos::VERB_EXTREME) << std::endl;
505  const RCP<const LinearOpBase<double> > blocked_base = block2x2(op00_base,op01_base,op10_base,op11_base);
506 
507  const RCP<const RowStatLinearOpBase<double> > rowStatOp =
508  rcp_dynamic_cast<const RowStatLinearOpBase<double> >(blocked, true);
509 
510  // Get the inverse row sums
511 
512  const RCP<VectorBase<double> > inv_row_sums =
513  createMember<double>(blocked->range());
514  const RCP<VectorBase<double> > row_sums =
515  createMember<double>(blocked->range());
516 
517  rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
518  inv_row_sums.ptr());
519  rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
520  row_sums.ptr());
521 
522  TEST_FLOATING_EQUALITY(
523  sum<double>(*inv_row_sums),
524  as<double>((1.0/(numRows*2.0+numCols*8.0))*numRows + (1.0/(numRows*9.0+numCols*0.0))*numCols),
525  as<double>(10.0 * ST::eps())
526  );
527  TEST_FLOATING_EQUALITY(
528  sum<double>(*row_sums),
529  as<double>((numRows*2.0+numCols*8.0)*numRows + (numRows*9.0+numCols*0.0)*numCols),
530  as<double>(10.0 * ST::eps())
531  );
532 
533  {
534  const RCP<VectorBase<double> > left_scale = createMember<double>(blocked_base->range());
535  const RCP<VectorBase<double> > right_scale = createMember<double>(blocked_base->domain());
536 
537  put_scalar(7.0,left_scale.ptr());
538  put_scalar(-4.0,right_scale.ptr());
539 
540  rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleLeft(*left_scale);
541 
542  {
543  LinearOpTester<double> tester;
544  tester.set_all_error_tol(1e-10);
545  tester.show_all_tests(true);
546  tester.dump_all(false);
547  tester.num_random_vectors(2);
548  const RCP<const LinearOpBase<double> > left_op = Thyra::diagonal(left_scale);
549  const RCP<const LinearOpBase<double> > ref_op = multiply(left_op,blocked_base);
550 
551  updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
552  }
553 
554  rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleRight(*right_scale);
555 
556  {
557  LinearOpTester<double> tester;
558  tester.set_all_error_tol(1e-10);
559  tester.show_all_tests(true);
560  tester.dump_all(false);
561  tester.num_random_vectors(5);
562  const RCP<const LinearOpBase<double> > left_op = Thyra::diagonal(left_scale);
563  const RCP<const LinearOpBase<double> > right_op = Thyra::diagonal(right_scale);
564  const RCP<const LinearOpBase<double> > ref_op = multiply(left_op,blocked_base,right_op);
565 
566  updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
567  }
568  }
569 }
570 
571 
572 TEUCHOS_UNIT_TEST( EpetraLinearOp, rectangular )
573 {
574  using Teuchos::null;
575  using Teuchos::inOutArg;
576  using Teuchos::updateSuccess;
577 
578  const RCP<const Epetra_Comm> comm = getEpetraComm();
579  const int numProcs = comm->NumProc();
580 
581  const int numLocalRows = g_localDim;
582  const int numRows = numLocalRows * comm->NumProc();
583  const int numCols = numLocalRows / 2;
584 
585  const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
586 
587  const RCP<const LinearOpBase<double> > epetraOp = epetraLinearOp(epetraCrsM);
588 
589  LinearOpTester<double> linearOpTester;
590  linearOpTester.check_adjoint(numProcs == 1);
591  linearOpTester.show_all_tests(g_show_all_tests);
592  linearOpTester.dump_all(g_dumpAll);
593  updateSuccess(linearOpTester.check(*epetraOp, inOutArg(out)), success);
594 
595  // NOTE: Above, it would seem the Epetra_CrsMatrix::Apply(...) does not work
596  // when doing and adjoint where the RowMap has empty processes.
597 
598 }
599 
600 
601 TEUCHOS_UNIT_TEST( EpetraLinearOp, blocked_op )
602 {
603 
604  if (Teuchos::GlobalMPISession::getNProc() > 2) {
605  out << "Skipping test if numProc > 2 since it fails for some reason\n";
606  return;
607  }
608 
609  using Teuchos::describe;
610 
611  // build sub operators
612  RCP<const LinearOpBase<double> > A00 =
613  epetraLinearOp(getEpetraMatrix(4,4,0));
614  RCP<const LinearOpBase<double> > A01 =
615  epetraLinearOp(getEpetraMatrix(4,3,1));
616  RCP<const LinearOpBase<double> > A02 =
617  epetraLinearOp(getEpetraMatrix(4,2,2));
618  RCP<const LinearOpBase<double> > A10 =
619  epetraLinearOp(getEpetraMatrix(3,4,3));
620  RCP<const LinearOpBase<double> > A11 =
621  epetraLinearOp(getEpetraMatrix(3,3,4));
622  RCP<const LinearOpBase<double> > A12 =
623  epetraLinearOp(getEpetraMatrix(3,2,5));
624  RCP<const LinearOpBase<double> > A20 =
625  epetraLinearOp(getEpetraMatrix(2,4,6));
626  RCP<const LinearOpBase<double> > A21 =
627  epetraLinearOp(getEpetraMatrix(2,3,8));
628  RCP<const LinearOpBase<double> > A22 =
629  epetraLinearOp(getEpetraMatrix(2,2,8));
630 
631  const Teuchos::EVerbosityLevel verbLevel =
632  (g_dumpAll ? Teuchos::VERB_HIGH : Teuchos::VERB_MEDIUM);
633 
634  out << "Sub operators built" << std::endl;
635 
636  {
637  // build composite operator
638  RCP<const LinearOpBase<double> > A =
639  block2x2<double>(
640  block2x2<double>(A00, A01, A10, A11), block2x1<double>(A02,A12),
641  block1x2<double>(A20, A21), A22
642  );
643 
644  out << "First composite operator built" << std::endl;
645 
646  // build vectors for use in apply
647  RCP<MultiVectorBase<double> > x = createMembers<double>(A->domain(), 3);
648  RCP<MultiVectorBase<double> > y = createMembers<double>(A->range(), 3);
649 
650  randomize(-1.0, 1.0, x.ptr());
651 
652  out << "A = \n" << describe(*A, verbLevel) << std::endl;
653  out << "x = \n" << describe(*x, verbLevel) << std::endl;
654  out << "y = \n" << describe(*y, verbLevel) << std::endl;
655 
656  // perform a matrix vector multiply
657  apply(*A, NOTRANS, *x, y.ptr());
658 
659  out << "First composite operator completed" << std::endl;
660  }
661 
662  {
663  RCP<const LinearOpBase<double> > A = block2x2<double>(
664  A11, block1x2<double>(A10, A12),
665  block2x1<double>(A01, A21), block2x2<double>(A00, A02, A20, A22)
666  );
667 
668  out << "Second composite operator built" << std::endl;
669 
670  // build vectors for use in apply
671  RCP<MultiVectorBase<double> > x = createMembers<double>(A->domain(), 3);
672  RCP<MultiVectorBase<double> > y = createMembers<double>(A->range(), 3);
673 
674  randomize(-1.0, 1.0, x.ptr());
675 
676  out << "A = \n" << describe(*A, verbLevel) << std::endl;
677  out << "x = \n" << describe(*x, verbLevel) << std::endl;
678  out << "y = \n" << describe(*y, verbLevel) << std::endl;
679 
680  // perform a matrix vector multiply
681  apply(*A, NOTRANS, *x, y.ptr());
682 
683  out << "Second composite operator completed" << std::endl;
684  }
685 
686  out << "Test complete" << std::endl;
687 
688 }
689 
690 
691 } // namespace Thyra
Thyra::EpetraLinearOp
Concrete LinearOpBase adapter subclass for Epetra_Operator object.
Definition: Thyra_EpetraLinearOp.hpp:46
Thyra_EpetraLinearOp.hpp
Thyra::getMyEpetraMatrix
RCP< Epetra_CrsMatrix > getMyEpetraMatrix(int numRows, int numCols, double shift=0.0)
Definition: EpetraLinearOp_UnitTests.cpp:273
Teuchos::DefaultComm::getComm
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Thyra::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)
Definition: EpetraOperatorWrapper_UnitTests.cpp:49
Teuchos::updateSuccess
void updateSuccess(const bool result, bool &success)
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos_DefaultComm.hpp
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos_UnitTestHarness.hpp
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Generated on Thu Nov 21 2024 09:21:28 for Thyra Package Browser (Single Doxygen Collection) by   doxygen 1.8.5