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