Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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:
141//
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
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
69
70namespace {
71
72
73
74} // namespace
75
76
77namespace Thyra {
78
79
80//
81// Unit Tests
82//
83
84
85TEUCHOS_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
251TEUCHOS_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
307RCP<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
335TEUCHOS_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
418TEUCHOS_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
487TEUCHOS_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();
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
TEUCHOS_UNIT_TEST(Comm, Issue1029)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Ptr< T > ptr() const
Concrete LinearOpBase adapter subclass for Epetra_Operator object.
NOTRANS
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Copy
void updateSuccess(const bool result, bool &success)
RCP< Epetra_CrsMatrix > getMyEpetraMatrix(int numRows, int numCols, double shift=0.0)