Teko  Version of the Day
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_EpetraExtDiagScaledMatProdTransformer.hpp"
56 #include "Thyra_EpetraExtDiagScalingTransformer.hpp"
57 #include "Thyra_EpetraExtAddTransformer.hpp"
58 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
59 #include "Thyra_DefaultMultipliedLinearOp.hpp"
60 #include "Thyra_DefaultZeroLinearOp.hpp"
61 #include "Thyra_DefaultProductMultiVector.hpp"
62 #include "Thyra_DefaultProductVectorSpace.hpp"
63 #include "Thyra_MultiVectorStdOps.hpp"
64 #include "Thyra_VectorStdOps.hpp"
65 #include "Thyra_SpmdVectorBase.hpp"
66 #include "Thyra_get_Epetra_Operator.hpp"
67 #include "Thyra_EpetraThyraWrappers.hpp"
68 #include "Thyra_EpetraOperatorWrapper.hpp"
69 #include "Thyra_EpetraLinearOp.hpp"
70 
71 // Teuchos includes
72 #include "Teuchos_Array.hpp"
73 
74 // Epetra includes
75 #include "Epetra_Operator.h"
76 #include "Epetra_CrsGraph.h"
77 #include "Epetra_CrsMatrix.h"
78 #include "Epetra_Vector.h"
79 #include "Epetra_Map.h"
80 
81 #include "EpetraExt_Transpose_RowMatrix.h"
82 #include "EpetraExt_MatrixMatrix.h"
83 
84 // Anasazi includes
85 #include "AnasaziBasicEigenproblem.hpp"
86 #include "AnasaziThyraAdapter.hpp"
87 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
88 #include "AnasaziBlockKrylovSchur.hpp"
89 #include "AnasaziStatusTestMaxIters.hpp"
90 
91 // Isorropia includes
92 #ifdef Teko_ENABLE_Isorropia
93 #include "Isorropia_EpetraProber.hpp"
94 #endif
95 
96 // Teko includes
97 #include "Teko_EpetraHelpers.hpp"
98 #include "Teko_EpetraOperatorWrapper.hpp"
99 #include "Teko_TpetraHelpers.hpp"
100 #include "Teko_TpetraOperatorWrapper.hpp"
101 
102 // Tpetra
103 #include "Thyra_TpetraLinearOp.hpp"
104 #include "Tpetra_CrsMatrix.hpp"
105 #include "Tpetra_Vector.hpp"
106 #include "Thyra_TpetraThyraWrappers.hpp"
107 #include "TpetraExt_MatrixMatrix.hpp"
108 #include "Tpetra_RowMatrixTransposer.hpp"
109 
110 #include <cmath>
111 
112 namespace Teko {
113 
114 using Teuchos::rcp;
115 using Teuchos::rcp_dynamic_cast;
116 using Teuchos::RCP;
117 #ifdef Teko_ENABLE_Isorropia
118 using Isorropia::Epetra::Prober;
119 #endif
120 
121 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
122 {
123  Teuchos::RCP<Teuchos::FancyOStream> os =
124  Teuchos::VerboseObjectBase::getDefaultOStream();
125 
126  //os->setShowProcRank(true);
127  //os->setOutputToRootOnly(-1);
128  return os;
129 }
130 
131 // distance function...not parallel...entirely internal to this cpp file
132 inline double dist(int dim,double * coords,int row,int col)
133 {
134  double value = 0.0;
135  for(int i=0;i<dim;i++)
136  value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
137 
138  // the distance between the two
139  return std::sqrt(value);
140 }
141 
142 // distance function...not parallel...entirely internal to this cpp file
143 inline double dist(double * x,double * y,double * z,int stride,int row,int col)
144 {
145  double value = 0.0;
146  if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
147  if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
148  if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
149 
150  // the distance between the two
151  return std::sqrt(value);
152 }
153 
172 RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil)
173 {
174  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
175  RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
176  stencil.MaxNumEntries()+1),true);
177 
178  // allocate an additional value for the diagonal, if neccessary
179  std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
180  std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
181 
182  // loop over all the rows
183  for(int j=0;j<gl->NumMyRows();j++) {
184  int row = gl->GRID(j);
185  double diagValue = 0.0;
186  int diagInd = -1;
187  int rowSz = 0;
188 
189  // extract a copy of this row...put it in rowData, rowIndicies
190  stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
191 
192  // loop over elements of row
193  for(int i=0;i<rowSz;i++) {
194  int col = rowInd[i];
195 
196  // is this a 0 entry masquerading as some thing else?
197  double value = rowData[i];
198  if(value==0) continue;
199 
200  // for nondiagonal entries
201  if(row!=col) {
202  double d = dist(dim,coords,row,col);
203  rowData[i] = -1.0/d;
204  diagValue += rowData[i];
205  }
206  else
207  diagInd = i;
208  }
209 
210  // handle diagonal entry
211  if(diagInd<0) { // diagonal not in row
212  rowData[rowSz] = -diagValue;
213  rowInd[rowSz] = row;
214  rowSz++;
215  }
216  else { // diagonal in row
217  rowData[diagInd] = -diagValue;
218  rowInd[diagInd] = row;
219  }
220 
221  // insert row data into graph Laplacian matrix
222  TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
223  }
224 
225  gl->FillComplete();
226 
227  return gl;
228 }
229 
230 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(int dim,ST * coords,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
231 {
232  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
233  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
234  stencil.getGlobalMaxNumRowEntries()+1));
235 
236  // allocate an additional value for the diagonal, if neccessary
237  std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
238  std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
239 
240  // loop over all the rows
241  for(LO j=0;j< (LO) gl->getNodeNumRows();j++) {
242  GO row = gl->getRowMap()->getGlobalElement(j);
243  ST diagValue = 0.0;
244  GO diagInd = -1;
245  size_t rowSz = 0;
246 
247  // extract a copy of this row...put it in rowData, rowIndicies
248  stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
249 
250  // loop over elements of row
251  for(size_t i=0;i<rowSz;i++) {
252  GO col = rowInd[i];
253 
254  // is this a 0 entry masquerading as some thing else?
255  ST value = rowData[i];
256  if(value==0) continue;
257 
258  // for nondiagonal entries
259  if(row!=col) {
260  ST d = dist(dim,coords,row,col);
261  rowData[i] = -1.0/d;
262  diagValue += rowData[i];
263  }
264  else
265  diagInd = i;
266  }
267 
268  // handle diagonal entry
269  if(diagInd<0) { // diagonal not in row
270  rowData[rowSz] = -diagValue;
271  rowInd[rowSz] = row;
272  rowSz++;
273  }
274  else { // diagonal in row
275  rowData[diagInd] = -diagValue;
276  rowInd[diagInd] = row;
277  }
278 
279  // insert row data into graph Laplacian matrix
280  gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
281  }
282 
283  gl->fillComplete();
284 
285  return gl;
286 }
287 
310 RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil)
311 {
312  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
313  RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
314  stencil.MaxNumEntries()+1),true);
315 
316  // allocate an additional value for the diagonal, if neccessary
317  std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
318  std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
319 
320  // loop over all the rows
321  for(int j=0;j<gl->NumMyRows();j++) {
322  int row = gl->GRID(j);
323  double diagValue = 0.0;
324  int diagInd = -1;
325  int rowSz = 0;
326 
327  // extract a copy of this row...put it in rowData, rowIndicies
328  stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
329 
330  // loop over elements of row
331  for(int i=0;i<rowSz;i++) {
332  int col = rowInd[i];
333 
334  // is this a 0 entry masquerading as some thing else?
335  double value = rowData[i];
336  if(value==0) continue;
337 
338  // for nondiagonal entries
339  if(row!=col) {
340  double d = dist(x,y,z,stride,row,col);
341  rowData[i] = -1.0/d;
342  diagValue += rowData[i];
343  }
344  else
345  diagInd = i;
346  }
347 
348  // handle diagonal entry
349  if(diagInd<0) { // diagonal not in row
350  rowData[rowSz] = -diagValue;
351  rowInd[rowSz] = row;
352  rowSz++;
353  }
354  else { // diagonal in row
355  rowData[diagInd] = -diagValue;
356  rowInd[diagInd] = row;
357  }
358 
359  // insert row data into graph Laplacian matrix
360  TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
361  }
362 
363  gl->FillComplete();
364 
365  return gl;
366 }
367 
368 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
369 {
370  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
371  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
372  stencil.getGlobalMaxNumRowEntries()+1),true);
373 
374  // allocate an additional value for the diagonal, if neccessary
375  std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
376  std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
377 
378  // loop over all the rows
379  for(LO j=0;j<(LO) gl->getNodeNumRows();j++) {
380  GO row = gl->getRowMap()->getGlobalElement(j);
381  ST diagValue = 0.0;
382  GO diagInd = -1;
383  size_t rowSz = 0;
384 
385  // extract a copy of this row...put it in rowData, rowIndicies
386  stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
387 
388  // loop over elements of row
389  for(size_t i=0;i<rowSz;i++) {
390  GO col = rowInd[i];
391 
392  // is this a 0 entry masquerading as some thing else?
393  ST value = rowData[i];
394  if(value==0) continue;
395 
396  // for nondiagonal entries
397  if(row!=col) {
398  ST d = dist(x,y,z,stride,row,col);
399  rowData[i] = -1.0/d;
400  diagValue += rowData[i];
401  }
402  else
403  diagInd = i;
404  }
405 
406  // handle diagonal entry
407  if(diagInd<0) { // diagonal not in row
408  rowData[rowSz] = -diagValue;
409  rowInd[rowSz] = row;
410  rowSz++;
411  }
412  else { // diagonal in row
413  rowData[diagInd] = -diagValue;
414  rowInd[diagInd] = row;
415  }
416 
417  // insert row data into graph Laplacian matrix
418  gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
419  }
420 
421  gl->fillComplete();
422 
423  return gl;
424 }
425 
441 void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
442 {
443  Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
444 }
445 
461 void applyTransposeOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
462 {
463  Thyra::apply(*A,Thyra::TRANS,*x,y.ptr(),alpha,beta);
464 }
465 
468 void update(double alpha,const MultiVector & x,double beta,MultiVector & y)
469 {
470  Teuchos::Array<double> scale;
471  Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > > vec;
472 
473  // build arrays needed for linear combo
474  scale.push_back(alpha);
475  vec.push_back(x.ptr());
476 
477  // compute linear combination
478  Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
479 }
480 
482 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
483 {
484  int rows = blockRowCount(blo);
485 
486  TEUCHOS_ASSERT(rows==blockColCount(blo));
487 
488  RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
489  RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
490 
491  // allocate new operator
492  BlockedLinearOp upper = createBlockedOp();
493 
494  // build new operator
495  upper->beginBlockFill(rows,rows);
496 
497  for(int i=0;i<rows;i++) {
498  // put zero operators on the diagonal
499  // this gurantees the vector space of
500  // the new operator are fully defined
501  RCP<const Thyra::LinearOpBase<double> > zed
502  = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
503  upper->setBlock(i,i,zed);
504 
505  for(int j=i+1;j<rows;j++) {
506  // get block i,j
507  LinearOp uij = blo->getBlock(i,j);
508 
509  // stuff it in U
510  if(uij!=Teuchos::null)
511  upper->setBlock(i,j,uij);
512  }
513  }
514  if(callEndBlockFill)
515  upper->endBlockFill();
516 
517  return upper;
518 }
519 
521 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
522 {
523  int rows = blockRowCount(blo);
524 
525  TEUCHOS_ASSERT(rows==blockColCount(blo));
526 
527  RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
528  RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
529 
530  // allocate new operator
531  BlockedLinearOp lower = createBlockedOp();
532 
533  // build new operator
534  lower->beginBlockFill(rows,rows);
535 
536  for(int i=0;i<rows;i++) {
537  // put zero operators on the diagonal
538  // this gurantees the vector space of
539  // the new operator are fully defined
540  RCP<const Thyra::LinearOpBase<double> > zed
541  = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
542  lower->setBlock(i,i,zed);
543 
544  for(int j=0;j<i;j++) {
545  // get block i,j
546  LinearOp uij = blo->getBlock(i,j);
547 
548  // stuff it in U
549  if(uij!=Teuchos::null)
550  lower->setBlock(i,j,uij);
551  }
552  }
553  if(callEndBlockFill)
554  lower->endBlockFill();
555 
556  return lower;
557 }
558 
578 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo)
579 {
580  int rows = blockRowCount(blo);
581 
582  TEUCHOS_ASSERT(rows==blockColCount(blo)); // assert that matrix is square
583 
584  RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
585  RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
586 
587  // allocate new operator
588  BlockedLinearOp zeroOp = createBlockedOp();
589 
590  // build new operator
591  zeroOp->beginBlockFill(rows,rows);
592 
593  for(int i=0;i<rows;i++) {
594  // put zero operators on the diagonal
595  // this gurantees the vector space of
596  // the new operator are fully defined
597  RCP<const Thyra::LinearOpBase<double> > zed
598  = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
599  zeroOp->setBlock(i,i,zed);
600  }
601 
602  return zeroOp;
603 }
604 
606 bool isZeroOp(const LinearOp op)
607 {
608  // if operator is null...then its zero!
609  if(op==Teuchos::null) return true;
610 
611  // try to cast it to a zero linear operator
612  const LinearOp test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(op);
613 
614  // if it works...then its zero...otherwise its null
615  return test!=Teuchos::null;
616 }
617 
626 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op)
627 {
628  bool isTpetra = false;
629  RCP<const Epetra_CrsMatrix> eCrsOp;
630  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
631 
632  try {
633  // get Epetra or Tpetra Operator
634  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
635  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
636 
637  // cast it to a CrsMatrix
638  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
639  if (!eOp.is_null()){
640  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
641  }
642  else if (!tOp.is_null()){
643  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
644  isTpetra = true;
645  }
646  else
647  throw std::logic_error("Neither Epetra nor Tpetra");
648  }
649  catch (std::exception & e) {
650  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
651 
652  *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
653  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
654  *out << " OR\n";
655  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
656  *out << std::endl;
657  *out << "*** THROWN EXCEPTION ***\n";
658  *out << e.what() << std::endl;
659  *out << "************************\n";
660 
661  throw e;
662  }
663 
664  if(!isTpetra){
665  // extract diagonal
666  const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
667  Epetra_Vector & diag = *ptrDiag;
668 
669  // compute absolute value row sum
670  diag.PutScalar(0.0);
671  for(int i=0;i<eCrsOp->NumMyRows();i++) {
672  double * values = 0;
673  int numEntries;
674  eCrsOp->ExtractMyRowView(i,numEntries,values);
675 
676  // build abs value row sum
677  for(int j=0;j<numEntries;j++)
678  diag[i] += std::abs(values[j]);
679  }
680 
681  // build Thyra diagonal operator
682  return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
683  }
684  else {
685  // extract diagonal
686  const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
687  Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
688 
689  // compute absolute value row sum
690  diag.putScalar(0.0);
691  for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
692  LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
693  std::vector<LO> indices(numEntries);
694  std::vector<ST> values(numEntries);
695  Teuchos::ArrayView<const LO> indices_av(indices);
696  Teuchos::ArrayView<const ST> values_av(values);
697  tCrsOp->getLocalRowView(i,indices_av,values_av);
698 
699  // build abs value row sum
700  for(LO j=0;j<numEntries;j++)
701  diag.sumIntoLocalValue(i,std::abs(values_av[j]));
702  }
703 
704  // build Thyra diagonal operator
705  return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),"absRowSum( " + op->getObjectLabel() + " ))");
706 
707  }
708 
709 }
710 
719 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op)
720 {
721  if(Teko::TpetraHelpers::isTpetraLinearOp(op)) {
722  ST scalar = 0.0;
723  bool transp = false;
724  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
725 
726  // extract diagonal
727  const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
728  Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
729 
730  // compute absolute value row sum
731  diag.putScalar(0.0);
732  for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
733  LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
734  std::vector<LO> indices(numEntries);
735  std::vector<ST> values(numEntries);
736  Teuchos::ArrayView<const LO> indices_av(indices);
737  Teuchos::ArrayView<const ST> values_av(values);
738  tCrsOp->getLocalRowView(i,indices_av,values_av);
739 
740  // build abs value row sum
741  for(LO j=0;j<numEntries;j++)
742  diag.sumIntoLocalValue(i,std::abs(values_av[j]));
743  }
744  diag.scale(scalar);
745  diag.reciprocal(diag); // invert entries
746 
747  // build Thyra diagonal operator
748  return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),"absRowSum( " + op->getObjectLabel() + " ))");
749 
750  }
751  else{
752  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op,true);
753  RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
754 
755  // extract diagonal
756  const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
757  Epetra_Vector & diag = *ptrDiag;
758 
759  // compute absolute value row sum
760  diag.PutScalar(0.0);
761  for(int i=0;i<eCrsOp->NumMyRows();i++) {
762  double * values = 0;
763  int numEntries;
764  eCrsOp->ExtractMyRowView(i,numEntries,values);
765 
766  // build abs value row sum
767  for(int j=0;j<numEntries;j++)
768  diag[i] += std::abs(values[j]);
769  }
770  diag.Reciprocal(diag); // invert entries
771 
772  // build Thyra diagonal operator
773  return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
774  }
775 
776 }
777 
785 ModifiableLinearOp getLumpedMatrix(const LinearOp & op)
786 {
787  RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
788  RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
789 
790  // set to all ones
791  Thyra::assign(ones.ptr(),1.0);
792 
793  // compute lumped diagonal
794  // Thyra::apply(*op,Thyra::NONCONJ_ELE,*ones,&*diag);
795  Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
796 
797  return rcp(new Thyra::DefaultDiagonalLinearOp<ST>(diag));
798 }
799 
808 ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op)
809 {
810  RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
811  RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
812 
813  // set to all ones
814  Thyra::assign(ones.ptr(),1.0);
815 
816  // compute lumped diagonal
817  Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
818  Thyra::reciprocal(*diag,diag.ptr());
819 
820  return rcp(new Thyra::DefaultDiagonalLinearOp<ST>(diag));
821 }
822 
834 const ModifiableLinearOp getDiagonalOp(const LinearOp & op)
835 {
836  bool isTpetra = false;
837  RCP<const Epetra_CrsMatrix> eCrsOp;
838  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
839 
840  try {
841  // get Epetra or Tpetra Operator
842  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
843  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
844 
845  // cast it to a CrsMatrix
846  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
847  if (!eOp.is_null()){
848  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
849  }
850  else if (!tOp.is_null()){
851  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
852  isTpetra = true;
853  }
854  else
855  throw std::logic_error("Neither Epetra nor Tpetra");
856  }
857  catch (std::exception & e) {
858  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
859 
860  *out << "Teko: getDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
861  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
862  *out << " OR\n";
863  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
864  *out << std::endl;
865  *out << "*** THROWN EXCEPTION ***\n";
866  *out << e.what() << std::endl;
867  *out << "************************\n";
868 
869  throw e;
870  }
871 
872  if(!isTpetra){
873  // extract diagonal
874  const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
875  TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
876 
877  // build Thyra diagonal operator
878  return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
879  }
880  else {
881  // extract diagonal
882  const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
883  tCrsOp->getLocalDiagCopy(*diag);
884 
885  // build Thyra diagonal operator
886  return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
887 
888  }
889 }
890 
891 const MultiVector getDiagonal(const LinearOp & op)
892 {
893  bool isTpetra = false;
894  RCP<const Epetra_CrsMatrix> eCrsOp;
895  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
896 
897  try {
898  // get Epetra or Tpetra Operator
899  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
900  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
901 
902  // cast it to a CrsMatrix
903  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
904  if (!eOp.is_null()){
905  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
906  }
907  else if (!tOp.is_null()){
908  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
909  isTpetra = true;
910  }
911  else
912  throw std::logic_error("Neither Epetra nor Tpetra");
913  }
914  catch (std::exception & e) {
915  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
916 
917  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
918  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
919  *out << eOp;
920  *out << tOp;
921 
922  *out << "Teko: getDiagonal requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
923  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
924  *out << " OR\n";
925  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
926  *out << std::endl;
927  *out << "*** THROWN EXCEPTION ***\n";
928  *out << e.what() << std::endl;
929  *out << "************************\n";
930 
931  throw e;
932  }
933 
934  if(!isTpetra){
935  // extract diagonal
936  const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
937  TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
938 
939  return Thyra::create_Vector(diag,Thyra::create_VectorSpace(Teuchos::rcpFromRef(eCrsOp->RowMap())));
940  }
941  else {
942  // extract diagonal
943  const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
944  tCrsOp->getLocalDiagCopy(*diag);
945 
946  // build Thyra diagonal operator
947  return Thyra::createVector<ST,LO,GO,NT>(diag,Thyra::createVectorSpace<ST,LO,GO,NT>(tCrsOp->getRowMap()));
948 
949  }
950 }
951 
952 const MultiVector getDiagonal(const Teko::LinearOp & A,const DiagonalType & dt)
953 {
954  LinearOp diagOp = Teko::getDiagonalOp(A,dt);
955 
956  Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
957  Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag();
958  return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
959 }
960 
972 const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op)
973 {
974  if (Teko::TpetraHelpers::isTpetraLinearOp(op)){
975  ST scalar = 0.0;
976  bool transp = false;
977  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = Teko::TpetraHelpers::getTpetraCrsMatrix(op, &scalar, &transp);
978 
979  // extract diagonal
980  const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
981  diag->scale(scalar);
982  tCrsOp->getLocalDiagCopy(*diag);
983  diag->reciprocal(*diag);
984 
985  // build Thyra diagonal operator
986  return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
987 
988  }
989  else {
990  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op,true);
991  RCP<const Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
992 
993  // extract diagonal
994  const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
995  TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
996  diag->Reciprocal(*diag);
997 
998  // build Thyra diagonal operator
999  return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
1000  }
1001 }
1002 
1015 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr)
1016 {
1017  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1018  bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1019  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1020 
1021  if(isTpetral && isTpetram && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1022 
1023  // Get left and right Tpetra crs operators
1024  ST scalarl = 0.0;
1025  bool transpl = false;
1026  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1027  ST scalarm = 0.0;
1028  bool transpm = false;
1029  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1030  ST scalarr = 0.0;
1031  bool transpr = false;
1032  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1033 
1034  // Build output operator
1035  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1036  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1037 
1038  // Do explicit matrix-matrix multiply
1039  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1040  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1041  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1042  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1043  explicitCrsOp->resumeFill();
1044  explicitCrsOp->scale(scalarl*scalarm*scalarr);
1045  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1046  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1047  return tExplicitOp;
1048 
1049  } else if (isTpetral && !isTpetram && isTpetrar){ // Assume that the middle operator is diagonal
1050 
1051  // Get left and right Tpetra crs operators
1052  ST scalarl = 0.0;
1053  bool transpl = false;
1054  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1055  ST scalarr = 0.0;
1056  bool transpr = false;
1057  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1058 
1059  // Cast middle operator as DiagonalLinearOp and extract diagonal as Vector
1060  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm,true);
1061  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),true);
1062  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1063  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()));
1064 
1065  // Do the diagonal scaling
1066  tCrsOplm->rightScale(*diagPtr);
1067 
1068  // Build output operator
1069  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1070  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1071 
1072  // Do explicit matrix-matrix multiply
1073  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1074  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1075  explicitCrsOp->resumeFill();
1076  explicitCrsOp->scale(scalarl*scalarr);
1077  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1078  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1079  return tExplicitOp;
1080 
1081  } else { // Assume Epetra and we can use transformers
1082 
1083  // build implicit multiply
1084  const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1085 
1086  // build transformer
1087  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1088  Thyra::epetraExtDiagScaledMatProdTransformer();
1089 
1090  // build operator and multiply
1091  const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1092  prodTrans->transform(*implicitOp,explicitOp.ptr());
1093  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1094  " * " + opm->getObjectLabel() +
1095  " * " + opr->getObjectLabel() + " )");
1096 
1097  return explicitOp;
1098 
1099  }
1100 }
1101 
1116 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
1117  const ModifiableLinearOp & destOp)
1118 {
1119  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1120  bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1121  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1122 
1123  if(isTpetral && isTpetram && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1124 
1125  // Get left and right Tpetra crs operators
1126  ST scalarl = 0.0;
1127  bool transpl = false;
1128  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1129  ST scalarm = 0.0;
1130  bool transpm = false;
1131  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1132  ST scalarr = 0.0;
1133  bool transpr = false;
1134  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1135 
1136  // Build output operator
1137  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1138  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1139 
1140  // Do explicit matrix-matrix multiply
1141  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1142  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1143  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1144  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1145  explicitCrsOp->resumeFill();
1146  explicitCrsOp->scale(scalarl*scalarm*scalarr);
1147  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1148  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1149  return tExplicitOp;
1150 
1151  } else if (isTpetral && !isTpetram && isTpetrar){ // Assume that the middle operator is diagonal
1152 
1153  // Get left and right Tpetra crs operators
1154  ST scalarl = 0.0;
1155  bool transpl = false;
1156  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1157  ST scalarr = 0.0;
1158  bool transpr = false;
1159  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1160 
1161  // Cast middle operator as DiagonalLinearOp and extract diagonal as Vector
1162  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm,true);
1163  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),true);
1164  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1165  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()));
1166 
1167  // Do the diagonal scaling
1168  tCrsOplm->rightScale(*diagPtr);
1169 
1170  // Build output operator
1171  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1172  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1173 
1174  // Do explicit matrix-matrix multiply
1175  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1176  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1177  explicitCrsOp->resumeFill();
1178  explicitCrsOp->scale(scalarl*scalarr);
1179  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1180  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1181  return tExplicitOp;
1182 
1183  } else { // Assume Epetra and we can use transformers
1184 
1185  // build implicit multiply
1186  const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1187 
1188  // build transformer
1189  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1190  Thyra::epetraExtDiagScaledMatProdTransformer();
1191 
1192  // build operator destination operator
1193  ModifiableLinearOp explicitOp;
1194 
1195  // if neccessary build a operator to put the explicit multiply into
1196  if(destOp==Teuchos::null)
1197  explicitOp = prodTrans->createOutputOp();
1198  else
1199  explicitOp = destOp;
1200 
1201  // perform multiplication
1202  prodTrans->transform(*implicitOp,explicitOp.ptr());
1203 
1204  // label it
1205  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1206  " * " + opm->getObjectLabel() +
1207  " * " + opr->getObjectLabel() + " )");
1208 
1209  return explicitOp;
1210 
1211  }
1212 }
1213 
1224 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr)
1225 {
1226  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1227  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1228 
1229  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1230  // Get left and right Tpetra crs operators
1231  ST scalarl = 0.0;
1232  bool transpl = false;
1233  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1234  ST scalarr = 0.0;
1235  bool transpr = false;
1236  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1237 
1238  // Build output operator
1239  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1240  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1241 
1242  // Do explicit matrix-matrix multiply
1243  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1244  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1245  explicitCrsOp->resumeFill();
1246  explicitCrsOp->scale(scalarl*scalarr);
1247  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1248  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1249  return tExplicitOp;
1250 
1251  } else if (isTpetral && !isTpetrar){ // Assume that the right operator is diagonal
1252 
1253  // Get left Tpetra crs operator
1254  ST scalarl = 0.0;
1255  bool transpl = false;
1256  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1257 
1258  // Cast right operator as DiagonalLinearOp and extract diagonal as Vector
1259  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr);
1260  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),true);
1261  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1262  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()));
1263 
1264  explicitCrsOp->rightScale(*diagPtr);
1265  explicitCrsOp->resumeFill();
1266  explicitCrsOp->scale(scalarl);
1267  explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1268 
1269  return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1270 
1271  } else if (!isTpetral && isTpetrar){ // Assume that the left operator is diagonal
1272 
1273  // Get right Tpetra crs operator
1274  ST scalarr = 0.0;
1275  bool transpr = false;
1276  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1277 
1278  // Cast leftt operator as DiagonalLinearOp and extract diagonal as Vector
1279  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl);
1280  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),true);
1281  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1282  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()));
1283 
1284  explicitCrsOp->leftScale(*diagPtr);
1285  explicitCrsOp->resumeFill();
1286  explicitCrsOp->scale(scalarr);
1287  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1288 
1289  return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1290 
1291  } else { // Assume Epetra and we can use transformers
1292 
1293  // build implicit multiply
1294  const LinearOp implicitOp = Thyra::multiply(opl,opr);
1295 
1296  // build a scaling transformer
1297  RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1298  = Thyra::epetraExtDiagScalingTransformer();
1299 
1300  // check to see if a scaling transformer works: if not use the
1301  // DiagScaledMatrixProduct transformer
1302  if(not prodTrans->isCompatible(*implicitOp))
1303  prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1304 
1305  // build operator and multiply
1306  const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1307  prodTrans->transform(*implicitOp,explicitOp.ptr());
1308  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1309  " * " + opr->getObjectLabel() + " )");
1310 
1311  return explicitOp;
1312  }
1313 }
1314 
1328 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
1329  const ModifiableLinearOp & destOp)
1330 {
1331  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1332  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1333 
1334  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix multiply
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 scalarr = 0.0;
1341  bool transpr = false;
1342  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1343 
1344  // Build output operator
1345  RCP<Thyra::LinearOpBase<ST> > explicitOp;
1346  if(destOp!=Teuchos::null)
1347  explicitOp = destOp;
1348  else
1349  explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1350  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1351 
1352  // Do explicit matrix-matrix multiply
1353  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1354  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1355  explicitCrsOp->resumeFill();
1356  explicitCrsOp->scale(scalarl*scalarr);
1357  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1358  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1359  return tExplicitOp;
1360 
1361  } else if (isTpetral && !isTpetrar){ // Assume that the right operator is diagonal
1362 
1363  // Get left Tpetra crs operator
1364  ST scalarl = 0.0;
1365  bool transpl = false;
1366  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1367 
1368  // Cast right operator as DiagonalLinearOp and extract diagonal as Vector
1369  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr);
1370  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),true);
1371  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1372 
1373  // Scale by the diagonal operator
1374  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()));
1375  explicitCrsOp->rightScale(*diagPtr);
1376  explicitCrsOp->resumeFill();
1377  explicitCrsOp->scale(scalarl);
1378  explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1379  return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1380 
1381  } else if (!isTpetral && isTpetrar){ // Assume that the left operator is diagonal
1382 
1383  // Get right Tpetra crs operator
1384  ST scalarr = 0.0;
1385  bool transpr = false;
1386  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1387 
1388  // Cast leftt operator as DiagonalLinearOp and extract diagonal as Vector
1389  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl,true);
1390  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),true);
1391  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1392 
1393  // Scale by the diagonal operator
1394  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()));
1395  explicitCrsOp->leftScale(*diagPtr);
1396  explicitCrsOp->resumeFill();
1397  explicitCrsOp->scale(scalarr);
1398  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1399  return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1400 
1401  } else { // Assume Epetra and we can use transformers
1402 
1403  // build implicit multiply
1404  const LinearOp implicitOp = Thyra::multiply(opl,opr);
1405 
1406  // build a scaling transformer
1407 
1408  RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1409  = Thyra::epetraExtDiagScalingTransformer();
1410 
1411  // check to see if a scaling transformer works: if not use the
1412  // DiagScaledMatrixProduct transformer
1413  if(not prodTrans->isCompatible(*implicitOp))
1414  prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1415 
1416  // build operator destination operator
1417  ModifiableLinearOp explicitOp;
1418 
1419  // if neccessary build a operator to put the explicit multiply into
1420  if(destOp==Teuchos::null)
1421  explicitOp = prodTrans->createOutputOp();
1422  else
1423  explicitOp = destOp;
1424 
1425  // perform multiplication
1426  prodTrans->transform(*implicitOp,explicitOp.ptr());
1427 
1428  // label it
1429  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1430  " * " + opr->getObjectLabel() + " )");
1431 
1432  return explicitOp;
1433  }
1434 }
1435 
1446 const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr)
1447 {
1448  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1449  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1450 
1451  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix add
1452 
1453  // Get left and right Tpetra crs operators
1454  ST scalarl = 0.0;
1455  bool transpl = false;
1456  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1457  ST scalarr = 0.0;
1458  bool transpr = false;
1459  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1460 
1461  // Build output operator
1462  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1463  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1464 
1465  // Do explicit matrix-matrix add
1466  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1467  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1468  return tExplicitOp;
1469 
1470  }else{//Assume Epetra
1471  // build implicit add
1472  const LinearOp implicitOp = Thyra::add(opl,opr);
1473 
1474  // build transformer
1475  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1476  Thyra::epetraExtAddTransformer();
1477 
1478  // build operator and add
1479  const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1480  prodTrans->transform(*implicitOp,explicitOp.ptr());
1481  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1482  " + " + opr->getObjectLabel() + " )");
1483 
1484  return explicitOp;
1485  }
1486 }
1487 
1500 const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
1501  const ModifiableLinearOp & destOp)
1502 {
1503  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1504  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1505 
1506  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix add
1507 
1508  // Get left and right Tpetra crs operators
1509  ST scalarl = 0.0;
1510  bool transpl = false;
1511  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1512  ST scalarr = 0.0;
1513  bool transpr = false;
1514  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1515 
1516  // Build output operator
1517  RCP<Thyra::LinearOpBase<ST> > explicitOp;
1518  if(destOp!=Teuchos::null)
1519  explicitOp = destOp;
1520  else
1521  explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1522  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1523 
1524  // Do explicit matrix-matrix add
1525  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1526  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1527  return tExplicitOp;
1528 
1529  }else{ // Assume Epetra
1530 
1531  // build implicit add
1532  const LinearOp implicitOp = Thyra::add(opl,opr);
1533 
1534  // build transformer
1535  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1536  Thyra::epetraExtAddTransformer();
1537 
1538  // build or reuse destination operator
1539  RCP<Thyra::LinearOpBase<double> > explicitOp;
1540  if(destOp!=Teuchos::null)
1541  explicitOp = destOp;
1542  else
1543  explicitOp = prodTrans->createOutputOp();
1544 
1545  // perform add
1546  prodTrans->transform(*implicitOp,explicitOp.ptr());
1547  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1548  " + " + opr->getObjectLabel() + " )");
1549 
1550  return explicitOp;
1551  }
1552 }
1553 
1558 const ModifiableLinearOp explicitSum(const LinearOp & op,
1559  const ModifiableLinearOp & destOp)
1560 {
1561  // convert operators to Epetra_CrsMatrix
1562  const RCP<const Epetra_CrsMatrix> epetraOp =
1563  rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*op), true);
1564 
1565  if(destOp==Teuchos::null) {
1566  Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(new Epetra_CrsMatrix(*epetraOp));
1567 
1568  return Thyra::nonconstEpetraLinearOp(epetraDest);
1569  }
1570 
1571  const RCP<Epetra_CrsMatrix> epetraDest =
1572  rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp), true);
1573 
1574  EpetraExt::MatrixMatrix::Add(*epetraOp,false,1.0,*epetraDest,1.0);
1575 
1576  return destOp;
1577 }
1578 
1579 const LinearOp explicitTranspose(const LinearOp & op)
1580 {
1581  if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
1582 
1583  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op,true);
1584  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
1585 
1586  Tpetra::RowMatrixTransposer<ST,LO,GO,NT> transposer(tCrsOp);
1587  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > transOp = transposer.createTranspose();
1588 
1589  return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(transOp->getDomainMap()),transOp);
1590 
1591  } else {
1592 
1593  RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
1594  TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
1595  "Teko::explicitTranspose Not an Epetra_Operator");
1596  RCP<const Epetra_RowMatrix> eRowMatrixOp
1597  = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(eOp);
1598  TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
1599  "Teko::explicitTranspose Not an Epetra_RowMatrix");
1600 
1601  // we now have a delete transpose operator
1602  EpetraExt::RowMatrix_Transpose tranposeOp;
1603  Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
1604 
1605  // this copy is because of a poor implementation of the EpetraExt::Transform
1606  // implementation
1607  Teuchos::RCP<Epetra_CrsMatrix> crsMat
1608  = Teuchos::rcp(new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
1609 
1610  return Thyra::epetraLinearOp(crsMat);
1611  }
1612 }
1613 
1614 const double frobeniusNorm(const LinearOp & op)
1615 {
1616  if(Teko::TpetraHelpers::isTpetraLinearOp(op)){
1617  const RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
1618  const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
1619  return crsOp->getFrobeniusNorm();
1620  } else {
1621  const RCP<const Epetra_Operator> epOp = Thyra::get_Epetra_Operator(*op);
1622  const RCP<const Epetra_CrsMatrix> crsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(epOp,true);
1623  return crsOp->NormFrobenius();
1624  }
1625 }
1626 
1627 const LinearOp buildDiagonal(const MultiVector & src,const std::string & lbl)
1628 {
1629  RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range());
1630  Thyra::copy(*src->col(0),dst.ptr());
1631 
1632  return Thyra::diagonal<double>(dst,lbl);
1633 }
1634 
1635 const LinearOp buildInvDiagonal(const MultiVector & src,const std::string & lbl)
1636 {
1637  const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
1638  RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range());
1639  Thyra::reciprocal<double>(*srcV,dst.ptr());
1640 
1641  return Thyra::diagonal<double>(dst,lbl);
1642 }
1643 
1645 BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvv)
1646 {
1647  Teuchos::Array<MultiVector> mvA;
1648  Teuchos::Array<VectorSpace> vsA;
1649 
1650  // build arrays of multi vectors and vector spaces
1651  std::vector<MultiVector>::const_iterator itr;
1652  for(itr=mvv.begin();itr!=mvv.end();++itr) {
1653  mvA.push_back(*itr);
1654  vsA.push_back((*itr)->range());
1655  }
1656 
1657  // construct the product vector space
1658  const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
1659  = Thyra::productVectorSpace<double>(vsA);
1660 
1661  return Thyra::defaultProductMultiVector<double>(vs,mvA);
1662 }
1663 
1674 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
1675  const std::vector<int> & indices,
1676  const VectorSpace & vs,
1677  double onValue, double offValue)
1678 
1679 {
1680  using Teuchos::RCP;
1681 
1682  // create a new vector
1683  RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
1684  Thyra::put_scalar<double>(offValue,v.ptr()); // fill it with "off" values
1685 
1686  // set on values
1687  for(std::size_t i=0;i<indices.size();i++)
1688  Thyra::set_ele<double>(indices[i],onValue,v.ptr());
1689 
1690  return v;
1691 }
1692 
1717 double computeSpectralRad(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
1718  bool isHermitian,int numBlocks,int restart,int verbosity)
1719 {
1720  typedef Thyra::LinearOpBase<double> OP;
1721  typedef Thyra::MultiVectorBase<double> MV;
1722 
1723  int startVectors = 1;
1724 
1725  // construct an initial guess
1726  const RCP<MV> ivec = Thyra::createMember(A->domain());
1727  Thyra::randomize(-1.0,1.0,ivec.ptr());
1728 
1729  RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1730  = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1731  eigProb->setNEV(1);
1732  eigProb->setHermitian(isHermitian);
1733 
1734  // set the problem up
1735  if(not eigProb->setProblem()) {
1736  // big time failure!
1737  return Teuchos::ScalarTraits<double>::nan();
1738  }
1739 
1740  // we want largert magnitude eigenvalue
1741  std::string which("LM"); // largest magnitude
1742 
1743  // Create the parameter list for the eigensolver
1744  // verbosity+=Anasazi::TimingDetails;
1745  Teuchos::ParameterList MyPL;
1746  MyPL.set( "Verbosity", verbosity );
1747  MyPL.set( "Which", which );
1748  MyPL.set( "Block Size", startVectors );
1749  MyPL.set( "Num Blocks", numBlocks );
1750  MyPL.set( "Maximum Restarts", restart );
1751  MyPL.set( "Convergence Tolerance", tol );
1752 
1753  // build status test manager
1754  // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
1755  // = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(numBlocks*(restart+1)));
1756 
1757  // Create the Block Krylov Schur solver
1758  // This takes as inputs the eigenvalue problem and the solver parameters
1759  Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1760 
1761  // Solve the eigenvalue problem, and save the return code
1762  Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1763 
1764  if(solverreturn==Anasazi::Unconverged) {
1765  double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
1766  double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
1767 
1768  return -std::abs(std::sqrt(real*real+comp*comp));
1769 
1770  // cout << "Anasazi::BlockKrylovSchur::solve() did not converge!" << std::endl;
1771  // return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
1772  }
1773  else { // solverreturn==Anasazi::Converged
1774  double real = eigProb->getSolution().Evals.begin()->realpart;
1775  double comp = eigProb->getSolution().Evals.begin()->imagpart;
1776 
1777  return std::abs(std::sqrt(real*real+comp*comp));
1778 
1779  // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
1780  // return std::abs(eigProb->getSolution().Evals.begin()->realpart);
1781  }
1782 }
1783 
1807 double computeSmallestMagEig(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
1808  bool isHermitian,int numBlocks,int restart,int verbosity)
1809 {
1810  typedef Thyra::LinearOpBase<double> OP;
1811  typedef Thyra::MultiVectorBase<double> MV;
1812 
1813  int startVectors = 1;
1814 
1815  // construct an initial guess
1816  const RCP<MV> ivec = Thyra::createMember(A->domain());
1817  Thyra::randomize(-1.0,1.0,ivec.ptr());
1818 
1819  RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1820  = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1821  eigProb->setNEV(1);
1822  eigProb->setHermitian(isHermitian);
1823 
1824  // set the problem up
1825  if(not eigProb->setProblem()) {
1826  // big time failure!
1827  return Teuchos::ScalarTraits<double>::nan();
1828  }
1829 
1830  // we want largert magnitude eigenvalue
1831  std::string which("SM"); // smallest magnitude
1832 
1833  // Create the parameter list for the eigensolver
1834  Teuchos::ParameterList MyPL;
1835  MyPL.set( "Verbosity", verbosity );
1836  MyPL.set( "Which", which );
1837  MyPL.set( "Block Size", startVectors );
1838  MyPL.set( "Num Blocks", numBlocks );
1839  MyPL.set( "Maximum Restarts", restart );
1840  MyPL.set( "Convergence Tolerance", tol );
1841 
1842  // build status test manager
1843  // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
1844  // = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(10));
1845 
1846  // Create the Block Krylov Schur solver
1847  // This takes as inputs the eigenvalue problem and the solver parameters
1848  Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1849 
1850  // Solve the eigenvalue problem, and save the return code
1851  Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1852 
1853  if(solverreturn==Anasazi::Unconverged) {
1854  // cout << "Anasazi::BlockKrylovSchur::solve() did not converge! " << std::endl;
1855  return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
1856  }
1857  else { // solverreturn==Anasazi::Converged
1858  // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
1859  return std::abs(eigProb->getSolution().Evals.begin()->realpart);
1860  }
1861 }
1862 
1871 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt)
1872 {
1873  switch(dt) {
1874  case Diagonal:
1875  return getDiagonalOp(A);
1876  case Lumped:
1877  return getLumpedMatrix(A);
1878  case AbsRowSum:
1879  return getAbsRowSumMatrix(A);
1880  case NotDiag:
1881  default:
1882  TEUCHOS_TEST_FOR_EXCEPT(true);
1883  };
1884 
1885  return Teuchos::null;
1886 }
1887 
1896 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const Teko::DiagonalType & dt)
1897 {
1898  switch(dt) {
1899  case Diagonal:
1900  return getInvDiagonalOp(A);
1901  case Lumped:
1902  return getInvLumpedMatrix(A);
1903  case AbsRowSum:
1904  return getAbsRowSumInvMatrix(A);
1905  case NotDiag:
1906  default:
1907  TEUCHOS_TEST_FOR_EXCEPT(true);
1908  };
1909 
1910  return Teuchos::null;
1911 }
1912 
1919 std::string getDiagonalName(const DiagonalType & dt)
1920 {
1921  switch(dt) {
1922  case Diagonal:
1923  return "Diagonal";
1924  case Lumped:
1925  return "Lumped";
1926  case AbsRowSum:
1927  return "AbsRowSum";
1928  case NotDiag:
1929  return "NotDiag";
1930  case BlkDiag:
1931  return "BlkDiag";
1932  };
1933 
1934  return "<error>";
1935 }
1936 
1945 DiagonalType getDiagonalType(std::string name)
1946 {
1947  if(name=="Diagonal")
1948  return Diagonal;
1949  if(name=="Lumped")
1950  return Lumped;
1951  if(name=="AbsRowSum")
1952  return AbsRowSum;
1953  if(name=="BlkDiag")
1954  return BlkDiag;
1955 
1956  return NotDiag;
1957 }
1958 
1959 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,const LinearOp & Op){
1960 #ifdef Teko_ENABLE_Isorropia
1961  Teuchos::ParameterList probeList;
1962  Prober prober(G,probeList,true);
1963  Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(new Epetra_CrsMatrix(Copy,*G));
1965  prober.probe(Mwrap,*Mat);
1966  return Thyra::epetraLinearOp(Mat);
1967 #else
1968  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,"Probe requires Isorropia");
1969 #endif
1970 }
1971 
1972 double norm_1(const MultiVector & v,std::size_t col)
1973 {
1974  Teuchos::Array<double> n(v->domain()->dim());
1975  Thyra::norms_1<double>(*v,n);
1976 
1977  return n[col];
1978 }
1979 
1980 double norm_2(const MultiVector & v,std::size_t col)
1981 {
1982  Teuchos::Array<double> n(v->domain()->dim());
1983  Thyra::norms_2<double>(*v,n);
1984 
1985  return n[col];
1986 }
1987 
1988 void putScalar(const ModifiableLinearOp & op,double scalar)
1989 {
1990  try {
1991  // get Epetra_Operator
1992  RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
1993 
1994  // cast it to a CrsMatrix
1995  RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
1996 
1997  eCrsOp->PutScalar(scalar);
1998  }
1999  catch (std::exception & e) {
2000  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2001 
2002  *out << "Teko: putScalar requires an Epetra_CrsMatrix\n";
2003  *out << " Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
2004  *out << " OR\n";
2005  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
2006  *out << std::endl;
2007  *out << "*** THROWN EXCEPTION ***\n";
2008  *out << e.what() << std::endl;
2009  *out << "************************\n";
2010 
2011  throw e;
2012  }
2013 }
2014 
2015 void clipLower(MultiVector & v,double lowerBound)
2016 {
2017  using Teuchos::RCP;
2018  using Teuchos::rcp_dynamic_cast;
2019 
2020  // cast so entries are accessible
2021  // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2022  // = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
2023 
2024  for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2025  RCP<Thyra::SpmdVectorBase<double> > spmdVec
2026  = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2027 
2028  Teuchos::ArrayRCP<double> values;
2029  // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2030  spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2031  for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2032  if(values[j]<lowerBound)
2033  values[j] = lowerBound;
2034  }
2035  }
2036 }
2037 
2038 void clipUpper(MultiVector & v,double upperBound)
2039 {
2040  using Teuchos::RCP;
2041  using Teuchos::rcp_dynamic_cast;
2042 
2043  // cast so entries are accessible
2044  // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2045  // = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
2046  for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2047  RCP<Thyra::SpmdVectorBase<double> > spmdVec
2048  = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2049 
2050  Teuchos::ArrayRCP<double> values;
2051  // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2052  spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2053  for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2054  if(values[j]>upperBound)
2055  values[j] = upperBound;
2056  }
2057  }
2058 }
2059 
2060 void replaceValue(MultiVector & v,double currentValue,double newValue)
2061 {
2062  using Teuchos::RCP;
2063  using Teuchos::rcp_dynamic_cast;
2064 
2065  // cast so entries are accessible
2066  // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2067  // = rcp_dynamic_cast<Thyra::SpmdMultiVectorBase<double> >(v,true);
2068  for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2069  RCP<Thyra::SpmdVectorBase<double> > spmdVec
2070  = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2071 
2072  Teuchos::ArrayRCP<double> values;
2073  // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2074  spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2075  for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2076  if(values[j]==currentValue)
2077  values[j] = newValue;
2078  }
2079  }
2080 }
2081 
2082 void columnAverages(const MultiVector & v,std::vector<double> & averages)
2083 {
2084  averages.resize(v->domain()->dim());
2085 
2086  // sum over each column
2087  Thyra::sums<double>(*v,averages);
2088 
2089  // build averages
2090  Thyra::Ordinal rows = v->range()->dim();
2091  for(std::size_t i=0;i<averages.size();i++)
2092  averages[i] = averages[i]/rows;
2093 }
2094 
2095 double average(const MultiVector & v)
2096 {
2097  Thyra::Ordinal rows = v->range()->dim();
2098  Thyra::Ordinal cols = v->domain()->dim();
2099 
2100  std::vector<double> averages;
2101  columnAverages(v,averages);
2102 
2103  double sum = 0.0;
2104  for(std::size_t i=0;i<averages.size();i++)
2105  sum += averages[i] * rows;
2106 
2107  return sum/(rows*cols);
2108 }
2109 
2110 }
Specifies that a block diagonal approximation is to be used.
Specifies that just the diagonal is used.
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
Specifies that row sum is used to form a diagonal.
void scale(const double alpha, MultiVector &x)
Scale a multivector by a constant.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...
For user convenience, if Teko recieves this value, exceptions will be thrown.
DiagonalType
Type describing the type of diagonal to construct.
int blockColCount(const BlockedLinearOp &blo)
Get the column count in a block linear operator.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
Specifies that the diagonal entry is .