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" 65 #include "Thyra_DefaultBlockedLinearOp.hpp" 67 #include "Teuchos_UnitTestHarness.hpp" 68 #include "Teuchos_DefaultComm.hpp" 88 using Teuchos::inOutArg;
89 using Teuchos::updateSuccess;
90 using Teuchos::rcp_dynamic_cast;
91 typedef ScalarTraits<double> ST;
95 const RCP<const Epetra_Comm> comm = getEpetraComm();
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);
102 const RCP<ScaledLinearOpBase<double> > scaledOp =
103 rcp_dynamic_cast<ScaledLinearOpBase<double> >(epetraOp,
true);
107 const double two = 2.0;
109 const RCP<VectorBase<double> > rhs_vec =
110 createMember<double>(epetraOp->domain());
111 assign<double>(rhs_vec.ptr(), two);
113 const RCP<VectorBase<double> > lhs_orig_vec =
114 createMember<double>(epetraOp->range());
116 apply<double>(*epetraOp, Thyra::NOTRANS, *rhs_vec, lhs_orig_vec.ptr());
119 out <<
"epetraOp = " << *epetraOp;
120 out <<
"rhs_vec = " << *rhs_vec;
121 out <<
"lhs_orig_vec = " << *lhs_orig_vec;
126 const double three = 3.0;
128 const RCP<VectorBase<double> > row_scaling =
129 createMember<double>(epetraOp->range());
130 assign<double>(row_scaling.ptr(), three);
132 scaledOp->scaleLeft(*row_scaling);
135 out <<
"row_scaling = " << *row_scaling;
136 out <<
"epetraOp left scaled = " << *epetraOp;
141 const RCP<VectorBase<double> > lhs_left_scaled_vec =
142 createMember<double>(epetraOp->range());
144 apply<double>(*epetraOp,
NOTRANS, *rhs_vec, lhs_left_scaled_vec.ptr());
147 out <<
"lhs_left_scaled_vec = " << *lhs_left_scaled_vec;
150 TEST_FLOATING_EQUALITY(
151 three * sum<double>(*lhs_orig_vec),
152 sum<double>(*lhs_left_scaled_vec),
153 as<double>(10.0 * ST::eps())
158 const RCP<VectorBase<double> > inv_row_scaling =
159 createMember<double>(epetraOp->range());
160 reciprocal<double>(*row_scaling, inv_row_scaling.ptr());
162 scaledOp->scaleLeft(*inv_row_scaling);
165 out <<
"inv_row_scaling = " << *row_scaling;
166 out <<
"epetraOp left scaled back to orig = " << *epetraOp;
169 const RCP<VectorBase<double> > lhs_orig2_vec =
170 createMember<double>(epetraOp->range());
172 apply<double>(*epetraOp,
NOTRANS, *rhs_vec, lhs_orig2_vec.ptr());
175 out <<
"lhs_orig2_vec = " << *lhs_orig2_vec;
178 TEST_FLOATING_EQUALITY(
179 sum<double>(*lhs_orig_vec),
180 sum<double>(*lhs_orig2_vec),
181 as<double>(10.0 * ST::eps())
188 const double four = 4.0;
190 const RCP<VectorBase<double> > col_scaling =
191 createMember<double>(epetraOp->domain());
192 assign<double>(col_scaling.ptr(), four);
194 scaledOp->scaleRight(*col_scaling);
197 out <<
"col_scaling = " << *col_scaling;
198 out <<
"epetraOp right scaled = " << *epetraOp;
203 const RCP<VectorBase<double> > lhs_right_scaled_vec =
204 createMember<double>(epetraOp->range());
206 apply<double>(*epetraOp,
NOTRANS, *rhs_vec, lhs_right_scaled_vec.ptr());
209 out <<
"lhs_right_scaled_vec = " << *lhs_right_scaled_vec;
212 TEST_FLOATING_EQUALITY(
213 four * sum<double>(*lhs_orig_vec),
214 sum<double>(*lhs_right_scaled_vec),
215 as<double>(10.0 * ST::eps())
220 const RCP<VectorBase<double> > inv_col_scaling =
221 createMember<double>(epetraOp->domain());
222 reciprocal<double>(*col_scaling, inv_col_scaling.ptr());
224 scaledOp->scaleRight(*inv_col_scaling);
227 out <<
"inv_col_scaling = " << *col_scaling;
228 out <<
"epetraOp right scaled back to orig = " << *epetraOp;
231 const RCP<VectorBase<double> > lhs_orig3_vec =
232 createMember<double>(epetraOp->range());
234 apply<double>(*epetraOp,
NOTRANS, *rhs_vec, lhs_orig3_vec.ptr());
237 out <<
"lhs_orig3_vec = " << *lhs_orig3_vec;
240 TEST_FLOATING_EQUALITY(
241 sum<double>(*lhs_orig_vec),
242 sum<double>(*lhs_orig3_vec),
243 as<double>(10.0 * ST::eps())
254 using Teuchos::inOutArg;
255 using Teuchos::updateSuccess;
256 using Teuchos::rcp_dynamic_cast;
257 typedef ScalarTraits<double> ST;
261 const RCP<const Epetra_Comm> comm = getEpetraComm();
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);
268 const RCP<LinearOpBase<double> > epetraOp = nonconstEpetraLinearOp(epetraCrsM);
270 const RCP<RowStatLinearOpBase<double> > rowStatOp =
271 rcp_dynamic_cast<RowStatLinearOpBase<double> >(epetraOp,
true);
274 out <<
"epetraOp = " << *epetraOp;
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());
284 rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
286 rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
290 out <<
"inv_row_sums = " << *inv_row_sums;
291 out <<
"row_sums = " << *row_sums;
294 TEST_FLOATING_EQUALITY(
295 sum<double>(*inv_row_sums),
296 as<double>((1.0 / (two * numCols)) * numRows),
297 as<double>(10.0 * ST::eps())
300 TEST_FLOATING_EQUALITY(
301 sum<double>(*row_sums),
302 as<double>(two * numCols * numRows),
303 as<double>(10.0 * ST::eps())
309 const RCP<const Epetra_Comm> comm = getEpetraComm();
311 const Epetra_Map rowMap(numRows, 0, *comm);
312 const Epetra_Map domainMap(numCols, numCols, 0, *comm);
314 const RCP<Epetra_CrsMatrix> epetraCrsM =
315 rcp(
new Epetra_CrsMatrix(Copy, rowMap,domainMap,0));
317 Array<double> rowEntries(numCols);
318 Array<int> columnIndices(numCols);
319 for (
int j = 0; j < numCols; ++j)
320 columnIndices[j] = j;
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;
328 epetraCrsM->InsertMyValues( i, numCols, &rowEntries[0], &columnIndices[0] );
331 epetraCrsM->FillComplete(domainMap,rowMap);
338 using Teuchos::inOutArg;
339 using Teuchos::updateSuccess;
340 using Teuchos::rcp_dynamic_cast;
345 const RCP<const Epetra_Comm> comm = getEpetraComm();
347 const int numRows = numLocalRows * comm->NumProc();
348 const int numCols = numLocalRows ;
350 out <<
"numRows = " << numRows <<
", numCols = " << numCols << std::endl;
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);
364 epetraCrsM00->PutScalar(2.0);
365 epetraCrsM01->PutScalar(-8.0);
366 epetraCrsM10->PutScalar(-9.0);
367 epetraCrsM11->PutScalar(3.0);
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);
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);
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);
382 const RCP<VectorBase<double> > left_scale = createMember<double>(blocked_base->range());
383 const RCP<VectorBase<double> > right_scale = createMember<double>(blocked_base->domain());
385 put_scalar(7.0,left_scale.ptr());
386 put_scalar(-4.0,right_scale.ptr());
388 rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleLeft(*left_scale);
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);
399 updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
402 rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleRight(*right_scale);
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);
414 updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
421 using Teuchos::inOutArg;
422 using Teuchos::updateSuccess;
423 using Teuchos::rcp_dynamic_cast;
424 typedef ScalarTraits<double> ST;
428 const RCP<const Epetra_Comm> comm = getEpetraComm();
430 const int numRows = numLocalRows * comm->NumProc();
431 const int numCols = numLocalRows / 2;
433 out <<
"numRows = " << numRows <<
", numCols = " << numCols << std::endl;
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);
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);
449 const RCP<const LinearOpBase<double> > blocked = block2x2(op00,op01,op10,op11);
451 const RCP<const RowStatLinearOpBase<double> > rowStatOp =
452 rcp_dynamic_cast<
const RowStatLinearOpBase<double> >(blocked,
true);
455 out <<
"epetraOp = " << *blocked;
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());
465 rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
467 rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
471 out <<
"inv_row_sums = " << *inv_row_sums;
472 out <<
"row_sums = " << *row_sums;
475 TEST_FLOATING_EQUALITY(
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())
480 TEST_FLOATING_EQUALITY(
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())
490 using Teuchos::inOutArg;
491 using Teuchos::updateSuccess;
492 using Teuchos::rcp_dynamic_cast;
493 typedef ScalarTraits<double> ST;
497 const RCP<const Epetra_Comm> comm = getEpetraComm();
498 const RCP<const Teuchos::Comm<Ordinal> > tComm =
499 Teuchos::DefaultComm<Ordinal>::getComm();
500 const int numLocalRows = 4;
501 const int numRows = numLocalRows * comm->NumProc();
502 const int numCols = numLocalRows / 2;
504 out <<
"numRows = " << numRows <<
", numCols = " << numCols << std::endl;
507 const RCP<Epetra_CrsMatrix> epetraCrsM00_base =
getMyEpetraMatrix(numRows, numRows);
508 epetraCrsM00->PutScalar(2.0);
509 epetraCrsM00_base->PutScalar(2.0);
511 const RCP<LinearOpBase<double> > op00 = nonconstEpetraLinearOp(epetraCrsM00);
512 const RCP<const LinearOpBase<double> > op00_base = epetraLinearOp(epetraCrsM00_base);
514 RCP<const Thyra::VectorSpaceBase<double> > vs_0 = op00->range();
515 RCP<const Thyra::VectorSpaceBase<double> > vs_1 = Thyra::locallyReplicatedDefaultSpmdVectorSpace<double>(tComm,numCols);
517 RCP<Thyra::MultiVectorBase<double> > vec_01 = Thyra::createMembers(vs_0,numCols);
518 RCP<Thyra::MultiVectorBase<double> > vec_10t = Thyra::createMembers(op00->domain(),numCols);
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);
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);
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());
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());
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);
541 const RCP<const RowStatLinearOpBase<double> > rowStatOp =
542 rcp_dynamic_cast<
const RowStatLinearOpBase<double> >(blocked,
true);
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());
551 rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
553 rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM,
556 TEST_FLOATING_EQUALITY(
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())
561 TEST_FLOATING_EQUALITY(
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())
568 const RCP<VectorBase<double> > left_scale = createMember<double>(blocked_base->range());
569 const RCP<VectorBase<double> > right_scale = createMember<double>(blocked_base->domain());
571 put_scalar(7.0,left_scale.ptr());
572 put_scalar(-4.0,right_scale.ptr());
574 rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleLeft(*left_scale);
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);
585 updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
588 rcp_dynamic_cast<ScaledLinearOpBase<double> >(blocked)->scaleRight(*right_scale);
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);
600 updateSuccess(tester.compare(*ref_op, *blocked, ptrFromRef(out)), success);
609 using Teuchos::inOutArg;
610 using Teuchos::updateSuccess;
612 const RCP<const Epetra_Comm> comm = getEpetraComm();
613 const int numProcs = comm->NumProc();
616 const int numRows = numLocalRows * comm->NumProc();
617 const int numCols = numLocalRows / 2;
619 const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
621 const RCP<const LinearOpBase<double> > epetraOp = epetraLinearOp(epetraCrsM);
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);
638 if (Teuchos::GlobalMPISession::getNProc() > 2) {
639 out <<
"Skipping test if numProc > 2 since it fails for some reason\n";
643 using Teuchos::describe;
646 RCP<const LinearOpBase<double> > A00 =
647 epetraLinearOp(getEpetraMatrix(4,4,0));
648 RCP<const LinearOpBase<double> > A01 =
649 epetraLinearOp(getEpetraMatrix(4,3,1));
650 RCP<const LinearOpBase<double> > A02 =
651 epetraLinearOp(getEpetraMatrix(4,2,2));
652 RCP<const LinearOpBase<double> > A10 =
653 epetraLinearOp(getEpetraMatrix(3,4,3));
654 RCP<const LinearOpBase<double> > A11 =
655 epetraLinearOp(getEpetraMatrix(3,3,4));
656 RCP<const LinearOpBase<double> > A12 =
657 epetraLinearOp(getEpetraMatrix(3,2,5));
658 RCP<const LinearOpBase<double> > A20 =
659 epetraLinearOp(getEpetraMatrix(2,4,6));
660 RCP<const LinearOpBase<double> > A21 =
661 epetraLinearOp(getEpetraMatrix(2,3,8));
662 RCP<const LinearOpBase<double> > A22 =
663 epetraLinearOp(getEpetraMatrix(2,2,8));
665 const Teuchos::EVerbosityLevel verbLevel =
666 (g_dumpAll ? Teuchos::VERB_HIGH : Teuchos::VERB_MEDIUM);
668 out <<
"Sub operators built" << std::endl;
672 RCP<const LinearOpBase<double> > A =
674 block2x2<double>(A00, A01, A10, A11), block2x1<double>(A02,A12),
675 block1x2<double>(A20, A21), A22
678 out <<
"First composite operator built" << std::endl;
681 RCP<MultiVectorBase<double> > x = createMembers<double>(A->domain(), 3);
682 RCP<MultiVectorBase<double> > y = createMembers<double>(A->range(), 3);
684 randomize(-1.0, 1.0, x.ptr());
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;
691 apply(*A,
NOTRANS, *x, y.ptr());
693 out <<
"First composite operator completed" << std::endl;
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)
702 out <<
"Second composite operator built" << std::endl;
705 RCP<MultiVectorBase<double> > x = createMembers<double>(A->domain(), 3);
706 RCP<MultiVectorBase<double> > y = createMembers<double>(A->range(), 3);
708 randomize(-1.0, 1.0, x.ptr());
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;
715 apply(*A,
NOTRANS, *x, y.ptr());
717 out <<
"Second composite operator completed" << std::endl;
720 out <<
"Test complete" << std::endl;
Concrete LinearOpBase adapter subclass for Epetra_Operator object.
RCP< Epetra_CrsMatrix > getMyEpetraMatrix(int numRows, int numCols, double shift=0.0)
TEUCHOS_UNIT_TEST(EpetraOperatorWrapper, basic)