Teko Version of the Day
Loading...
Searching...
No Matches
Teko_Utilities.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_Config.h"
48#include "Teko_Utilities.hpp"
49
50// Thyra includes
51#include "Thyra_MultiVectorStdOps.hpp"
52#include "Thyra_ZeroLinearOpBase.hpp"
53#include "Thyra_DefaultDiagonalLinearOp.hpp"
54#include "Thyra_DefaultAddedLinearOp.hpp"
55#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
56#include "Thyra_DefaultMultipliedLinearOp.hpp"
57#include "Thyra_DefaultZeroLinearOp.hpp"
58#include "Thyra_DefaultProductMultiVector.hpp"
59#include "Thyra_DefaultProductVectorSpace.hpp"
60#include "Thyra_MultiVectorStdOps.hpp"
61#include "Thyra_VectorStdOps.hpp"
62#include "Thyra_SpmdVectorBase.hpp"
63#include <utility>
64
65#ifdef TEKO_HAVE_EPETRA
66#include "Thyra_EpetraExtDiagScaledMatProdTransformer.hpp"
67#include "Thyra_EpetraExtDiagScalingTransformer.hpp"
68#include "Thyra_EpetraExtAddTransformer.hpp"
69#include "Thyra_get_Epetra_Operator.hpp"
70#include "Thyra_EpetraThyraWrappers.hpp"
71#include "Thyra_EpetraOperatorWrapper.hpp"
72#include "Thyra_EpetraLinearOp.hpp"
73#endif
74
75// Teuchos includes
76#include "Teuchos_Array.hpp"
77
78// Epetra includes
79#ifdef TEKO_HAVE_EPETRA
80#include "Epetra_Operator.h"
81#include "Epetra_CrsGraph.h"
82#include "Epetra_CrsMatrix.h"
83#include "Epetra_Vector.h"
84#include "Epetra_Map.h"
85
86#include "EpetraExt_Transpose_RowMatrix.h"
87#include "EpetraExt_MatrixMatrix.h"
88
89#include "Teko_EpetraHelpers.hpp"
90#include "Teko_EpetraOperatorWrapper.hpp"
91#endif
92
93// Anasazi includes
94#include "AnasaziBasicEigenproblem.hpp"
95#include "AnasaziThyraAdapter.hpp"
96#include "AnasaziBlockKrylovSchurSolMgr.hpp"
97#include "AnasaziBlockKrylovSchur.hpp"
98#include "AnasaziStatusTestMaxIters.hpp"
99
100// Isorropia includes
101#ifdef Teko_ENABLE_Isorropia
102#include "Isorropia_EpetraProber.hpp"
103#endif
104
105// Teko includes
106#include "Teko_TpetraHelpers.hpp"
107#include "Teko_TpetraOperatorWrapper.hpp"
108
109// Tpetra
110#include "Thyra_TpetraLinearOp.hpp"
111#include "Tpetra_CrsMatrix.hpp"
112#include "Tpetra_Vector.hpp"
113#include "Thyra_TpetraThyraWrappers.hpp"
114#include "TpetraExt_MatrixMatrix.hpp"
115#include "Tpetra_RowMatrixTransposer.hpp"
116
117#include <cmath>
118
119namespace Teko {
120
121using Teuchos::rcp;
122using Teuchos::rcp_dynamic_cast;
123using Teuchos::RCP;
124#ifdef Teko_ENABLE_Isorropia
125using Isorropia::Epetra::Prober;
126#endif
127
128const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
129{
130 Teuchos::RCP<Teuchos::FancyOStream> os =
131 Teuchos::VerboseObjectBase::getDefaultOStream();
132
133 //os->setShowProcRank(true);
134 //os->setOutputToRootOnly(-1);
135 return os;
136}
137
138// distance function...not parallel...entirely internal to this cpp file
139inline double dist(int dim,double * coords,int row,int col)
140{
141 double value = 0.0;
142 for(int i=0;i<dim;i++)
143 value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
144
145 // the distance between the two
146 return std::sqrt(value);
147}
148
149// distance function...not parallel...entirely internal to this cpp file
150inline double dist(double * x,double * y,double * z,int stride,int row,int col)
151{
152 double value = 0.0;
153 if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
154 if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
155 if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
156
157 // the distance between the two
158 return std::sqrt(value);
159}
160
179#ifdef TEKO_HAVE_EPETRA
180RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil)
181{
182 // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
183 RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
184 stencil.MaxNumEntries()+1),true);
185
186 // allocate an additional value for the diagonal, if neccessary
187 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
188 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
189
190 // loop over all the rows
191 for(int j=0;j<gl->NumMyRows();j++) {
192 int row = gl->GRID(j);
193 double diagValue = 0.0;
194 int diagInd = -1;
195 int rowSz = 0;
196
197 // extract a copy of this row...put it in rowData, rowIndicies
198 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
199
200 // loop over elements of row
201 for(int i=0;i<rowSz;i++) {
202 int col = rowInd[i];
203
204 // is this a 0 entry masquerading as some thing else?
205 double value = rowData[i];
206 if(value==0) continue;
207
208 // for nondiagonal entries
209 if(row!=col) {
210 double d = dist(dim,coords,row,col);
211 rowData[i] = -1.0/d;
212 diagValue += rowData[i];
213 }
214 else
215 diagInd = i;
216 }
217
218 // handle diagonal entry
219 if(diagInd<0) { // diagonal not in row
220 rowData[rowSz] = -diagValue;
221 rowInd[rowSz] = row;
222 rowSz++;
223 }
224 else { // diagonal in row
225 rowData[diagInd] = -diagValue;
226 rowInd[diagInd] = row;
227 }
228
229 // insert row data into graph Laplacian matrix
230 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
231 }
232
233 gl->FillComplete();
234
235 return gl;
236}
237#endif
238
239RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(int dim,ST * coords,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
240{
241 // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
242 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
243 stencil.getGlobalMaxNumRowEntries()+1));
244
245 // allocate an additional value for the diagonal, if neccessary
246 auto rowInd = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
247 auto rowData = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
248
249 // loop over all the rows
250 for(LO j=0;j< (LO) gl->getLocalNumRows();j++) {
251 GO row = gl->getRowMap()->getGlobalElement(j);
252 ST diagValue = 0.0;
253 GO diagInd = -1;
254 size_t rowSz = 0;
255
256 // extract a copy of this row...put it in rowData, rowIndicies
257 stencil.getGlobalRowCopy(row,rowInd,rowData,rowSz);
258
259 // loop over elements of row
260 for(size_t i=0;i<rowSz;i++) {
261 GO col = rowInd(i);
262
263 // is this a 0 entry masquerading as some thing else?
264 ST value = rowData[i];
265 if(value==0) continue;
266
267 // for nondiagonal entries
268 if(row!=col) {
269 ST d = dist(dim,coords,row,col);
270 rowData[i] = -1.0/d;
271 diagValue += rowData(i);
272 }
273 else
274 diagInd = i;
275 }
276
277 // handle diagonal entry
278 if(diagInd<0) { // diagonal not in row
279 rowData(rowSz) = -diagValue;
280 rowInd(rowSz) = row;
281 rowSz++;
282 }
283 else { // diagonal in row
284 rowData(diagInd) = -diagValue;
285 rowInd(diagInd) = row;
286 }
287
288 // insert row data into graph Laplacian matrix
289 gl->replaceGlobalValues(row,rowInd,rowData);
290 }
291
292 gl->fillComplete();
293
294 return gl;
295}
296
319#ifdef TEKO_HAVE_EPETRA
320RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil)
321{
322 // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
323 RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
324 stencil.MaxNumEntries()+1),true);
325
326 // allocate an additional value for the diagonal, if neccessary
327 std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
328 std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
329
330 // loop over all the rows
331 for(int j=0;j<gl->NumMyRows();j++) {
332 int row = gl->GRID(j);
333 double diagValue = 0.0;
334 int diagInd = -1;
335 int rowSz = 0;
336
337 // extract a copy of this row...put it in rowData, rowIndicies
338 stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
339
340 // loop over elements of row
341 for(int i=0;i<rowSz;i++) {
342 int col = rowInd[i];
343
344 // is this a 0 entry masquerading as some thing else?
345 double value = rowData[i];
346 if(value==0) continue;
347
348 // for nondiagonal entries
349 if(row!=col) {
350 double d = dist(x,y,z,stride,row,col);
351 rowData[i] = -1.0/d;
352 diagValue += rowData[i];
353 }
354 else
355 diagInd = i;
356 }
357
358 // handle diagonal entry
359 if(diagInd<0) { // diagonal not in row
360 rowData[rowSz] = -diagValue;
361 rowInd[rowSz] = row;
362 rowSz++;
363 }
364 else { // diagonal in row
365 rowData[diagInd] = -diagValue;
366 rowInd[diagInd] = row;
367 }
368
369 // insert row data into graph Laplacian matrix
370 TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
371 }
372
373 gl->FillComplete();
374
375 return gl;
376}
377#endif
378
379RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
380{
381 // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
382 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
383 stencil.getGlobalMaxNumRowEntries()+1),true);
384
385 // allocate an additional value for the diagonal, if neccessary
386 auto rowInd = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_global_inds_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
387 auto rowData = typename Tpetra::CrsMatrix<ST,LO,GO,NT>::nonconst_values_host_view_type(Kokkos::ViewAllocateWithoutInitializing("rowIndices"),stencil.getGlobalMaxNumRowEntries()+1);
388
389 // loop over all the rows
390 for(LO j=0;j<(LO) gl->getLocalNumRows();j++) {
391 GO row = gl->getRowMap()->getGlobalElement(j);
392 ST diagValue = 0.0;
393 GO diagInd = -1;
394 size_t rowSz = 0;
395
396 // extract a copy of this row...put it in rowData, rowIndicies
397 stencil.getGlobalRowCopy(row,rowInd,rowData,rowSz);
398
399 // loop over elements of row
400 for(size_t i=0;i<rowSz;i++) {
401 GO col = rowInd(i);
402
403 // is this a 0 entry masquerading as some thing else?
404 ST value = rowData[i];
405 if(value==0) continue;
406
407 // for nondiagonal entries
408 if(row!=col) {
409 ST d = dist(x,y,z,stride,row,col);
410 rowData[i] = -1.0/d;
411 diagValue += rowData(i);
412 }
413 else
414 diagInd = i;
415 }
416
417 // handle diagonal entry
418 if(diagInd<0) { // diagonal not in row
419 rowData(rowSz) = -diagValue;
420 rowInd(rowSz) = row;
421 rowSz++;
422 }
423 else { // diagonal in row
424 rowData(diagInd) = -diagValue;
425 rowInd(diagInd) = row;
426 }
427
428 // insert row data into graph Laplacian matrix
429 gl->replaceGlobalValues(row,rowInd,rowData);
430 }
431
432 gl->fillComplete();
433
434 return gl;
435}
436
452void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
453{
454 Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
455}
456
472void applyTransposeOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
473{
474 Thyra::apply(*A,Thyra::TRANS,*x,y.ptr(),alpha,beta);
475}
476
479void update(double alpha,const MultiVector & x,double beta,MultiVector & y)
480{
481 Teuchos::Array<double> scale;
482 Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > > vec;
483
484 // build arrays needed for linear combo
485 scale.push_back(alpha);
486 vec.push_back(x.ptr());
487
488 // compute linear combination
489 Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
490}
491
493BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
494{
495 int rows = blockRowCount(blo);
496
497 TEUCHOS_ASSERT(rows==blockColCount(blo));
498
499 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
500 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
501
502 // allocate new operator
503 BlockedLinearOp upper = createBlockedOp();
504
505 // build new operator
506 upper->beginBlockFill(rows,rows);
507
508 for(int i=0;i<rows;i++) {
509 // put zero operators on the diagonal
510 // this gurantees the vector space of
511 // the new operator are fully defined
512 RCP<const Thyra::LinearOpBase<double> > zed
513 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
514 upper->setBlock(i,i,zed);
515
516 for(int j=i+1;j<rows;j++) {
517 // get block i,j
518 LinearOp uij = blo->getBlock(i,j);
519
520 // stuff it in U
521 if(uij!=Teuchos::null)
522 upper->setBlock(i,j,uij);
523 }
524 }
525 if(callEndBlockFill)
526 upper->endBlockFill();
527
528 return upper;
529}
530
532BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
533{
534 int rows = blockRowCount(blo);
535
536 TEUCHOS_ASSERT(rows==blockColCount(blo));
537
538 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
539 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
540
541 // allocate new operator
542 BlockedLinearOp lower = createBlockedOp();
543
544 // build new operator
545 lower->beginBlockFill(rows,rows);
546
547 for(int i=0;i<rows;i++) {
548 // put zero operators on the diagonal
549 // this gurantees the vector space of
550 // the new operator are fully defined
551 RCP<const Thyra::LinearOpBase<double> > zed
552 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
553 lower->setBlock(i,i,zed);
554
555 for(int j=0;j<i;j++) {
556 // get block i,j
557 LinearOp uij = blo->getBlock(i,j);
558
559 // stuff it in U
560 if(uij!=Teuchos::null)
561 lower->setBlock(i,j,uij);
562 }
563 }
564 if(callEndBlockFill)
565 lower->endBlockFill();
566
567 return lower;
568}
569
589BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo)
590{
591 int rows = blockRowCount(blo);
592
593 TEUCHOS_ASSERT(rows==blockColCount(blo)); // assert that matrix is square
594
595 RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
596 RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
597
598 // allocate new operator
599 BlockedLinearOp zeroOp = createBlockedOp();
600
601 // build new operator
602 zeroOp->beginBlockFill(rows,rows);
603
604 for(int i=0;i<rows;i++) {
605 // put zero operators on the diagonal
606 // this gurantees the vector space of
607 // the new operator are fully defined
608 RCP<const Thyra::LinearOpBase<double> > zed
609 = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
610 zeroOp->setBlock(i,i,zed);
611 }
612
613 return zeroOp;
614}
615
617bool isZeroOp(const LinearOp op)
618{
619 // if operator is null...then its zero!
620 if(op==Teuchos::null) return true;
621
622 // try to cast it to a zero linear operator
623 LinearOp test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(op);
624
625 // if it works...then its zero...otherwise its null
626 if(test!=Teuchos::null) return true;
627
628 // See if the operator is a wrapped zero op
629 ST scalar = 0.0;
630 Thyra::EOpTransp transp = Thyra::NOTRANS;
631 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
632 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
633 test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(wrapped_op);
634 return test!=Teuchos::null;
635}
636
637std::pair<ModifiableLinearOp, bool> getAbsRowSumMatrixEpetra(const LinearOp &op) {
638#ifndef TEKO_HAVE_EPETRA
639 return std::make_pair(ModifiableLinearOp{}, false);
640#else
641 RCP<const Epetra_CrsMatrix> eCrsOp;
642
643 const auto eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op);
644
645 if(!eOp) {
646 return std::make_pair(ModifiableLinearOp{}, false);
647 }
648
649 eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(), true);
650
651
652 // extract diagonal
653 const auto ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
654 Epetra_Vector &diag = *ptrDiag;
655
656 // compute absolute value row sum
657 diag.PutScalar(0.0);
658 for (int i = 0; i < eCrsOp->NumMyRows(); i++) {
659 double *values = 0;
660 int numEntries;
661 eCrsOp->ExtractMyRowView(i, numEntries, values);
662
663 // build abs value row sum
664 for (int j = 0; j < numEntries; j++)
665 diag[i] += std::abs(values[j]);
666 }
667
668 // build Thyra diagonal operator
669 return std::make_pair(
670 Teko::Epetra::thyraDiagOp(ptrDiag, eCrsOp->RowMap(),
671 "absRowSum( " + op->getObjectLabel() + " )"),
672 true);
673#endif
674}
675
676std::pair<ModifiableLinearOp, bool> getAbsRowSumMatrixTpetra(const LinearOp &op) {
677 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT>> tCrsOp;
678
679 const auto tOp =
680 rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT>>(op);
681
682 tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST, LO, GO, NT>>(
683 tOp->getConstTpetraOperator(), true);
684
685 // extract diagonal
686 const auto ptrDiag =
687 Tpetra::createVector<ST, LO, GO, NT>(tCrsOp->getRowMap());
688 auto &diag = *ptrDiag;
689
690 // compute absolute value row sum
691 diag.putScalar(0.0);
692 for (LO i = 0; i < (LO)tCrsOp->getLocalNumRows(); i++) {
693 auto numEntries = tCrsOp->getNumEntriesInLocalRow(i);
694 typename Tpetra::CrsMatrix<ST, LO, GO, NT>::local_inds_host_view_type
695 indices;
696 typename Tpetra::CrsMatrix<ST, LO, GO, NT>::values_host_view_type values;
697 tCrsOp->getLocalRowView(i, indices, values);
698
699 // build abs value row sum
700 for (size_t j = 0; j < numEntries; j++)
701 diag.sumIntoLocalValue(i, std::abs(values(j)));
702 }
703
704 // build Thyra diagonal operator
705 return std::make_pair(Teko::TpetraHelpers::thyraDiagOp(
706 ptrDiag, *tCrsOp->getRowMap(),
707 "absRowSum( " + op->getObjectLabel() + " ))"),
708 true);
709}
710
719ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op)
720{
721 try {
722 auto eResult = getAbsRowSumMatrixEpetra(op);
723 if (eResult.second) {
724 return eResult.first;
725 }
726
727 auto tResult = getAbsRowSumMatrixTpetra(op);
728 if (tResult.second) {
729 return tResult.first;
730 } else {
731 throw std::logic_error("Neither Epetra nor Tpetra");
732 }
733 } catch (std::exception &e) {
734 auto out = Teuchos::VerboseObjectBase::getDefaultOStream();
735
736 *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix or a "
737 "Tpetra::CrsMatrix\n";
738 *out << " Could not extract an Epetra_Operator or a Tpetra_Operator "
739 "from a \""
740 << op->description() << std::endl;
741 *out << " OR\n";
742 *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or "
743 "a Tpetra_Operator to a Tpetra::CrsMatrix\n";
744 *out << std::endl;
745 *out << "*** THROWN EXCEPTION ***\n";
746 *out << e.what() << std::endl;
747 *out << "************************\n";
748
749 throw e;
750 }
751}
752
761ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op)
762{
763 // if this is a blocked operator, extract diagonals block by block
764 // FIXME: this does not add in values from off-diagonal blocks
765 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_op = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
766 if(blocked_op != Teuchos::null){
767 int numRows = blocked_op->productRange()->numBlocks();
768 TEUCHOS_ASSERT(blocked_op->productDomain()->numBlocks() == numRows);
769 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_diag = Thyra::defaultBlockedLinearOp<double>();
770 blocked_diag->beginBlockFill(numRows,numRows);
771 for(int r = 0; r < numRows; ++r){
772 for(int c = 0; c < numRows; ++c){
773 if(r==c)
774 blocked_diag->setNonconstBlock(r,c,getAbsRowSumInvMatrix(blocked_op->getBlock(r,c)));
775 else
776 blocked_diag->setBlock(r,c,Thyra::zero<double>(blocked_op->getBlock(r,c)->range(),blocked_op->getBlock(r,c)->domain()));
777 }
778 }
779 blocked_diag->endBlockFill();
780 return blocked_diag;
781 }
782
783 if(Teko::TpetraHelpers::isTpetraLinearOp(op)) {
784 ST scalar = 0.0;
785 bool transp = false;
786 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
787
788 // extract diagonal
789 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
790 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
791
792 // compute absolute value row sum
793 diag.putScalar(0.0);
794 for(LO i=0;i<(LO) tCrsOp->getLocalNumRows();i++) {
795 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
796 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::local_inds_host_view_type indices;
797 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::values_host_view_type values;
798 tCrsOp->getLocalRowView(i,indices,values);
799
800 // build abs value row sum
801 for(LO j=0;j<numEntries;j++)
802 diag.sumIntoLocalValue(i,std::abs(values(j)));
803 }
804 diag.scale(scalar);
805 diag.reciprocal(diag); // invert entries
806
807 // build Thyra diagonal operator
808 return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),"absRowSum( " + op->getObjectLabel() + " ))");
809
810 }
811 else{
812#ifdef TEKO_HAVE_EPETRA
813 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op,true);
814 RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
815
816 // extract diagonal
817 const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
818 Epetra_Vector & diag = *ptrDiag;
819
820 // compute absolute value row sum
821 diag.PutScalar(0.0);
822 for(int i=0;i<eCrsOp->NumMyRows();i++) {
823 double * values = 0;
824 int numEntries;
825 eCrsOp->ExtractMyRowView(i,numEntries,values);
826
827 // build abs value row sum
828 for(int j=0;j<numEntries;j++)
829 diag[i] += std::abs(values[j]);
830 }
831 diag.Reciprocal(diag); // invert entries
832
833 // build Thyra diagonal operator
834 return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
835#else
836 throw std::logic_error("getAbsRowSumInvMatrix is trying to use Epetra "
837 "code, but TEKO_HAVE_EPETRA is disabled!");
838#endif
839 }
840
841}
842
850ModifiableLinearOp getLumpedMatrix(const LinearOp & op)
851{
852 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
853 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
854
855 // set to all ones
856 Thyra::assign(ones.ptr(),1.0);
857
858 // compute lumped diagonal
859 // Thyra::apply(*op,Thyra::NONCONJ_ELE,*ones,&*diag);
860 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
861
862 return rcp(new Thyra::DefaultDiagonalLinearOp<ST>(diag));
863}
864
873ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op)
874{
875 RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
876 RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
877
878 // set to all ones
879 Thyra::assign(ones.ptr(),1.0);
880
881 // compute lumped diagonal
882 Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
883 Thyra::reciprocal(*diag,diag.ptr());
884
885 return rcp(new Thyra::DefaultDiagonalLinearOp<ST>(diag));
886}
887
888const std::pair<ModifiableLinearOp, bool> getDiagonalOpEpetra(const LinearOp &op) {
889#ifndef TEKO_HAVE_EPETRA
890 return std::make_pair(ModifiableLinearOp{}, false);
891#else
892 RCP<const Epetra_CrsMatrix> eCrsOp;
893
894 const auto eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp>(op);
895 if (!eOp) {
896 return std::make_pair(ModifiableLinearOp{}, false);
897 }
898
899 eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(), true);
900
901 // extract diagonal
902 const auto diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
903 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
904
905 // build Thyra diagonal operator
906 return std::make_pair(
907 Teko::Epetra::thyraDiagOp(diag, eCrsOp->RowMap(),
908 "inv(diag( " + op->getObjectLabel() + " ))"),
909 true);
910#endif
911}
912
913const std::pair<ModifiableLinearOp, bool> getDiagonalOpTpetra(const LinearOp &op) {
914 RCP<const Tpetra::CrsMatrix<ST, LO, GO, NT>> tCrsOp;
915
916 const auto tOp =
917 rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT>>(op);
918 if (!tOp) {
919 return std::make_pair(ModifiableLinearOp{}, false);
920 }
921
922 tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST, LO, GO, NT>>(
923 tOp->getConstTpetraOperator(), true);
924
925 // extract diagonal
926 const auto diag = Tpetra::createVector<ST, LO, GO, NT>(tCrsOp->getRowMap());
927 tCrsOp->getLocalDiagCopy(*diag);
928
929 // build Thyra diagonal operator
930 return std::make_pair(Teko::TpetraHelpers::thyraDiagOp(
931 diag, *tCrsOp->getRowMap(),
932 "inv(diag( " + op->getObjectLabel() + " ))"),
933 true);
934}
935
947const ModifiableLinearOp getDiagonalOp(const LinearOp & op)
948{
949 try {
950 // get Epetra or Tpetra Operator
951 const auto eDiagOp = getDiagonalOpEpetra(op);
952
953 if (eDiagOp.second) {
954 return eDiagOp.first;
955 }
956
957 const auto tDiagOp = getDiagonalOpTpetra(op);
958 if (tDiagOp.second) {
959 return tDiagOp.first;
960 }
961 else
962 throw std::logic_error("Neither Epetra nor Tpetra");
963 }
964 catch (std::exception & e) {
965 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
966
967 *out << "Teko: getDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
968 *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
969 *out << " OR\n";
970 *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
971 *out << std::endl;
972 *out << "*** THROWN EXCEPTION ***\n";
973 *out << e.what() << std::endl;
974 *out << "************************\n";
975
976 throw e;
977 }
978}
979
980const MultiVector getDiagonal(const LinearOp & op)
981{
982 try {
983 // get Epetra or Tpetra Operator
984 auto diagOp = getDiagonalOpEpetra(op);
985
986 if (!diagOp.second) {
987 diagOp = getDiagonalOpTpetra(op);
988
989 if(!diagOp.second){
990 throw std::logic_error("Neither Epetra nor Tpetra");
991 }
992 }
993
994 Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
995 Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp.first)->getDiag();
996 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
997 }
998 catch (std::exception & e) {
999 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
1000
1001 *out << "Teko: getDiagonal requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
1002 *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
1003 *out << " OR\n";
1004 *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
1005 *out << std::endl;
1006 *out << "*** THROWN EXCEPTION ***\n";
1007 *out << e.what() << std::endl;
1008 *out << "************************\n";
1009
1010 throw e;
1011 }
1012}
1013
1014const MultiVector getDiagonal(const Teko::LinearOp & A,const DiagonalType & dt)
1015{
1016 LinearOp diagOp = Teko::getDiagonalOp(A,dt);
1017
1018 Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
1019 Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag();
1020 return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
1021}
1022
1034const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op)
1035{
1036 // if this is a diagonal linear op already, just take the reciprocal
1037 auto diagonal_op = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double>>(op);
1038 if(diagonal_op != Teuchos::null){
1039 auto diag = diagonal_op->getDiag();
1040 auto inv_diag = diag->clone_v();
1041 Thyra::reciprocal(*diag,inv_diag.ptr());
1042 return rcp(new Thyra::DefaultDiagonalLinearOp<double>(inv_diag));
1043 }
1044
1045 // if this is a blocked operator, extract diagonals block by block
1046 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_op = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
1047 if(blocked_op != Teuchos::null){
1048 int numRows = blocked_op->productRange()->numBlocks();
1049 TEUCHOS_ASSERT(blocked_op->productDomain()->numBlocks() == numRows);
1050 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_diag = Thyra::defaultBlockedLinearOp<double>();
1051 blocked_diag->beginBlockFill(numRows,numRows);
1052 for(int r = 0; r < numRows; ++r){
1053 for(int c = 0; c < numRows; ++c){
1054 if(r==c)
1055 blocked_diag->setNonconstBlock(r,c,getInvDiagonalOp(blocked_op->getBlock(r,c)));
1056 else
1057 blocked_diag->setBlock(r,c,Thyra::zero<double>(blocked_op->getBlock(r,c)->range(),blocked_op->getBlock(r,c)->domain()));
1058 }
1059 }
1060 blocked_diag->endBlockFill();
1061 return blocked_diag;
1062 }
1063
1064 if (Teko::TpetraHelpers::isTpetraLinearOp(op)){
1065 ST scalar = 0.0;
1066 bool transp = false;
1067 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
1068
1069 // extract diagonal
1070 const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
1071 diag->scale(scalar);
1072 tCrsOp->getLocalDiagCopy(*diag);
1073 diag->reciprocal(*diag);
1074
1075 // build Thyra diagonal operator
1076 return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
1077
1078 }
1079 else {
1080#ifdef TEKO_HAVE_EPETRA
1081 RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op,true);
1082 RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
1083
1084 // extract diagonal
1085 const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
1086 TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
1087 diag->Reciprocal(*diag);
1088
1089 // build Thyra diagonal operator
1090 return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
1091#else
1092 throw std::logic_error("getInvDiagonalOp is trying to use Epetra "
1093 "code, but TEKO_HAVE_EPETRA is disabled!");
1094#endif
1095 }
1096}
1097
1110const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr)
1111{
1112 // if this is a blocked operator, multiply block by block
1113 // it is possible that not every factor in the product is blocked and these situations are handled separately
1114
1115 bool isBlockedL = isPhysicallyBlockedLinearOp(opl);
1116 bool isBlockedM = isPhysicallyBlockedLinearOp(opm);
1117 bool isBlockedR = isPhysicallyBlockedLinearOp(opr);
1118
1119 // all factors blocked
1120 if((isBlockedL && isBlockedM && isBlockedR)){
1121
1122 double scalarl = 0.0;
1123 bool transpl = false;
1124 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1125 double scalarm = 0.0;
1126 bool transpm = false;
1127 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opm = getPhysicallyBlockedLinearOp(opm,&scalarm,&transpm);
1128 double scalarr = 0.0;
1129 bool transpr = false;
1130 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1131 double scalar = scalarl*scalarm*scalarr;
1132
1133 int numRows = blocked_opl->productRange()->numBlocks();
1134 int numCols = blocked_opr->productDomain()->numBlocks();
1135 int numMiddle = blocked_opm->productRange()->numBlocks();
1136
1137 // Assume that the middle block is block nxn and that it's diagonal. Otherwise use the two argument explicitMultiply twice
1138 TEUCHOS_ASSERT(blocked_opm->productDomain()->numBlocks() == numMiddle);
1139 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1140 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1141
1142 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1143 blocked_product->beginBlockFill(numRows,numCols);
1144 for(int r = 0; r < numRows; ++r){
1145 for(int c = 0; c < numCols; ++c){
1146 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),blocked_opm->getBlock(0,0),blocked_opr->getBlock(0,c));
1147 for(int m = 1; m < numMiddle; ++m){
1148 LinearOp product_m = explicitMultiply(blocked_opl->getBlock(r,m),blocked_opm->getBlock(m,m),blocked_opr->getBlock(m,c));
1149 product_rc = explicitAdd(product_rc,product_m);
1150 }
1151 blocked_product->setBlock(r,c,product_rc);
1152 }
1153 }
1154 blocked_product->endBlockFill();
1155 return Thyra::scale<double>(scalar,blocked_product.getConst());
1156 }
1157
1158 // left and right factors blocked
1159 if(isBlockedL && !isBlockedM && isBlockedR){
1160 double scalarl = 0.0;
1161 bool transpl = false;
1162 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1163 double scalarr = 0.0;
1164 bool transpr = false;
1165 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1166 double scalar = scalarl*scalarr;
1167
1168 int numRows = blocked_opl->productRange()->numBlocks();
1169 int numCols = blocked_opr->productDomain()->numBlocks();
1170 int numMiddle = 1;
1171
1172 // Assume that the middle block is 1x1 diagonal. Left must be rx1, right 1xc
1173 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1174 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1175
1176 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1177 blocked_product->beginBlockFill(numRows,numCols);
1178 for(int r = 0; r < numRows; ++r){
1179 for(int c = 0; c < numCols; ++c){
1180 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),opm,blocked_opr->getBlock(0,c));
1181 blocked_product->setBlock(r,c,product_rc);
1182 }
1183 }
1184 blocked_product->endBlockFill();
1185 return Thyra::scale<double>(scalar,blocked_product.getConst());
1186 }
1187
1188 // only right factor blocked
1189 if(!isBlockedL && !isBlockedM && isBlockedR){
1190 double scalarr = 0.0;
1191 bool transpr = false;
1192 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1193 double scalar = scalarr;
1194
1195 int numRows = 1;
1196 int numCols = blocked_opr->productDomain()->numBlocks();
1197 int numMiddle = 1;
1198
1199 // Assume that the middle block is 1x1 diagonal, left is 1x1. Right must be 1xc
1200 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1201
1202 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1203 blocked_product->beginBlockFill(numRows,numCols);
1204 for(int c = 0; c < numCols; ++c){
1205 LinearOp product_c = explicitMultiply(opl,opm,blocked_opr->getBlock(0,c));
1206 blocked_product->setBlock(0,c,product_c);
1207 }
1208 blocked_product->endBlockFill();
1209 return Thyra::scale<double>(scalar,blocked_product.getConst());
1210 }
1211
1212 //TODO: three more cases (only non-blocked - blocked - non-blocked not possible)
1213
1214 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1215 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1216 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1217
1218 if(isTpetral && isTpetram && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1219
1220 // Get left and right Tpetra crs operators
1221 ST scalarl = 0.0;
1222 bool transpl = false;
1223 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1224 ST scalarm = 0.0;
1225 bool transpm = false;
1226 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1227 ST scalarr = 0.0;
1228 bool transpr = false;
1229 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1230
1231 // Build output operator
1232 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1233 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1234
1235 // Do explicit matrix-matrix multiply
1236 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1237 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1238 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1239 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1240 explicitCrsOp->resumeFill();
1241 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1242 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1243 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1244 return tExplicitOp;
1245
1246 } else if (isTpetral && !isTpetram && isTpetrar){ // Assume that the middle operator is diagonal
1247
1248 // Get left and right Tpetra crs operators
1249 ST scalarl = 0.0;
1250 bool transpl = false;
1251 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1252 ST scalarr = 0.0;
1253 bool transpr = false;
1254 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1255
1256 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr;
1257
1258 // Cast middle operator as DiagonalLinearOp and extract diagonal as Vector
1259 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm);
1260 if(dOpm != Teuchos::null){
1261 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),true);
1262 diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1263 }
1264 // If it's not diagonal, maybe it's zero
1265 else if(rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<ST> >(opm) != Teuchos::null){
1266 diagPtr = rcp(new Tpetra::Vector<ST,LO,GO,NT>(tCrsOpl->getDomainMap()));
1267 }
1268 else
1269 TEUCHOS_ASSERT(false);
1270
1271 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1272
1273 // Do the diagonal scaling
1274 tCrsOplm->rightScale(*diagPtr);
1275
1276 // Build output operator
1277 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1278 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1279
1280 // Do explicit matrix-matrix multiply
1281 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1282 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1283 explicitCrsOp->resumeFill();
1284 explicitCrsOp->scale(scalarl*scalarr);
1285 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1286 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1287 return tExplicitOp;
1288
1289 } else { // Assume Epetra and we can use transformers
1290#ifdef TEKO_HAVE_EPETRA
1291 // build implicit multiply
1292 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1293
1294 // build transformer
1295 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1296 Thyra::epetraExtDiagScaledMatProdTransformer();
1297
1298 // build operator and multiply
1299 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1300 prodTrans->transform(*implicitOp,explicitOp.ptr());
1301 explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1302 " * " + opm->getObjectLabel() +
1303 " * " + opr->getObjectLabel() + " )");
1304
1305 return explicitOp;
1306#else
1307 throw std::logic_error("explicitMultiply is trying to use Epetra "
1308 "code, but TEKO_HAVE_EPETRA is disabled!");
1309#endif
1310 }
1311}
1312
1327const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
1328 const ModifiableLinearOp & destOp)
1329{
1330 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1331 bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1332 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1333
1334 if(isTpetral && isTpetram && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1335
1336 // Get left and right Tpetra crs operators
1337 ST scalarl = 0.0;
1338 bool transpl = false;
1339 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1340 ST scalarm = 0.0;
1341 bool transpm = false;
1342 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opm, &scalarm, &transpm);
1343 ST scalarr = 0.0;
1344 bool transpr = false;
1345 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1346
1347 // Build output operator
1348 auto tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(destOp);
1349 if(tExplicitOp.is_null())
1350 tExplicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1351
1352 // Do explicit matrix-matrix multiply
1353 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1354 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1355 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1356 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1357 explicitCrsOp->resumeFill();
1358 explicitCrsOp->scale(scalarl*scalarm*scalarr);
1359 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1360 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1361 return tExplicitOp;
1362
1363 } else if (isTpetral && !isTpetram && isTpetrar){ // Assume that the middle operator is diagonal
1364
1365 // Get left and right Tpetra crs operators
1366 ST scalarl = 0.0;
1367 bool transpl = false;
1368 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1369 ST scalarr = 0.0;
1370 bool transpr = false;
1371 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1372
1373 // Cast middle operator as DiagonalLinearOp and extract diagonal as Vector
1374 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm,true);
1375 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),true);
1376 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1377 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1378
1379 // Do the diagonal scaling
1380 tCrsOplm->rightScale(*diagPtr);
1381
1382 // Build output operator
1383 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1384 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1385
1386 // Do explicit matrix-matrix multiply
1387 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1388 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1389 explicitCrsOp->resumeFill();
1390 explicitCrsOp->scale(scalarl*scalarr);
1391 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1392 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1393 return tExplicitOp;
1394
1395 } else { // Assume Epetra and we can use transformers
1396#ifdef TEKO_HAVE_EPETRA
1397 // build implicit multiply
1398 const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1399
1400 // build transformer
1401 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1402 Thyra::epetraExtDiagScaledMatProdTransformer();
1403
1404 // build operator destination operator
1405 ModifiableLinearOp explicitOp;
1406
1407 // if neccessary build a operator to put the explicit multiply into
1408 if(destOp==Teuchos::null)
1409 explicitOp = prodTrans->createOutputOp();
1410 else
1411 explicitOp = destOp;
1412
1413 // perform multiplication
1414 prodTrans->transform(*implicitOp,explicitOp.ptr());
1415
1416 // label it
1417 explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1418 " * " + opm->getObjectLabel() +
1419 " * " + opr->getObjectLabel() + " )");
1420
1421 return explicitOp;
1422#else
1423 throw std::logic_error("explicitMultiply is trying to use Epetra "
1424 "code, but TEKO_HAVE_EPETRA is disabled!");
1425#endif
1426 }
1427}
1428
1439const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr)
1440{
1441 // if this is a blocked operator, multiply block by block
1442 // it is possible that not every factor in the product is blocked and these situations are handled separately
1443
1444 bool isBlockedL = isPhysicallyBlockedLinearOp(opl);
1445 bool isBlockedR = isPhysicallyBlockedLinearOp(opr);
1446
1447 // both factors blocked
1448 if((isBlockedL && isBlockedR)){
1449 double scalarl = 0.0;
1450 bool transpl = false;
1451 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1452 double scalarr = 0.0;
1453 bool transpr = false;
1454 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1455 double scalar = scalarl*scalarr;
1456
1457 int numRows = blocked_opl->productRange()->numBlocks();
1458 int numCols = blocked_opr->productDomain()->numBlocks();
1459 int numMiddle = blocked_opl->productDomain()->numBlocks();
1460
1461 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1462
1463 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1464 blocked_product->beginBlockFill(numRows,numCols);
1465 for(int r = 0; r < numRows; ++r){
1466 for(int c = 0; c < numCols; ++c){
1467 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),blocked_opr->getBlock(0,c));
1468 for(int m = 1; m < numMiddle; ++m){
1469 LinearOp product_m = explicitMultiply(blocked_opl->getBlock(r,m),blocked_opr->getBlock(m,c));
1470 product_rc = explicitAdd(product_rc,product_m);
1471 }
1472 blocked_product->setBlock(r,c,Thyra::scale(scalar,product_rc));
1473 }
1474 }
1475 blocked_product->endBlockFill();
1476 return blocked_product;
1477 }
1478
1479 // only left factor blocked
1480 if((isBlockedL && !isBlockedR)){
1481 double scalarl = 0.0;
1482 bool transpl = false;
1483 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1484 double scalar = scalarl;
1485
1486 int numRows = blocked_opl->productRange()->numBlocks();
1487 int numCols = 1;
1488 int numMiddle = 1;
1489
1490 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1491
1492 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1493 blocked_product->beginBlockFill(numRows,numCols);
1494 for(int r = 0; r < numRows; ++r){
1495 LinearOp product_r = explicitMultiply(blocked_opl->getBlock(r,0),opr);
1496 blocked_product->setBlock(r,0,Thyra::scale(scalar,product_r));
1497 }
1498 blocked_product->endBlockFill();
1499 return blocked_product;
1500 }
1501
1502 // only right factor blocked
1503 if((!isBlockedL && isBlockedR)){
1504 double scalarr = 0.0;
1505 bool transpr = false;
1506 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1507 double scalar = scalarr;
1508
1509 int numRows = 1;
1510 int numCols = blocked_opr->productDomain()->numBlocks();
1511 int numMiddle = 1;
1512
1513 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1514
1515 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1516 blocked_product->beginBlockFill(numRows,numCols);
1517 for(int c = 0; c < numCols; ++c){
1518 LinearOp product_c = explicitMultiply(opl,blocked_opr->getBlock(0,c));
1519 blocked_product->setBlock(0,c,Thyra::scale(scalar,product_c));
1520 }
1521 blocked_product->endBlockFill();
1522 return blocked_product;
1523 }
1524
1525 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1526 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1527
1528 if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1529 // Get left and right Tpetra crs operators
1530 ST scalarl = 0.0;
1531 bool transpl = false;
1532 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1533 ST scalarr = 0.0;
1534 bool transpr = false;
1535 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1536
1537 // Build output operator
1538 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1539 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1540
1541 // Do explicit matrix-matrix multiply
1542 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1543 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1544 explicitCrsOp->resumeFill();
1545 explicitCrsOp->scale(scalarl*scalarr);
1546 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1547 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1548 return tExplicitOp;
1549
1550 } else if (isTpetral && !isTpetrar){ // Assume that the right operator is diagonal
1551
1552 // Get left Tpetra crs operator
1553 ST scalarl = 0.0;
1554 bool transpl = false;
1555 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1556
1557 // Cast right operator as DiagonalLinearOp and extract diagonal as Vector
1558 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr,true);
1559 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),true);
1560 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1561 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1562
1563 explicitCrsOp->rightScale(*diagPtr);
1564 explicitCrsOp->resumeFill();
1565 explicitCrsOp->scale(scalarl);
1566 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1567
1568 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1569
1570 } else if (!isTpetral && isTpetrar){ // Assume that the left operator is diagonal
1571
1572 // Get right Tpetra crs operator
1573 ST scalarr = 0.0;
1574 bool transpr = false;
1575 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1576
1577 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr;
1578
1579 // Cast left operator as DiagonalLinearOp and extract diagonal as Vector
1580 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl);
1581 if(dOpl != Teuchos::null){
1582 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),true);
1583 diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1584 }
1585 // If it's not diagonal, maybe it's zero
1586 else if(rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<ST> >(opl) != Teuchos::null){
1587 diagPtr = rcp(new Tpetra::Vector<ST,LO,GO,NT>(tCrsOpr->getRangeMap()));
1588 }
1589 else
1590 TEUCHOS_ASSERT(false);
1591
1592 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpr, Tpetra::Import<LO,GO,NT>(tCrsOpr->getRowMap(),tCrsOpr->getRowMap()));
1593
1594 explicitCrsOp->leftScale(*diagPtr);
1595 explicitCrsOp->resumeFill();
1596 explicitCrsOp->scale(scalarr);
1597 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1598
1599 return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1600
1601 } else { // Assume Epetra and we can use transformers
1602#ifdef TEKO_HAVE_EPETRA
1603 // build implicit multiply
1604 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1605
1606 // build a scaling transformer
1607 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1608 = Thyra::epetraExtDiagScalingTransformer();
1609
1610 // check to see if a scaling transformer works: if not use the
1611 // DiagScaledMatrixProduct transformer
1612 if(not prodTrans->isCompatible(*implicitOp))
1613 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1614
1615 // build operator and multiply
1616 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1617 prodTrans->transform(*implicitOp,explicitOp.ptr());
1618 explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1619 " * " + opr->getObjectLabel() + " )");
1620
1621 return explicitOp;
1622#else
1623 throw std::logic_error("explicitMultiply is trying to use Epetra "
1624 "code, but TEKO_HAVE_EPETRA is disabled!");
1625#endif
1626 }
1627}
1628
1642const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
1643 const ModifiableLinearOp & destOp)
1644{
1645 // if this is a blocked operator, multiply block by block
1646 // it is possible that not every factor in the product is blocked and these situations are handled separately
1647
1648 bool isBlockedL = isPhysicallyBlockedLinearOp(opl);
1649 bool isBlockedR = isPhysicallyBlockedLinearOp(opr);
1650
1651 // both factors blocked
1652 if((isBlockedL && isBlockedR)){
1653 double scalarl = 0.0;
1654 bool transpl = false;
1655 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1656 double scalarr = 0.0;
1657 bool transpr = false;
1658 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1659 double scalar = scalarl*scalarr;
1660
1661 int numRows = blocked_opl->productRange()->numBlocks();
1662 int numCols = blocked_opr->productDomain()->numBlocks();
1663 int numMiddle = blocked_opl->productDomain()->numBlocks();
1664
1665 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1666
1667 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1668 blocked_product->beginBlockFill(numRows,numCols);
1669 for(int r = 0; r < numRows; ++r){
1670 for(int c = 0; c < numCols; ++c){
1671
1672 LinearOp product_rc = explicitMultiply(blocked_opl->getBlock(r,0),blocked_opr->getBlock(0,c));
1673 for(int m = 1; m < numMiddle; ++m){
1674 LinearOp product_m = explicitMultiply(blocked_opl->getBlock(r,m),blocked_opr->getBlock(m,c));
1675 product_rc = explicitAdd(product_rc,product_m);
1676 }
1677 blocked_product->setBlock(r,c,Thyra::scale(scalar,product_rc));
1678 }
1679 }
1680 blocked_product->endBlockFill();
1681 return blocked_product;
1682 }
1683
1684 // only left factor blocked
1685 if((isBlockedL && !isBlockedR)){
1686 double scalarl = 0.0;
1687 bool transpl = false;
1688 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl,&scalarl,&transpl);
1689 double scalar = scalarl;
1690
1691 int numRows = blocked_opl->productRange()->numBlocks();
1692 int numCols = 1;
1693 int numMiddle = 1;
1694
1695 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == numMiddle);
1696
1697 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1698 blocked_product->beginBlockFill(numRows,numCols);
1699 for(int r = 0; r < numRows; ++r){
1700 LinearOp product_r = explicitMultiply(blocked_opl->getBlock(r,0),opr);
1701 blocked_product->setBlock(r,0,Thyra::scale(scalar,product_r));
1702 }
1703 blocked_product->endBlockFill();
1704 return blocked_product;
1705 }
1706
1707 // only right factor blocked
1708 if((!isBlockedL && isBlockedR)){
1709 double scalarr = 0.0;
1710 bool transpr = false;
1711 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr,&scalarr,&transpr);
1712 double scalar = scalarr;
1713
1714 int numRows = 1;
1715 int numCols = blocked_opr->productDomain()->numBlocks();
1716 int numMiddle = 1;
1717
1718 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numMiddle);
1719
1720 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_product = Thyra::defaultBlockedLinearOp<double>();
1721 blocked_product->beginBlockFill(numRows,numCols);
1722 for(int c = 0; c < numCols; ++c){
1723 LinearOp product_c = explicitMultiply(opl,blocked_opr->getBlock(0,c));
1724 blocked_product->setBlock(0,c,Thyra::scale(scalar,product_c));
1725 }
1726 blocked_product->endBlockFill();
1727 return blocked_product;
1728 }
1729
1730 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1731 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1732
1733 if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix multiply
1734
1735 // Get left and right Tpetra crs operators
1736 ST scalarl = 0.0;
1737 bool transpl = false;
1738 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1739 ST scalarr = 0.0;
1740 bool transpr = false;
1741 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1742
1743 // Build output operator
1744 RCP<Thyra::LinearOpBase<ST> > explicitOp;
1745 if(destOp!=Teuchos::null)
1746 explicitOp = destOp;
1747 else
1748 explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1749 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1750
1751 // Do explicit matrix-matrix multiply
1752 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1753 Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1754 explicitCrsOp->resumeFill();
1755 explicitCrsOp->scale(scalarl*scalarr);
1756 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1757 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1758 return tExplicitOp;
1759
1760 } else if (isTpetral && !isTpetrar){ // Assume that the right operator is diagonal
1761
1762 // Get left Tpetra crs operator
1763 ST scalarl = 0.0;
1764 bool transpl = false;
1765 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1766
1767 // Cast right operator as DiagonalLinearOp and extract diagonal as Vector
1768 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr);
1769 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),true);
1770 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1771
1772 // Scale by the diagonal operator
1773 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpl, Tpetra::Import<LO,GO,NT>(tCrsOpl->getRowMap(),tCrsOpl->getRowMap()));
1774 explicitCrsOp->rightScale(*diagPtr);
1775 explicitCrsOp->resumeFill();
1776 explicitCrsOp->scale(scalarl);
1777 explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1778 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1779
1780 } else if (!isTpetral && isTpetrar){ // Assume that the left operator is diagonal
1781
1782 // Get right Tpetra crs operator
1783 ST scalarr = 0.0;
1784 bool transpr = false;
1785 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1786
1787 // Cast leftt operator as DiagonalLinearOp and extract diagonal as Vector
1788 RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl,true);
1789 RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),true);
1790 RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1791
1792 // Scale by the diagonal operator
1793 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::importAndFillCompleteCrsMatrix<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tCrsOpr, Tpetra::Import<LO,GO,NT>(tCrsOpr->getRowMap(),tCrsOpr->getRowMap()));
1794 explicitCrsOp->leftScale(*diagPtr);
1795 explicitCrsOp->resumeFill();
1796 explicitCrsOp->scale(scalarr);
1797 explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1798 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1799
1800 } else { // Assume Epetra and we can use transformers
1801#ifdef TEKO_HAVE_EPETRA
1802 // build implicit multiply
1803 const LinearOp implicitOp = Thyra::multiply(opl,opr);
1804
1805 // build a scaling transformer
1806
1807 RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1808 = Thyra::epetraExtDiagScalingTransformer();
1809
1810 // check to see if a scaling transformer works: if not use the
1811 // DiagScaledMatrixProduct transformer
1812 if(not prodTrans->isCompatible(*implicitOp))
1813 prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1814
1815 // build operator destination operator
1816 ModifiableLinearOp explicitOp;
1817
1818 // if neccessary build a operator to put the explicit multiply into
1819 if(destOp==Teuchos::null)
1820 explicitOp = prodTrans->createOutputOp();
1821 else
1822 explicitOp = destOp;
1823
1824 // perform multiplication
1825 prodTrans->transform(*implicitOp,explicitOp.ptr());
1826
1827 // label it
1828 explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1829 " * " + opr->getObjectLabel() + " )");
1830
1831 return explicitOp;
1832#else
1833 throw std::logic_error("explicitMultiply is trying to use Epetra "
1834 "code, but TEKO_HAVE_EPETRA is disabled!");
1835#endif
1836 }
1837}
1838
1849const LinearOp explicitAdd(const LinearOp & opl_in,const LinearOp & opr_in)
1850{
1851 // if both blocked, add block by block
1852 if(isPhysicallyBlockedLinearOp(opl_in) && isPhysicallyBlockedLinearOp(opr_in)){
1853 double scalarl = 0.0;
1854 bool transpl = false;
1855 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl_in, &scalarl, &transpl);
1856
1857 double scalarr = 0.0;
1858 bool transpr = false;
1859 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr_in, &scalarr, &transpr);
1860
1861 int numRows = blocked_opl->productRange()->numBlocks();
1862 int numCols = blocked_opl->productDomain()->numBlocks();
1863 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numRows);
1864 TEUCHOS_ASSERT(blocked_opr->productDomain()->numBlocks() == numCols);
1865
1866 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_sum = Thyra::defaultBlockedLinearOp<double>();
1867 blocked_sum->beginBlockFill(numRows,numCols);
1868 for(int r = 0; r < numRows; ++r)
1869 for(int c = 0; c < numCols; ++c)
1870 blocked_sum->setBlock(r,c,explicitAdd(Thyra::scale(scalarl,blocked_opl->getBlock(r,c)),Thyra::scale(scalarr,blocked_opr->getBlock(r,c))));
1871 blocked_sum->endBlockFill();
1872 return blocked_sum;
1873 }
1874
1875 // if only one is blocked, it must be 1x1
1876 LinearOp opl = opl_in;
1877 LinearOp opr = opr_in;
1878 if(isPhysicallyBlockedLinearOp(opl_in)){
1879 double scalarl = 0.0;
1880 bool transpl = false;
1881 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl_in, &scalarl, &transpl);
1882 TEUCHOS_ASSERT(blocked_opl->productRange()->numBlocks() == 1);
1883 TEUCHOS_ASSERT(blocked_opl->productDomain()->numBlocks() == 1);
1884 opl = Thyra::scale(scalarl,blocked_opl->getBlock(0,0));
1885 }
1886 if(isPhysicallyBlockedLinearOp(opr_in)){
1887 double scalarr = 0.0;
1888 bool transpr = false;
1889 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr_in, &scalarr, &transpr);
1890 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == 1);
1891 TEUCHOS_ASSERT(blocked_opr->productDomain()->numBlocks() == 1);
1892 opr = Thyra::scale(scalarr,blocked_opr->getBlock(0,0));
1893 }
1894
1895 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1896 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1897
1898 // if one of the operators in the sum is a thyra zero op
1899 if(isZeroOp(opl)){
1900 if(isZeroOp(opr))
1901 return opr; // return a zero op if both are zero
1902 if(isTpetrar){ // if other op is tpetra, replace this with a zero crs matrix
1903 ST scalar = 0.0;
1904 bool transp = false;
1905 auto crs_op = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalar, &transp);
1906 auto zero_crs = Tpetra::createCrsMatrix<ST,LO,GO,NT>(crs_op->getRowMap());
1907 zero_crs->fillComplete();
1908 opl = Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getDomainMap()),zero_crs);
1909 isTpetral = true;
1910 } else
1911 return opr->clone();
1912 }
1913 if(isZeroOp(opr)){
1914 if(isTpetral){ // if other op is tpetra, replace this with a zero crs matrix
1915 ST scalar = 0.0;
1916 bool transp = false;
1917 auto crs_op = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalar, &transp);
1918 auto zero_crs = Tpetra::createCrsMatrix<ST,LO,GO,NT>(crs_op->getRowMap());
1919 zero_crs->fillComplete();
1920 opr = Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(crs_op->getDomainMap()),zero_crs);
1921 isTpetrar = true;
1922 } else
1923 return opl->clone();
1924 }
1925
1926 if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix add
1927
1928 // Get left and right Tpetra crs operators
1929 ST scalarl = 0.0;
1930 bool transpl = false;
1931 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1932 ST scalarr = 0.0;
1933 bool transpr = false;
1934 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1935
1936 // Build output operator
1937 RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1938 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1939
1940 // Do explicit matrix-matrix add
1941 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1942 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1943 return tExplicitOp;
1944
1945 }else{//Assume Epetra
1946#ifdef TEKO_HAVE_EPETRA
1947 // build implicit add
1948 const LinearOp implicitOp = Thyra::add(opl,opr);
1949
1950 // build transformer
1951 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1952 Thyra::epetraExtAddTransformer();
1953
1954 // build operator and add
1955 const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1956 prodTrans->transform(*implicitOp,explicitOp.ptr());
1957 explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1958 " + " + opr->getObjectLabel() + " )");
1959
1960 return explicitOp;
1961#else
1962 throw std::logic_error("explicitAdd is trying to use Epetra "
1963 "code, but TEKO_HAVE_EPETRA is disabled!");
1964#endif
1965 }
1966}
1967
1980const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
1981 const ModifiableLinearOp & destOp)
1982{
1983 // if blocked, add block by block
1984 if(isPhysicallyBlockedLinearOp(opl)){
1985 TEUCHOS_ASSERT(isPhysicallyBlockedLinearOp(opr));
1986
1987 double scalarl = 0.0;
1988 bool transpl = false;
1989 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opl = getPhysicallyBlockedLinearOp(opl, &scalarl, &transpl);
1990
1991 double scalarr = 0.0;
1992 bool transpr = false;
1993 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_opr = getPhysicallyBlockedLinearOp(opr, &scalarr, &transpr);
1994
1995 int numRows = blocked_opl->productRange()->numBlocks();
1996 int numCols = blocked_opl->productDomain()->numBlocks();
1997 TEUCHOS_ASSERT(blocked_opr->productRange()->numBlocks() == numRows);
1998 TEUCHOS_ASSERT(blocked_opr->productDomain()->numBlocks() == numCols);
1999
2000 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_sum = Teuchos::rcp_dynamic_cast<Thyra::DefaultBlockedLinearOp<double>>(destOp);
2001 if(blocked_sum.is_null()) {
2002 // take care of the null case, this means we must alllocate memory
2003 blocked_sum = Thyra::defaultBlockedLinearOp<double>();
2004
2005 blocked_sum->beginBlockFill(numRows,numCols);
2006 for(int r = 0; r < numRows; ++r) {
2007 for(int c = 0; c < numCols; ++c) {
2008 auto block = explicitAdd(Thyra::scale(scalarl,blocked_opl->getBlock(r,c)),Thyra::scale(scalarr,blocked_opr->getBlock(r,c)),Teuchos::null);
2009 blocked_sum->setNonconstBlock(r,c,block);
2010 }
2011 }
2012 blocked_sum->endBlockFill();
2013
2014 }
2015 else {
2016 // in this case memory can be reused
2017 for(int r = 0; r < numRows; ++r)
2018 for(int c = 0; c < numCols; ++c)
2019 explicitAdd(Thyra::scale(scalarl,blocked_opl->getBlock(r,c)),Thyra::scale(scalarr,blocked_opr->getBlock(r,c)),blocked_sum->getNonconstBlock(r,c));
2020 }
2021
2022 return blocked_sum;
2023 }
2024
2025 bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
2026 bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
2027
2028 if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix add
2029
2030 // Get left and right Tpetra crs operators
2031 ST scalarl = 0.0;
2032 bool transpl = false;
2033 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
2034 ST scalarr = 0.0;
2035 bool transpr = false;
2036 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
2037
2038 // Build output operator
2039 RCP<Thyra::LinearOpBase<ST> > explicitOp;
2040 if(destOp!=Teuchos::null)
2041 explicitOp = destOp;
2042 else
2043 explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
2044 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
2045
2046 // Do explicit matrix-matrix add
2047 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
2048 tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
2049 return tExplicitOp;
2050
2051 }else{ // Assume Epetra
2052#ifdef TEKO_HAVE_EPETRA
2053 // build implicit add
2054 const LinearOp implicitOp = Thyra::add(opl,opr);
2055
2056 // build transformer
2057 const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
2058 Thyra::epetraExtAddTransformer();
2059
2060 // build or reuse destination operator
2061 RCP<Thyra::LinearOpBase<double> > explicitOp;
2062 if(destOp!=Teuchos::null)
2063 explicitOp = destOp;
2064 else
2065 explicitOp = prodTrans->createOutputOp();
2066
2067 // perform add
2068 prodTrans->transform(*implicitOp,explicitOp.ptr());
2069 explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
2070 " + " + opr->getObjectLabel() + " )");
2071
2072 return explicitOp;
2073#else
2074 throw std::logic_error("explicitAdd is trying to use Epetra "
2075 "code, but TEKO_HAVE_EPETRA is disabled!");
2076#endif
2077 }
2078}
2079
2084const ModifiableLinearOp explicitSum(const LinearOp & op,
2085 const ModifiableLinearOp & destOp)
2086{
2087#ifdef TEKO_HAVE_EPETRA
2088 // convert operators to Epetra_CrsMatrix
2089 const RCP<const Epetra_CrsMatrix> epetraOp =
2090 rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*op), true);
2091
2092 if(destOp==Teuchos::null) {
2093 Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(new Epetra_CrsMatrix(*epetraOp));
2094
2095 return Thyra::nonconstEpetraLinearOp(epetraDest);
2096 }
2097
2098 const RCP<Epetra_CrsMatrix> epetraDest =
2099 rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp), true);
2100
2101 EpetraExt::MatrixMatrix::Add(*epetraOp,false,1.0,*epetraDest,1.0);
2102
2103 return destOp;
2104#else
2105 throw std::logic_error("explicitSum is trying to use Epetra "
2106 "code, but TEKO_HAVE_EPETRA is disabled!");
2107#endif
2108}
2109
2110const LinearOp explicitTranspose(const LinearOp & op)
2111{
2112 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2113
2114 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op,true);
2115 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
2116
2117 Tpetra::RowMatrixTransposer<ST,LO,GO,NT> transposer(tCrsOp);
2118 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > transOp = transposer.createTranspose();
2119
2120 return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getRangeMap()),
2121 Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getDomainMap()),
2122 transOp);
2123
2124 } else {
2125#ifdef TEKO_HAVE_EPETRA
2126 RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
2127 TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
2128 "Teko::explicitTranspose Not an Epetra_Operator");
2129 RCP<const Epetra_RowMatrix> eRowMatrixOp
2130 = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(eOp);
2131 TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
2132 "Teko::explicitTranspose Not an Epetra_RowMatrix");
2133
2134 // we now have a delete transpose operator
2135 EpetraExt::RowMatrix_Transpose tranposeOp;
2136 Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
2137
2138 // this copy is because of a poor implementation of the EpetraExt::Transform
2139 // implementation
2140 Teuchos::RCP<Epetra_CrsMatrix> crsMat
2141 = Teuchos::rcp(new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
2142
2143 return Thyra::epetraLinearOp(crsMat);
2144#else
2145 throw std::logic_error("explicitTranspose is trying to use Epetra "
2146 "code, but TEKO_HAVE_EPETRA is disabled!");
2147#endif
2148 }
2149}
2150
2151double frobeniusNorm(const LinearOp & op_in)
2152{
2153 LinearOp op;
2154 double scalar = 1.0;
2155
2156 // if blocked, must be 1x1
2157 if(isPhysicallyBlockedLinearOp(op_in)){
2158 bool transp = false;
2159 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > blocked_op = getPhysicallyBlockedLinearOp(op_in,&scalar,&transp);
2160 TEUCHOS_ASSERT(blocked_op->productRange()->numBlocks() == 1);
2161 TEUCHOS_ASSERT(blocked_op->productDomain()->numBlocks() == 1);
2162 op = blocked_op->getBlock(0,0);
2163 } else
2164 op = op_in;
2165
2166 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2167 const RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
2168 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
2169 return crsOp->getFrobeniusNorm();
2170 } else {
2171#ifdef TEKO_HAVE_EPETRA
2172 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
2173 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,true);
2174 return crsOp->NormFrobenius();
2175#else
2176 throw std::logic_error("frobeniusNorm is trying to use Epetra "
2177 "code, but TEKO_HAVE_EPETRA is disabled!");
2178#endif
2179 }
2180}
2181
2182double oneNorm(const LinearOp & op)
2183{
2184 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2185 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"One norm not currently implemented for Tpetra matrices");
2186
2187 } else {
2188#ifdef TEKO_HAVE_EPETRA
2189 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
2190 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,true);
2191 return crsOp->NormOne();
2192#else
2193 throw std::logic_error("oneNorm is trying to use Epetra "
2194 "code, but TEKO_HAVE_EPETRA is disabled!");
2195#endif
2196 }
2197}
2198
2199double infNorm(const LinearOp & op)
2200{
2201 if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
2202 ST scalar = 0.0;
2203 bool transp = false;
2204 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
2205
2206 // extract diagonal
2207 const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
2208 Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
2209
2210 // compute absolute value row sum
2211 diag.putScalar(0.0);
2212 for(LO i=0;i<(LO) tCrsOp->getLocalNumRows();i++) {
2213 LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
2214 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::local_inds_host_view_type indices;
2215 typename Tpetra::CrsMatrix<ST,LO,GO,NT>::values_host_view_type values;
2216 tCrsOp->getLocalRowView(i,indices,values);
2217
2218 // build abs value row sum
2219 for(LO j=0;j<numEntries;j++)
2220 diag.sumIntoLocalValue(i,std::abs(values(j)));
2221 }
2222 return diag.normInf()*scalar;
2223
2224 } else {
2225#ifdef TEKO_HAVE_EPETRA
2226 const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
2227 const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,true);
2228 return crsOp->NormInf();
2229#else
2230 throw std::logic_error("infNorm is trying to use Epetra "
2231 "code, but TEKO_HAVE_EPETRA is disabled!");
2232#endif
2233 }
2234}
2235
2236const LinearOp buildDiagonal(const MultiVector & src,const std::string & lbl)
2237{
2238 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range());
2239 Thyra::copy(*src->col(0),dst.ptr());
2240
2241 return Thyra::diagonal<double>(dst,lbl);
2242}
2243
2244const LinearOp buildInvDiagonal(const MultiVector & src,const std::string & lbl)
2245{
2246 const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
2247 RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range());
2248 Thyra::reciprocal<double>(*srcV,dst.ptr());
2249
2250 return Thyra::diagonal<double>(dst,lbl);
2251}
2252
2254BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvv)
2255{
2256 Teuchos::Array<MultiVector> mvA;
2257 Teuchos::Array<VectorSpace> vsA;
2258
2259 // build arrays of multi vectors and vector spaces
2260 std::vector<MultiVector>::const_iterator itr;
2261 for(itr=mvv.begin();itr!=mvv.end();++itr) {
2262 mvA.push_back(*itr);
2263 vsA.push_back((*itr)->range());
2264 }
2265
2266 // construct the product vector space
2267 const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
2268 = Thyra::productVectorSpace<double>(vsA);
2269
2270 return Thyra::defaultProductMultiVector<double>(vs,mvA);
2271}
2272
2283Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
2284 const std::vector<int> & indices,
2285 const VectorSpace & vs,
2286 double onValue, double offValue)
2287
2288{
2289 using Teuchos::RCP;
2290
2291 // create a new vector
2292 RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
2293 Thyra::put_scalar<double>(offValue,v.ptr()); // fill it with "off" values
2294
2295 // set on values
2296 for(std::size_t i=0;i<indices.size();i++)
2297 Thyra::set_ele<double>(indices[i],onValue,v.ptr());
2298
2299 return v;
2300}
2301
2326double computeSpectralRad(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
2327 bool isHermitian,int numBlocks,int restart,int verbosity)
2328{
2329 typedef Thyra::LinearOpBase<double> OP;
2330 typedef Thyra::MultiVectorBase<double> MV;
2331
2332 int startVectors = 1;
2333
2334 // construct an initial guess
2335 const RCP<MV> ivec = Thyra::createMember(A->domain());
2336 Thyra::randomize(-1.0,1.0,ivec.ptr());
2337
2338 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
2339 = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
2340 eigProb->setNEV(1);
2341 eigProb->setHermitian(isHermitian);
2342
2343 // set the problem up
2344 if(not eigProb->setProblem()) {
2345 // big time failure!
2346 return Teuchos::ScalarTraits<double>::nan();
2347 }
2348
2349 // we want largert magnitude eigenvalue
2350 std::string which("LM"); // largest magnitude
2351
2352 // Create the parameter list for the eigensolver
2353 // verbosity+=Anasazi::TimingDetails;
2354 Teuchos::ParameterList MyPL;
2355 MyPL.set( "Verbosity", verbosity );
2356 MyPL.set( "Which", which );
2357 MyPL.set( "Block Size", startVectors );
2358 MyPL.set( "Num Blocks", numBlocks );
2359 MyPL.set( "Maximum Restarts", restart );
2360 MyPL.set( "Convergence Tolerance", tol );
2361
2362 // build status test manager
2363 // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
2364 // = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(numBlocks*(restart+1)));
2365
2366 // Create the Block Krylov Schur solver
2367 // This takes as inputs the eigenvalue problem and the solver parameters
2368 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
2369
2370 // Solve the eigenvalue problem, and save the return code
2371 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
2372
2373 if(solverreturn==Anasazi::Unconverged) {
2374 double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
2375 double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
2376
2377 return -std::abs(std::sqrt(real*real+comp*comp));
2378
2379 // cout << "Anasazi::BlockKrylovSchur::solve() did not converge!" << std::endl;
2380 // return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
2381 }
2382 else { // solverreturn==Anasazi::Converged
2383 double real = eigProb->getSolution().Evals.begin()->realpart;
2384 double comp = eigProb->getSolution().Evals.begin()->imagpart;
2385
2386 return std::abs(std::sqrt(real*real+comp*comp));
2387
2388 // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
2389 // return std::abs(eigProb->getSolution().Evals.begin()->realpart);
2390 }
2391}
2392
2416double computeSmallestMagEig(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
2417 bool isHermitian,int numBlocks,int restart,int verbosity)
2418{
2419 typedef Thyra::LinearOpBase<double> OP;
2420 typedef Thyra::MultiVectorBase<double> MV;
2421
2422 int startVectors = 1;
2423
2424 // construct an initial guess
2425 const RCP<MV> ivec = Thyra::createMember(A->domain());
2426 Thyra::randomize(-1.0,1.0,ivec.ptr());
2427
2428 RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
2429 = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
2430 eigProb->setNEV(1);
2431 eigProb->setHermitian(isHermitian);
2432
2433 // set the problem up
2434 if(not eigProb->setProblem()) {
2435 // big time failure!
2436 return Teuchos::ScalarTraits<double>::nan();
2437 }
2438
2439 // we want largert magnitude eigenvalue
2440 std::string which("SM"); // smallest magnitude
2441
2442 // Create the parameter list for the eigensolver
2443 Teuchos::ParameterList MyPL;
2444 MyPL.set( "Verbosity", verbosity );
2445 MyPL.set( "Which", which );
2446 MyPL.set( "Block Size", startVectors );
2447 MyPL.set( "Num Blocks", numBlocks );
2448 MyPL.set( "Maximum Restarts", restart );
2449 MyPL.set( "Convergence Tolerance", tol );
2450
2451 // build status test manager
2452 // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
2453 // = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(10));
2454
2455 // Create the Block Krylov Schur solver
2456 // This takes as inputs the eigenvalue problem and the solver parameters
2457 Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
2458
2459 // Solve the eigenvalue problem, and save the return code
2460 Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
2461
2462 if(solverreturn==Anasazi::Unconverged) {
2463 // cout << "Anasazi::BlockKrylovSchur::solve() did not converge! " << std::endl;
2464 return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
2465 }
2466 else { // solverreturn==Anasazi::Converged
2467 // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
2468 return std::abs(eigProb->getSolution().Evals.begin()->realpart);
2469 }
2470}
2471
2480ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt)
2481{
2482 switch(dt) {
2483 case Diagonal:
2484 return getDiagonalOp(A);
2485 case Lumped:
2486 return getLumpedMatrix(A);
2487 case AbsRowSum:
2488 return getAbsRowSumMatrix(A);
2489 case NotDiag:
2490 default:
2491 TEUCHOS_TEST_FOR_EXCEPT(true);
2492 };
2493
2494 return Teuchos::null;
2495}
2496
2505ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const Teko::DiagonalType & dt)
2506{
2507 switch(dt) {
2508 case Diagonal:
2509 return getInvDiagonalOp(A);
2510 case Lumped:
2511 return getInvLumpedMatrix(A);
2512 case AbsRowSum:
2513 return getAbsRowSumInvMatrix(A);
2514 case NotDiag:
2515 default:
2516 TEUCHOS_TEST_FOR_EXCEPT(true);
2517 };
2518
2519 return Teuchos::null;
2520}
2521
2528std::string getDiagonalName(const DiagonalType & dt)
2529{
2530 switch(dt) {
2531 case Diagonal:
2532 return "Diagonal";
2533 case Lumped:
2534 return "Lumped";
2535 case AbsRowSum:
2536 return "AbsRowSum";
2537 case NotDiag:
2538 return "NotDiag";
2539 case BlkDiag:
2540 return "BlkDiag";
2541 };
2542
2543 return "<error>";
2544}
2545
2554DiagonalType getDiagonalType(std::string name)
2555{
2556 if(name=="Diagonal")
2557 return Diagonal;
2558 if(name=="Lumped")
2559 return Lumped;
2560 if(name=="AbsRowSum")
2561 return AbsRowSum;
2562 if(name=="BlkDiag")
2563 return BlkDiag;
2564
2565 return NotDiag;
2566}
2567
2568#ifdef TEKO_HAVE_EPETRA
2569LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,const LinearOp & Op){
2570#ifdef Teko_ENABLE_Isorropia
2571 Teuchos::ParameterList probeList;
2572 Prober prober(G,probeList,true);
2573 Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(new Epetra_CrsMatrix(Copy,*G));
2575 prober.probe(Mwrap,*Mat);
2576 return Thyra::epetraLinearOp(Mat);
2577#else
2578 (void)G;
2579 (void)Op;
2580 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,"Probe requires Isorropia");
2581#endif
2582}
2583#endif
2584
2585double norm_1(const MultiVector & v,std::size_t col)
2586{
2587 Teuchos::Array<double> n(v->domain()->dim());
2588 Thyra::norms_1<double>(*v,n);
2589
2590 return n[col];
2591}
2592
2593double norm_2(const MultiVector & v,std::size_t col)
2594{
2595 Teuchos::Array<double> n(v->domain()->dim());
2596 Thyra::norms_2<double>(*v,n);
2597
2598 return n[col];
2599}
2600
2601#ifdef TEKO_HAVE_EPETRA
2602void putScalar(const ModifiableLinearOp & op,double scalar)
2603{
2604 try {
2605 // get Epetra_Operator
2606 RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
2607
2608 // cast it to a CrsMatrix
2609 RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
2610
2611 eCrsOp->PutScalar(scalar);
2612 }
2613 catch (std::exception & e) {
2614 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2615
2616 *out << "Teko: putScalar requires an Epetra_CrsMatrix\n";
2617 *out << " Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
2618 *out << " OR\n";
2619 *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
2620 *out << std::endl;
2621 *out << "*** THROWN EXCEPTION ***\n";
2622 *out << e.what() << std::endl;
2623 *out << "************************\n";
2624
2625 throw e;
2626 }
2627}
2628#endif
2629
2630void clipLower(MultiVector & v,double lowerBound)
2631{
2632 using Teuchos::RCP;
2633 using Teuchos::rcp_dynamic_cast;
2634
2635 // cast so entries are accessible
2636 // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2637 // = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
2638
2639 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2640 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2641 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2642
2643 Teuchos::ArrayRCP<double> values;
2644 // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2645 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2646 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2647 if(values[j]<lowerBound)
2648 values[j] = lowerBound;
2649 }
2650 }
2651}
2652
2653void clipUpper(MultiVector & v,double upperBound)
2654{
2655 using Teuchos::RCP;
2656 using Teuchos::rcp_dynamic_cast;
2657
2658 // cast so entries are accessible
2659 // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2660 // = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
2661 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2662 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2663 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2664
2665 Teuchos::ArrayRCP<double> values;
2666 // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2667 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2668 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2669 if(values[j]>upperBound)
2670 values[j] = upperBound;
2671 }
2672 }
2673}
2674
2675void replaceValue(MultiVector & v,double currentValue,double newValue)
2676{
2677 using Teuchos::RCP;
2678 using Teuchos::rcp_dynamic_cast;
2679
2680 // cast so entries are accessible
2681 // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2682 // = rcp_dynamic_cast<Thyra::SpmdMultiVectorBase<double> >(v,true);
2683 for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2684 RCP<Thyra::SpmdVectorBase<double> > spmdVec
2685 = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2686
2687 Teuchos::ArrayRCP<double> values;
2688 // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2689 spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2690 for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2691 if(values[j]==currentValue)
2692 values[j] = newValue;
2693 }
2694 }
2695}
2696
2697void columnAverages(const MultiVector & v,std::vector<double> & averages)
2698{
2699 averages.resize(v->domain()->dim());
2700
2701 // sum over each column
2702 Thyra::sums<double>(*v,averages);
2703
2704 // build averages
2705 Thyra::Ordinal rows = v->range()->dim();
2706 for(std::size_t i=0;i<averages.size();i++)
2707 averages[i] = averages[i]/rows;
2708}
2709
2710double average(const MultiVector & v)
2711{
2712 Thyra::Ordinal rows = v->range()->dim();
2713 Thyra::Ordinal cols = v->domain()->dim();
2714
2715 std::vector<double> averages;
2716 columnAverages(v,averages);
2717
2718 double sum = 0.0;
2719 for(std::size_t i=0;i<averages.size();i++)
2720 sum += averages[i] * rows;
2721
2722 return sum/(rows*cols);
2723}
2724
2725bool isPhysicallyBlockedLinearOp(const LinearOp & op)
2726{
2727 // See if the operator is a PBLO
2728 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
2729 if (!pblo.is_null())
2730 return true;
2731
2732 // See if the operator is a wrapped PBLO
2733 ST scalar = 0.0;
2734 Thyra::EOpTransp transp = Thyra::NOTRANS;
2735 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
2736 Thyra::unwrap(op, &scalar, &transp, &wrapped_op);
2737 pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(wrapped_op);
2738 if (!pblo.is_null())
2739 return true;
2740
2741 return false;
2742}
2743
2744RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > getPhysicallyBlockedLinearOp(const LinearOp & op, ST *scalar, bool *transp)
2745{
2746 // If the operator is a TpetraLinearOp
2747 RCP<const Thyra::PhysicallyBlockedLinearOpBase<double> > pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(op);
2748 if(!pblo.is_null()){
2749 *scalar = 1.0;
2750 *transp = false;
2751 return pblo;
2752 }
2753
2754 // If the operator is a wrapped TpetraLinearOp
2755 RCP<const Thyra::LinearOpBase<ST> > wrapped_op;
2756 Thyra::EOpTransp eTransp = Thyra::NOTRANS;
2757 Thyra::unwrap(op, scalar, &eTransp, &wrapped_op);
2758 pblo = rcp_dynamic_cast<const Thyra::PhysicallyBlockedLinearOpBase<double> >(wrapped_op,true);
2759 if(!pblo.is_null()){
2760 *transp = true;
2761 if(eTransp == Thyra::NOTRANS)
2762 *transp = false;
2763 return pblo;
2764 }
2765
2766 return Teuchos::null;
2767}
2768
2769
2770}
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
DiagonalType
Type describing the type of diagonal to construct.
@ BlkDiag
Specifies that a block diagonal approximation is to be used.
@ NotDiag
For user convenience, if Teko recieves this value, exceptions will be thrown.
@ AbsRowSum
Specifies that the diagonal entry is .
@ Diagonal
Specifies that just the diagonal is used.
@ Lumped
Specifies that row sum is used to form a diagonal.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...