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 "Epetra/Teko_EpetraHelpers.hpp"
98 #include "Epetra/Teko_EpetraOperatorWrapper.hpp"
99 #include "Tpetra/Teko_TpetraHelpers.hpp"
100 #include "Tpetra/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 
109 #include <cmath>
110 
111 namespace Teko {
112 
113 using Teuchos::rcp;
114 using Teuchos::rcp_dynamic_cast;
115 using Teuchos::RCP;
116 #ifdef Teko_ENABLE_Isorropia
117 using Isorropia::Epetra::Prober;
118 #endif
119 
120 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
121 {
122  Teuchos::RCP<Teuchos::FancyOStream> os =
123  Teuchos::VerboseObjectBase::getDefaultOStream();
124 
125  //os->setShowProcRank(true);
126  //os->setOutputToRootOnly(-1);
127  return os;
128 }
129 
130 // distance function...not parallel...entirely internal to this cpp file
131 inline double dist(int dim,double * coords,int row,int col)
132 {
133  double value = 0.0;
134  for(int i=0;i<dim;i++)
135  value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
136 
137  // the distance between the two
138  return std::sqrt(value);
139 }
140 
141 // distance function...not parallel...entirely internal to this cpp file
142 inline double dist(double * x,double * y,double * z,int stride,int row,int col)
143 {
144  double value = 0.0;
145  if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
146  if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
147  if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
148 
149  // the distance between the two
150  return std::sqrt(value);
151 }
152 
171 RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil)
172 {
173  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
174  RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
175  stencil.MaxNumEntries()+1),true);
176 
177  // allocate an additional value for the diagonal, if neccessary
178  std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
179  std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
180 
181  // loop over all the rows
182  for(int j=0;j<gl->NumMyRows();j++) {
183  int row = gl->GRID(j);
184  double diagValue = 0.0;
185  int diagInd = -1;
186  int rowSz = 0;
187 
188  // extract a copy of this row...put it in rowData, rowIndicies
189  stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
190 
191  // loop over elements of row
192  for(int i=0;i<rowSz;i++) {
193  int col = rowInd[i];
194 
195  // is this a 0 entry masquerading as some thing else?
196  double value = rowData[i];
197  if(value==0) continue;
198 
199  // for nondiagonal entries
200  if(row!=col) {
201  double d = dist(dim,coords,row,col);
202  rowData[i] = -1.0/d;
203  diagValue += rowData[i];
204  }
205  else
206  diagInd = i;
207  }
208 
209  // handle diagonal entry
210  if(diagInd<0) { // diagonal not in row
211  rowData[rowSz] = -diagValue;
212  rowInd[rowSz] = row;
213  rowSz++;
214  }
215  else { // diagonal in row
216  rowData[diagInd] = -diagValue;
217  rowInd[diagInd] = row;
218  }
219 
220  // insert row data into graph Laplacian matrix
221  TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
222  }
223 
224  gl->FillComplete();
225 
226  return gl;
227 }
228 
229 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(int dim,ST * coords,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
230 {
231  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
232  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
233  stencil.getGlobalMaxNumRowEntries()+1));
234 
235  // allocate an additional value for the diagonal, if neccessary
236  std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
237  std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
238 
239  // loop over all the rows
240  for(LO j=0;j< (LO) gl->getNodeNumRows();j++) {
241  GO row = gl->getRowMap()->getGlobalElement(j);
242  ST diagValue = 0.0;
243  GO diagInd = -1;
244  size_t rowSz = 0;
245 
246  // extract a copy of this row...put it in rowData, rowIndicies
247  stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
248 
249  // loop over elements of row
250  for(size_t i=0;i<rowSz;i++) {
251  GO col = rowInd[i];
252 
253  // is this a 0 entry masquerading as some thing else?
254  ST value = rowData[i];
255  if(value==0) continue;
256 
257  // for nondiagonal entries
258  if(row!=col) {
259  ST d = dist(dim,coords,row,col);
260  rowData[i] = -1.0/d;
261  diagValue += rowData[i];
262  }
263  else
264  diagInd = i;
265  }
266 
267  // handle diagonal entry
268  if(diagInd<0) { // diagonal not in row
269  rowData[rowSz] = -diagValue;
270  rowInd[rowSz] = row;
271  rowSz++;
272  }
273  else { // diagonal in row
274  rowData[diagInd] = -diagValue;
275  rowInd[diagInd] = row;
276  }
277 
278  // insert row data into graph Laplacian matrix
279  gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
280  }
281 
282  gl->fillComplete();
283 
284  return gl;
285 }
286 
309 RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil)
310 {
311  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
312  RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
313  stencil.MaxNumEntries()+1),true);
314 
315  // allocate an additional value for the diagonal, if neccessary
316  std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
317  std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
318 
319  // loop over all the rows
320  for(int j=0;j<gl->NumMyRows();j++) {
321  int row = gl->GRID(j);
322  double diagValue = 0.0;
323  int diagInd = -1;
324  int rowSz = 0;
325 
326  // extract a copy of this row...put it in rowData, rowIndicies
327  stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
328 
329  // loop over elements of row
330  for(int i=0;i<rowSz;i++) {
331  int col = rowInd[i];
332 
333  // is this a 0 entry masquerading as some thing else?
334  double value = rowData[i];
335  if(value==0) continue;
336 
337  // for nondiagonal entries
338  if(row!=col) {
339  double d = dist(x,y,z,stride,row,col);
340  rowData[i] = -1.0/d;
341  diagValue += rowData[i];
342  }
343  else
344  diagInd = i;
345  }
346 
347  // handle diagonal entry
348  if(diagInd<0) { // diagonal not in row
349  rowData[rowSz] = -diagValue;
350  rowInd[rowSz] = row;
351  rowSz++;
352  }
353  else { // diagonal in row
354  rowData[diagInd] = -diagValue;
355  rowInd[diagInd] = row;
356  }
357 
358  // insert row data into graph Laplacian matrix
359  TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
360  }
361 
362  gl->FillComplete();
363 
364  return gl;
365 }
366 
367 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > buildGraphLaplacian(ST * x,ST * y,ST * z,GO stride,const Tpetra::CrsMatrix<ST,LO,GO,NT> & stencil)
368 {
369  // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
370  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > gl = rcp(new Tpetra::CrsMatrix<ST,LO,GO,NT>(stencil.getRowMap(),stencil.getColMap(),
371  stencil.getGlobalMaxNumRowEntries()+1),true);
372 
373  // allocate an additional value for the diagonal, if neccessary
374  std::vector<ST> rowData(stencil.getGlobalMaxNumRowEntries()+1);
375  std::vector<GO> rowInd(stencil.getGlobalMaxNumRowEntries()+1);
376 
377  // loop over all the rows
378  for(LO j=0;j<(LO) gl->getNodeNumRows();j++) {
379  GO row = gl->getRowMap()->getGlobalElement(j);
380  ST diagValue = 0.0;
381  GO diagInd = -1;
382  size_t rowSz = 0;
383 
384  // extract a copy of this row...put it in rowData, rowIndicies
385  stencil.getGlobalRowCopy(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData),rowSz);
386 
387  // loop over elements of row
388  for(size_t i=0;i<rowSz;i++) {
389  GO col = rowInd[i];
390 
391  // is this a 0 entry masquerading as some thing else?
392  ST value = rowData[i];
393  if(value==0) continue;
394 
395  // for nondiagonal entries
396  if(row!=col) {
397  ST d = dist(x,y,z,stride,row,col);
398  rowData[i] = -1.0/d;
399  diagValue += rowData[i];
400  }
401  else
402  diagInd = i;
403  }
404 
405  // handle diagonal entry
406  if(diagInd<0) { // diagonal not in row
407  rowData[rowSz] = -diagValue;
408  rowInd[rowSz] = row;
409  rowSz++;
410  }
411  else { // diagonal in row
412  rowData[diagInd] = -diagValue;
413  rowInd[diagInd] = row;
414  }
415 
416  // insert row data into graph Laplacian matrix
417  gl->replaceGlobalValues(row,Teuchos::ArrayView<GO>(rowInd),Teuchos::ArrayView<ST>(rowData));
418  }
419 
420  gl->fillComplete();
421 
422  return gl;
423 }
424 
440 void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
441 {
442  // Thyra::apply(*A,Thyra::NONCONJ_ELE,*x,&*y,alpha,beta);
443  Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
444 }
445 
448 void update(double alpha,const MultiVector & x,double beta,MultiVector & y)
449 {
450  Teuchos::Array<double> scale;
451  Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > > vec;
452 
453  // build arrays needed for linear combo
454  scale.push_back(alpha);
455  vec.push_back(x.ptr());
456 
457  // compute linear combination
458  Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
459 }
460 
462 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
463 {
464  int rows = blockRowCount(blo);
465 
466  TEUCHOS_ASSERT(rows==blockColCount(blo));
467 
468  RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
469  RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
470 
471  // allocate new operator
472  BlockedLinearOp upper = createBlockedOp();
473 
474  // build new operator
475  upper->beginBlockFill(rows,rows);
476 
477  for(int i=0;i<rows;i++) {
478  // put zero operators on the diagonal
479  // this gurantees the vector space of
480  // the new operator are fully defined
481  RCP<const Thyra::LinearOpBase<double> > zed
482  = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
483  upper->setBlock(i,i,zed);
484 
485  for(int j=i+1;j<rows;j++) {
486  // get block i,j
487  LinearOp uij = blo->getBlock(i,j);
488 
489  // stuff it in U
490  if(uij!=Teuchos::null)
491  upper->setBlock(i,j,uij);
492  }
493  }
494  if(callEndBlockFill)
495  upper->endBlockFill();
496 
497  return upper;
498 }
499 
501 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
502 {
503  int rows = blockRowCount(blo);
504 
505  TEUCHOS_ASSERT(rows==blockColCount(blo));
506 
507  RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
508  RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
509 
510  // allocate new operator
511  BlockedLinearOp lower = createBlockedOp();
512 
513  // build new operator
514  lower->beginBlockFill(rows,rows);
515 
516  for(int i=0;i<rows;i++) {
517  // put zero operators on the diagonal
518  // this gurantees the vector space of
519  // the new operator are fully defined
520  RCP<const Thyra::LinearOpBase<double> > zed
521  = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
522  lower->setBlock(i,i,zed);
523 
524  for(int j=0;j<i;j++) {
525  // get block i,j
526  LinearOp uij = blo->getBlock(i,j);
527 
528  // stuff it in U
529  if(uij!=Teuchos::null)
530  lower->setBlock(i,j,uij);
531  }
532  }
533  if(callEndBlockFill)
534  lower->endBlockFill();
535 
536  return lower;
537 }
538 
558 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo)
559 {
560  int rows = blockRowCount(blo);
561 
562  TEUCHOS_ASSERT(rows==blockColCount(blo)); // assert that matrix is square
563 
564  RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
565  RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
566 
567  // allocate new operator
568  BlockedLinearOp zeroOp = createBlockedOp();
569 
570  // build new operator
571  zeroOp->beginBlockFill(rows,rows);
572 
573  for(int i=0;i<rows;i++) {
574  // put zero operators on the diagonal
575  // this gurantees the vector space of
576  // the new operator are fully defined
577  RCP<const Thyra::LinearOpBase<double> > zed
578  = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
579  zeroOp->setBlock(i,i,zed);
580  }
581 
582  return zeroOp;
583 }
584 
586 bool isZeroOp(const LinearOp op)
587 {
588  // if operator is null...then its zero!
589  if(op==Teuchos::null) return true;
590 
591  // try to cast it to a zero linear operator
592  const LinearOp test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(op);
593 
594  // if it works...then its zero...otherwise its null
595  return test!=Teuchos::null;
596 }
597 
606 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op)
607 {
608  bool isTpetra = false;
609  RCP<const Epetra_CrsMatrix> eCrsOp;
610  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
611 
612  try {
613  // get Epetra or Tpetra Operator
614  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
615  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
616 
617  // cast it to a CrsMatrix
618  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
619  if (!eOp.is_null()){
620  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
621  }
622  else if (!tOp.is_null()){
623  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
624  isTpetra = true;
625  }
626  else
627  throw std::logic_error("Neither Epetra nor Tpetra");
628  }
629  catch (std::exception & e) {
630  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
631 
632  *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
633  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
634  *out << " OR\n";
635  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
636  *out << std::endl;
637  *out << "*** THROWN EXCEPTION ***\n";
638  *out << e.what() << std::endl;
639  *out << "************************\n";
640 
641  throw e;
642  }
643 
644  if(!isTpetra){
645  // extract diagonal
646  const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
647  Epetra_Vector & diag = *ptrDiag;
648 
649  // compute absolute value row sum
650  diag.PutScalar(0.0);
651  for(int i=0;i<eCrsOp->NumMyRows();i++) {
652  double * values = 0;
653  int numEntries;
654  eCrsOp->ExtractMyRowView(i,numEntries,values);
655 
656  // build abs value row sum
657  for(int j=0;j<numEntries;j++)
658  diag[i] += std::abs(values[j]);
659  }
660 
661  // build Thyra diagonal operator
662  return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
663  }
664  else {
665  // extract diagonal
666  const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
667  Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
668 
669  // compute absolute value row sum
670  diag.putScalar(0.0);
671  for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
672  LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
673  std::vector<LO> indices(numEntries);
674  std::vector<ST> values(numEntries);
675  Teuchos::ArrayView<const LO> indices_av(indices);
676  Teuchos::ArrayView<const ST> values_av(values);
677  tCrsOp->getLocalRowView(i,indices_av,values_av);
678 
679  // build abs value row sum
680  for(LO j=0;j<numEntries;j++)
681  diag.sumIntoLocalValue(i,std::abs(values_av[j]));
682  }
683 
684  // build Thyra diagonal operator
685  return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),"absRowSum( " + op->getObjectLabel() + " ))");
686 
687  }
688 
689 }
690 
699 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op)
700 {
701  bool isTpetra = false;
702  RCP<const Epetra_CrsMatrix> eCrsOp;
703  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
704 
705  try {
706  // get Epetra or Tpetra Operator
707  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
708  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
709 
710  // cast it to a CrsMatrix
711  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
712  if (!eOp.is_null()){
713  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
714  }
715  else if (!tOp.is_null()){
716  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
717  isTpetra = true;
718  }
719  else
720  throw std::logic_error("Neither Epetra nor Tpetra");
721  }
722  catch (std::exception & e) {
723  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
724 
725  *out << "Teko: getAbsRowSumInvMatrix requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
726  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
727  *out << " OR\n";
728  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
729  *out << std::endl;
730  *out << "*** THROWN EXCEPTION ***\n";
731  *out << e.what() << std::endl;
732  *out << "************************\n";
733 
734  throw e;
735  }
736 
737  if(!isTpetra){
738  // extract diagonal
739  const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
740  Epetra_Vector & diag = *ptrDiag;
741 
742  // compute absolute value row sum
743  diag.PutScalar(0.0);
744  for(int i=0;i<eCrsOp->NumMyRows();i++) {
745  double * values = 0;
746  int numEntries;
747  eCrsOp->ExtractMyRowView(i,numEntries,values);
748 
749  // build abs value row sum
750  for(int j=0;j<numEntries;j++)
751  diag[i] += std::abs(values[j]);
752  }
753  diag.Reciprocal(diag); // invert entries
754 
755  // build Thyra diagonal operator
756  return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
757  }
758  else {
759  // extract diagonal
760  const RCP<Tpetra::Vector<ST,LO,GO,NT> > ptrDiag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
761  Tpetra::Vector<ST,LO,GO,NT> & diag = *ptrDiag;
762 
763  // compute absolute value row sum
764  diag.putScalar(0.0);
765  for(LO i=0;i<(LO) tCrsOp->getNodeNumRows();i++) {
766  LO numEntries = tCrsOp->getNumEntriesInLocalRow (i);
767  std::vector<LO> indices(numEntries);
768  std::vector<ST> values(numEntries);
769  Teuchos::ArrayView<const LO> indices_av(indices);
770  Teuchos::ArrayView<const ST> values_av(values);
771  tCrsOp->getLocalRowView(i,indices_av,values_av);
772 
773  // build abs value row sum
774  for(LO j=0;j<numEntries;j++)
775  diag.sumIntoLocalValue(i,std::abs(values_av[j]));
776  }
777  diag.reciprocal(diag); // invert entries
778 
779  // build Thyra diagonal operator
780  return Teko::TpetraHelpers::thyraDiagOp(ptrDiag,*tCrsOp->getRowMap(),"absRowSum( " + op->getObjectLabel() + " ))");
781 
782  }
783 
784 }
785 
793 ModifiableLinearOp getLumpedMatrix(const LinearOp & op)
794 {
795  RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
796  RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
797 
798  // set to all ones
799  Thyra::assign(ones.ptr(),1.0);
800 
801  // compute lumped diagonal
802  // Thyra::apply(*op,Thyra::NONCONJ_ELE,*ones,&*diag);
803  Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
804 
805  return rcp(new Thyra::DefaultDiagonalLinearOp<ST>(diag));
806 }
807 
816 ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op)
817 {
818  RCP<Thyra::VectorBase<ST> > ones = Thyra::createMember(op->domain());
819  RCP<Thyra::VectorBase<ST> > diag = Thyra::createMember(op->range());
820 
821  // set to all ones
822  Thyra::assign(ones.ptr(),1.0);
823 
824  // compute lumped diagonal
825  Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
826  Thyra::reciprocal(*diag,diag.ptr());
827 
828  return rcp(new Thyra::DefaultDiagonalLinearOp<ST>(diag));
829 }
830 
842 const ModifiableLinearOp getDiagonalOp(const LinearOp & op)
843 {
844  bool isTpetra = false;
845  RCP<const Epetra_CrsMatrix> eCrsOp;
846  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
847 
848  try {
849  // get Epetra or Tpetra Operator
850  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
851  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
852 
853  // cast it to a CrsMatrix
854  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
855  if (!eOp.is_null()){
856  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
857  }
858  else if (!tOp.is_null()){
859  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
860  isTpetra = true;
861  }
862  else
863  throw std::logic_error("Neither Epetra nor Tpetra");
864  }
865  catch (std::exception & e) {
866  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
867 
868  *out << "Teko: getDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
869  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
870  *out << " OR\n";
871  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
872  *out << std::endl;
873  *out << "*** THROWN EXCEPTION ***\n";
874  *out << e.what() << std::endl;
875  *out << "************************\n";
876 
877  throw e;
878  }
879 
880  if(!isTpetra){
881  // extract diagonal
882  const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
883  TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
884 
885  // build Thyra diagonal operator
886  return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
887  }
888  else {
889  // extract diagonal
890  const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
891  tCrsOp->getLocalDiagCopy(*diag);
892 
893  // build Thyra diagonal operator
894  return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
895 
896  }
897 }
898 
899 const MultiVector getDiagonal(const LinearOp & op)
900 {
901  bool isTpetra = false;
902  RCP<const Epetra_CrsMatrix> eCrsOp;
903  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
904 
905  try {
906  // get Epetra or Tpetra Operator
907  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
908  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
909 
910  // cast it to a CrsMatrix
911  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
912  if (!eOp.is_null()){
913  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
914  }
915  else if (!tOp.is_null()){
916  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
917  isTpetra = true;
918  }
919  else
920  throw std::logic_error("Neither Epetra nor Tpetra");
921  }
922  catch (std::exception & e) {
923  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
924 
925  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
926  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
927  *out << eOp;
928  *out << tOp;
929 
930  *out << "Teko: getDiagonal requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
931  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
932  *out << " OR\n";
933  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
934  *out << std::endl;
935  *out << "*** THROWN EXCEPTION ***\n";
936  *out << e.what() << std::endl;
937  *out << "************************\n";
938 
939  throw e;
940  }
941 
942  if(!isTpetra){
943  // extract diagonal
944  const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
945  TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
946 
947  return Thyra::create_Vector(diag,Thyra::create_VectorSpace(Teuchos::rcpFromRef(eCrsOp->RowMap())));
948  }
949  else {
950  // extract diagonal
951  const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
952  tCrsOp->getLocalDiagCopy(*diag);
953 
954  // build Thyra diagonal operator
955  return Thyra::createVector<ST,LO,GO,NT>(diag,Thyra::createVectorSpace<ST,LO,GO,NT>(tCrsOp->getRowMap()));
956 
957  }
958 }
959 
960 const MultiVector getDiagonal(const Teko::LinearOp & A,const DiagonalType & dt)
961 {
962  LinearOp diagOp = Teko::getDiagonalOp(A,dt);
963 
964  Teuchos::RCP<const Thyra::MultiVectorBase<double> > v =
965  Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag();
966  return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
967 }
968 
980 const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op)
981 {
982  bool isTpetra = false;
983  RCP<const Epetra_CrsMatrix> eCrsOp;
984  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOp;
985 
986  try {
987  // get Epetra or Tpetra Operator
988  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
989  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
990 
991  // cast it to a CrsMatrix
992  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
993  if (!eOp.is_null()){
994  eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp->epetra_op(),true);
995  }
996  else if (!tOp.is_null()){
997  tCrsOp = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator(),true);
998  isTpetra = true;
999  }
1000  else
1001  throw std::logic_error("Neither Epetra nor Tpetra");
1002  }
1003  catch (std::exception & e) {
1004  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
1005 
1006  RCP<const Thyra::EpetraLinearOp > eOp = rcp_dynamic_cast<const Thyra::EpetraLinearOp >(op);
1007  RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(op);
1008  *out << eOp;
1009  *out << tOp;
1010 
1011  *out << "Teko: getInvDiagonalOp requires an Epetra_CrsMatrix or a Tpetra::CrsMatrix\n";
1012  *out << " Could not extract an Epetra_Operator or a Tpetra_Operator from a \"" << op->description() << std::endl;
1013  *out << " OR\n";
1014  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix or a Tpetra_Operator to a Tpetra::CrsMatrix\n";
1015  *out << std::endl;
1016  *out << "*** THROWN EXCEPTION ***\n";
1017  *out << e.what() << std::endl;
1018  *out << "************************\n";
1019 
1020  throw e;
1021  }
1022 
1023  if(!isTpetra){
1024  // extract diagonal
1025  const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
1026  TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
1027  diag->Reciprocal(*diag);
1028 
1029  // build Thyra diagonal operator
1030  return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
1031  }
1032  else {
1033  // extract diagonal
1034  const RCP<Tpetra::Vector<ST,LO,GO,NT> > diag = Tpetra::createVector<ST,LO,GO,NT>(tCrsOp->getRowMap());
1035  tCrsOp->getLocalDiagCopy(*diag);
1036  diag->reciprocal(*diag);
1037 
1038  // build Thyra diagonal operator
1039  return Teko::TpetraHelpers::thyraDiagOp(diag,*tCrsOp->getRowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
1040 
1041  }
1042 }
1043 
1056 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr)
1057 {
1058  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1059  bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1060  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1061 
1062  if(isTpetral && isTpetram && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1063 
1064  // Get left and right Tpetra crs operators
1065  ST scalarl = 0.0;
1066  bool transpl = false;
1067  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1068  ST scalarm = 0.0;
1069  bool transpm = false;
1070  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1071  ST scalarr = 0.0;
1072  bool transpr = false;
1073  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1074 
1075  // Build output operator
1076  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1077  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1078 
1079  // Do explicit matrix-matrix multiply
1080  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1081  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1082  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1083  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1084  explicitCrsOp->resumeFill();
1085  explicitCrsOp->scale(scalarl*scalarm*scalarr);
1086  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1087  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1088  return tExplicitOp;
1089 
1090  } else if (isTpetral && !isTpetram && isTpetrar){ // Assume that the middle operator is diagonal
1091 
1092  // Get left and right Tpetra crs operators
1093  ST scalarl = 0.0;
1094  bool transpl = false;
1095  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1096  ST scalarr = 0.0;
1097  bool transpr = false;
1098  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1099 
1100  // Cast middle operator as DiagonalLinearOp and extract diagonal as Vector
1101  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm,true);
1102  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),true);
1103  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1104  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()));
1105 
1106  // Do the diagonal scaling
1107  tCrsOplm->rightScale(*diagPtr);
1108 
1109  // Build output operator
1110  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1111  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1112 
1113  // Do explicit matrix-matrix multiply
1114  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1115  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1116  explicitCrsOp->resumeFill();
1117  explicitCrsOp->scale(scalarl*scalarr);
1118  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1119  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1120  return tExplicitOp;
1121 
1122  } else { // Assume Epetra and we can use transformers
1123 
1124  // build implicit multiply
1125  const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1126 
1127  // build transformer
1128  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1129  Thyra::epetraExtDiagScaledMatProdTransformer();
1130 
1131  // build operator and multiply
1132  const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1133  prodTrans->transform(*implicitOp,explicitOp.ptr());
1134  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1135  " * " + opm->getObjectLabel() +
1136  " * " + opr->getObjectLabel() + " )");
1137 
1138  return explicitOp;
1139 
1140  }
1141 }
1142 
1157 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
1158  const ModifiableLinearOp & destOp)
1159 {
1160  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1161  bool isTpetram = Teko::TpetraHelpers::isTpetraLinearOp(opm);
1162  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1163 
1164  if(isTpetral && isTpetram && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1165 
1166  // Get left and right Tpetra crs operators
1167  ST scalarl = 0.0;
1168  bool transpl = false;
1169  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1170  ST scalarm = 0.0;
1171  bool transpm = false;
1172  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpm = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarm, &transpm);
1173  ST scalarr = 0.0;
1174  bool transpr = false;
1175  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1176 
1177  // Build output operator
1178  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1179  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1180 
1181  // Do explicit matrix-matrix multiply
1182  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOplm = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1183  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1184  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpm,transpm,*tCrsOplm);
1185  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1186  explicitCrsOp->resumeFill();
1187  explicitCrsOp->scale(scalarl*scalarm*scalarr);
1188  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1189  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1190  return tExplicitOp;
1191 
1192  } else if (isTpetral && !isTpetram && isTpetrar){ // Assume that the middle operator is diagonal
1193 
1194  // Get left and right Tpetra crs operators
1195  ST scalarl = 0.0;
1196  bool transpl = false;
1197  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1198  ST scalarr = 0.0;
1199  bool transpr = false;
1200  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1201 
1202  // Cast middle operator as DiagonalLinearOp and extract diagonal as Vector
1203  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpm = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opm,true);
1204  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpm->getDiag(),true);
1205  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1206  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()));
1207 
1208  // Do the diagonal scaling
1209  tCrsOplm->rightScale(*diagPtr);
1210 
1211  // Build output operator
1212  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1213  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1214 
1215  // Do explicit matrix-matrix multiply
1216  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1217  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOplm,false,*tCrsOpr,transpr,*explicitCrsOp);
1218  explicitCrsOp->resumeFill();
1219  explicitCrsOp->scale(scalarl*scalarr);
1220  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1221  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1222  return tExplicitOp;
1223 
1224  } else { // Assume Epetra and we can use transformers
1225 
1226  // build implicit multiply
1227  const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
1228 
1229  // build transformer
1230  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1231  Thyra::epetraExtDiagScaledMatProdTransformer();
1232 
1233  // build operator destination operator
1234  ModifiableLinearOp explicitOp;
1235 
1236  // if neccessary build a operator to put the explicit multiply into
1237  if(destOp==Teuchos::null)
1238  explicitOp = prodTrans->createOutputOp();
1239  else
1240  explicitOp = destOp;
1241 
1242  // perform multiplication
1243  prodTrans->transform(*implicitOp,explicitOp.ptr());
1244 
1245  // label it
1246  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1247  " * " + opm->getObjectLabel() +
1248  " * " + opr->getObjectLabel() + " )");
1249 
1250  return explicitOp;
1251 
1252  }
1253 }
1254 
1265 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr)
1266 {
1267  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1268  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1269 
1270  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so explicitly multiply them
1271  // Get left and right Tpetra crs operators
1272  ST scalarl = 0.0;
1273  bool transpl = false;
1274  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1275  ST scalarr = 0.0;
1276  bool transpr = false;
1277  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1278 
1279  // Build output operator
1280  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1281  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1282 
1283  // Do explicit matrix-matrix multiply
1284  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1285  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1286  explicitCrsOp->resumeFill();
1287  explicitCrsOp->scale(scalarl*scalarr);
1288  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1289  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1290  return tExplicitOp;
1291 
1292  } else if (isTpetral && !isTpetrar){ // Assume that the right operator is diagonal
1293 
1294  // Get left Tpetra crs operator
1295  ST scalarl = 0.0;
1296  bool transpl = false;
1297  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1298 
1299  // Cast right operator as DiagonalLinearOp and extract diagonal as Vector
1300  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr);
1301  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),true);
1302  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1303  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()));
1304 
1305  explicitCrsOp->rightScale(*diagPtr);
1306  explicitCrsOp->resumeFill();
1307  explicitCrsOp->scale(scalarl);
1308  explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1309 
1310  return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1311 
1312  } else if (!isTpetral && isTpetrar){ // Assume that the left operator is diagonal
1313 
1314  // Get right Tpetra crs operator
1315  ST scalarr = 0.0;
1316  bool transpr = false;
1317  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1318 
1319  // Cast leftt operator as DiagonalLinearOp and extract diagonal as Vector
1320  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl);
1321  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),true);
1322  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1323  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()));
1324 
1325  explicitCrsOp->leftScale(*diagPtr);
1326  explicitCrsOp->resumeFill();
1327  explicitCrsOp->scale(scalarr);
1328  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1329 
1330  return Thyra::constTpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1331 
1332  } else { // Assume Epetra and we can use transformers
1333 
1334  // build implicit multiply
1335  const LinearOp implicitOp = Thyra::multiply(opl,opr);
1336 
1337  // build a scaling transformer
1338  RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1339  = Thyra::epetraExtDiagScalingTransformer();
1340 
1341  // check to see if a scaling transformer works: if not use the
1342  // DiagScaledMatrixProduct transformer
1343  if(not prodTrans->isCompatible(*implicitOp))
1344  prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1345 
1346  // build operator and multiply
1347  const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1348  prodTrans->transform(*implicitOp,explicitOp.ptr());
1349  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1350  " * " + opr->getObjectLabel() + " )");
1351 
1352  return explicitOp;
1353  }
1354 }
1355 
1369 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
1370  const ModifiableLinearOp & destOp)
1371 {
1372  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1373  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1374 
1375  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix multiply
1376 
1377  // Get left and right Tpetra crs operators
1378  ST scalarl = 0.0;
1379  bool transpl = false;
1380  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1381  ST scalarr = 0.0;
1382  bool transpr = false;
1383  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1384 
1385  // Build output operator
1386  RCP<Thyra::LinearOpBase<ST> > explicitOp;
1387  if(destOp!=Teuchos::null)
1388  explicitOp = destOp;
1389  else
1390  explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1391  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1392 
1393  // Do explicit matrix-matrix multiply
1394  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::createCrsMatrix<ST,LO,GO,NT>(tCrsOpl->getRowMap());
1395  Tpetra::MatrixMatrix::Multiply<ST,LO,GO,NT>(*tCrsOpl,transpl,*tCrsOpr,transpr,*explicitCrsOp);
1396  explicitCrsOp->resumeFill();
1397  explicitCrsOp->scale(scalarl*scalarr);
1398  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpl->getRangeMap());
1399  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1400  return tExplicitOp;
1401 
1402  } else if (isTpetral && !isTpetrar){ // Assume that the right operator is diagonal
1403 
1404  // Get left Tpetra crs operator
1405  ST scalarl = 0.0;
1406  bool transpl = false;
1407  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1408 
1409  // Cast right operator as DiagonalLinearOp and extract diagonal as Vector
1410  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpr = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opr);
1411  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpr->getDiag(),true);
1412  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1413 
1414  // Scale by the diagonal operator
1415  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()));
1416  explicitCrsOp->rightScale(*diagPtr);
1417  explicitCrsOp->resumeFill();
1418  explicitCrsOp->scale(scalarl);
1419  explicitCrsOp->fillComplete(tCrsOpl->getDomainMap(),tCrsOpl->getRangeMap());
1420  return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1421 
1422  } else if (!isTpetral && isTpetrar){ // Assume that the left operator is diagonal
1423 
1424  // Get right Tpetra crs operator
1425  ST scalarr = 0.0;
1426  bool transpr = false;
1427  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1428 
1429  // Cast leftt operator as DiagonalLinearOp and extract diagonal as Vector
1430  RCP<const Thyra::DiagonalLinearOpBase<ST> > dOpl = rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<ST> >(opl,true);
1431  RCP<const Thyra::TpetraVector<ST,LO,GO,NT> > tPtr = rcp_dynamic_cast<const Thyra::TpetraVector<ST,LO,GO,NT> >(dOpl->getDiag(),true);
1432  RCP<const Tpetra::Vector<ST,LO,GO,NT> > diagPtr = rcp_dynamic_cast<const Tpetra::Vector<ST,LO,GO,NT> >(tPtr->getConstTpetraVector(),true);
1433 
1434  // Scale by the diagonal operator
1435  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()));
1436  explicitCrsOp->leftScale(*diagPtr);
1437  explicitCrsOp->resumeFill();
1438  explicitCrsOp->scale(scalarr);
1439  explicitCrsOp->fillComplete(tCrsOpr->getDomainMap(),tCrsOpr->getRangeMap());
1440  return Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1441 
1442  } else { // Assume Epetra and we can use transformers
1443 
1444  // build implicit multiply
1445  const LinearOp implicitOp = Thyra::multiply(opl,opr);
1446 
1447  // build a scaling transformer
1448 
1449  RCP<Thyra::LinearOpTransformerBase<double> > prodTrans
1450  = Thyra::epetraExtDiagScalingTransformer();
1451 
1452  // check to see if a scaling transformer works: if not use the
1453  // DiagScaledMatrixProduct transformer
1454  if(not prodTrans->isCompatible(*implicitOp))
1455  prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
1456 
1457  // build operator destination operator
1458  ModifiableLinearOp explicitOp;
1459 
1460  // if neccessary build a operator to put the explicit multiply into
1461  if(destOp==Teuchos::null)
1462  explicitOp = prodTrans->createOutputOp();
1463  else
1464  explicitOp = destOp;
1465 
1466  // perform multiplication
1467  prodTrans->transform(*implicitOp,explicitOp.ptr());
1468 
1469  // label it
1470  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1471  " * " + opr->getObjectLabel() + " )");
1472 
1473  return explicitOp;
1474  }
1475 }
1476 
1487 const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr)
1488 {
1489  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1490  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1491 
1492  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix add
1493 
1494  // Get left and right Tpetra crs operators
1495  ST scalarl = 0.0;
1496  bool transpl = false;
1497  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1498  ST scalarr = 0.0;
1499  bool transpr = false;
1500  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1501 
1502  // Build output operator
1503  RCP<Thyra::LinearOpBase<ST> > explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1504  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1505 
1506  // Do explicit matrix-matrix add
1507  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1508  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1509  return tExplicitOp;
1510 
1511  }else{//Assume Epetra
1512  // build implicit add
1513  const LinearOp implicitOp = Thyra::add(opl,opr);
1514 
1515  // build transformer
1516  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1517  Thyra::epetraExtAddTransformer();
1518 
1519  // build operator and add
1520  const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
1521  prodTrans->transform(*implicitOp,explicitOp.ptr());
1522  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1523  " + " + opr->getObjectLabel() + " )");
1524 
1525  return explicitOp;
1526  }
1527 }
1528 
1541 const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
1542  const ModifiableLinearOp & destOp)
1543 {
1544  bool isTpetral = Teko::TpetraHelpers::isTpetraLinearOp(opl);
1545  bool isTpetrar = Teko::TpetraHelpers::isTpetraLinearOp(opr);
1546 
1547  if(isTpetral && isTpetrar){ // Both operators are Tpetra matrices so use the explicit Tpetra matrix-matrix add
1548 
1549  // Get left and right Tpetra crs operators
1550  ST scalarl = 0.0;
1551  bool transpl = false;
1552  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpl = Teko::TpetraHelpers::getTpetraCrsMatrix(opl, &scalarl, &transpl);
1553  ST scalarr = 0.0;
1554  bool transpr = false;
1555  RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > tCrsOpr = Teko::TpetraHelpers::getTpetraCrsMatrix(opr, &scalarr, &transpr);
1556 
1557  // Build output operator
1558  RCP<Thyra::LinearOpBase<ST> > explicitOp;
1559  if(destOp!=Teuchos::null)
1560  explicitOp = destOp;
1561  else
1562  explicitOp = rcp(new Thyra::TpetraLinearOp<ST,LO,GO,NT>());
1563  RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tExplicitOp = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(explicitOp);
1564 
1565  // Do explicit matrix-matrix add
1566  RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > explicitCrsOp = Tpetra::MatrixMatrix::add<ST,LO,GO,NT>(scalarl,transpl,*tCrsOpl,scalarr,transpr,*tCrsOpr);
1567  tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(explicitCrsOp->getDomainMap()),explicitCrsOp);
1568  return tExplicitOp;
1569 
1570  }else{ // Assume Epetra
1571 
1572  // build implicit add
1573  const LinearOp implicitOp = Thyra::add(opl,opr);
1574 
1575  // build transformer
1576  const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
1577  Thyra::epetraExtAddTransformer();
1578 
1579  // build or reuse destination operator
1580  RCP<Thyra::LinearOpBase<double> > explicitOp;
1581  if(destOp!=Teuchos::null)
1582  explicitOp = destOp;
1583  else
1584  explicitOp = prodTrans->createOutputOp();
1585 
1586  // perform add
1587  prodTrans->transform(*implicitOp,explicitOp.ptr());
1588  explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
1589  " + " + opr->getObjectLabel() + " )");
1590 
1591  return explicitOp;
1592  }
1593 }
1594 
1599 const ModifiableLinearOp explicitSum(const LinearOp & op,
1600  const ModifiableLinearOp & destOp)
1601 {
1602  // convert operators to Epetra_CrsMatrix
1603  const RCP<const Epetra_CrsMatrix> epetraOp =
1604  rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*op), true);
1605 
1606  if(destOp==Teuchos::null) {
1607  Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(new Epetra_CrsMatrix(*epetraOp));
1608 
1609  return Thyra::nonconstEpetraLinearOp(epetraDest);
1610  }
1611 
1612  const RCP<Epetra_CrsMatrix> epetraDest =
1613  rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp), true);
1614 
1615  EpetraExt::MatrixMatrix::Add(*epetraOp,false,1.0,*epetraDest,1.0);
1616 
1617  return destOp;
1618 }
1619 
1620 const LinearOp explicitTranspose(const LinearOp & op)
1621 {
1622  RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
1623  TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
1624  "Teko::explicitTranspose Not an Epetra_Operator");
1625  RCP<const Epetra_RowMatrix> eRowMatrixOp
1626  = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(eOp);
1627  TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
1628  "Teko::explicitTranspose Not an Epetra_RowMatrix");
1629 
1630  // we now have a delete transpose operator
1631  EpetraExt::RowMatrix_Transpose tranposeOp;
1632  Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
1633 
1634  // this copy is because of a poor implementation of the EpetraExt::Transform
1635  // implementation
1636  Teuchos::RCP<Epetra_CrsMatrix> crsMat
1637  = Teuchos::rcp(new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
1638 
1639  return Thyra::epetraLinearOp(crsMat);
1640 }
1641 
1642 const LinearOp buildDiagonal(const MultiVector & src,const std::string & lbl)
1643 {
1644  RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range());
1645  Thyra::copy(*src->col(0),dst.ptr());
1646 
1647  return Thyra::diagonal<double>(dst,lbl);
1648 }
1649 
1650 const LinearOp buildInvDiagonal(const MultiVector & src,const std::string & lbl)
1651 {
1652  const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
1653  RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range());
1654  Thyra::reciprocal<double>(*srcV,dst.ptr());
1655 
1656  return Thyra::diagonal<double>(dst,lbl);
1657 }
1658 
1660 BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvv)
1661 {
1662  Teuchos::Array<MultiVector> mvA;
1663  Teuchos::Array<VectorSpace> vsA;
1664 
1665  // build arrays of multi vectors and vector spaces
1666  std::vector<MultiVector>::const_iterator itr;
1667  for(itr=mvv.begin();itr!=mvv.end();++itr) {
1668  mvA.push_back(*itr);
1669  vsA.push_back((*itr)->range());
1670  }
1671 
1672  // construct the product vector space
1673  const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
1674  = Thyra::productVectorSpace<double>(vsA);
1675 
1676  return Thyra::defaultProductMultiVector<double>(vs,mvA);
1677 }
1678 
1689 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
1690  const std::vector<int> & indices,
1691  const VectorSpace & vs,
1692  double onValue, double offValue)
1693 
1694 {
1695  using Teuchos::RCP;
1696 
1697  // create a new vector
1698  RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
1699  Thyra::put_scalar<double>(offValue,v.ptr()); // fill it with "off" values
1700 
1701  // set on values
1702  for(std::size_t i=0;i<indices.size();i++)
1703  Thyra::set_ele<double>(indices[i],onValue,v.ptr());
1704 
1705  return v;
1706 }
1707 
1732 double computeSpectralRad(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
1733  bool isHermitian,int numBlocks,int restart,int verbosity)
1734 {
1735  typedef Thyra::LinearOpBase<double> OP;
1736  typedef Thyra::MultiVectorBase<double> MV;
1737 
1738  int startVectors = 1;
1739 
1740  // construct an initial guess
1741  const RCP<MV> ivec = Thyra::createMember(A->domain());
1742  Thyra::randomize(-1.0,1.0,ivec.ptr());
1743 
1744  RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1745  = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1746  eigProb->setNEV(1);
1747  eigProb->setHermitian(isHermitian);
1748 
1749  // set the problem up
1750  if(not eigProb->setProblem()) {
1751  // big time failure!
1752  return Teuchos::ScalarTraits<double>::nan();
1753  }
1754 
1755  // we want largert magnitude eigenvalue
1756  std::string which("LM"); // largest magnitude
1757 
1758  // Create the parameter list for the eigensolver
1759  // verbosity+=Anasazi::TimingDetails;
1760  Teuchos::ParameterList MyPL;
1761  MyPL.set( "Verbosity", verbosity );
1762  MyPL.set( "Which", which );
1763  MyPL.set( "Block Size", startVectors );
1764  MyPL.set( "Num Blocks", numBlocks );
1765  MyPL.set( "Maximum Restarts", restart );
1766  MyPL.set( "Convergence Tolerance", tol );
1767 
1768  // build status test manager
1769  // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
1770  // = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(numBlocks*(restart+1)));
1771 
1772  // Create the Block Krylov Schur solver
1773  // This takes as inputs the eigenvalue problem and the solver parameters
1774  Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1775 
1776  // Solve the eigenvalue problem, and save the return code
1777  Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1778 
1779  if(solverreturn==Anasazi::Unconverged) {
1780  double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
1781  double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
1782 
1783  return -std::abs(std::sqrt(real*real+comp*comp));
1784 
1785  // cout << "Anasazi::BlockKrylovSchur::solve() did not converge!" << std::endl;
1786  // return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
1787  }
1788  else { // solverreturn==Anasazi::Converged
1789  double real = eigProb->getSolution().Evals.begin()->realpart;
1790  double comp = eigProb->getSolution().Evals.begin()->imagpart;
1791 
1792  return std::abs(std::sqrt(real*real+comp*comp));
1793 
1794  // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
1795  // return std::abs(eigProb->getSolution().Evals.begin()->realpart);
1796  }
1797 }
1798 
1822 double computeSmallestMagEig(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
1823  bool isHermitian,int numBlocks,int restart,int verbosity)
1824 {
1825  typedef Thyra::LinearOpBase<double> OP;
1826  typedef Thyra::MultiVectorBase<double> MV;
1827 
1828  int startVectors = 1;
1829 
1830  // construct an initial guess
1831  const RCP<MV> ivec = Thyra::createMember(A->domain());
1832  Thyra::randomize(-1.0,1.0,ivec.ptr());
1833 
1834  RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
1835  = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
1836  eigProb->setNEV(1);
1837  eigProb->setHermitian(isHermitian);
1838 
1839  // set the problem up
1840  if(not eigProb->setProblem()) {
1841  // big time failure!
1842  return Teuchos::ScalarTraits<double>::nan();
1843  }
1844 
1845  // we want largert magnitude eigenvalue
1846  std::string which("SM"); // smallest magnitude
1847 
1848  // Create the parameter list for the eigensolver
1849  Teuchos::ParameterList MyPL;
1850  MyPL.set( "Verbosity", verbosity );
1851  MyPL.set( "Which", which );
1852  MyPL.set( "Block Size", startVectors );
1853  MyPL.set( "Num Blocks", numBlocks );
1854  MyPL.set( "Maximum Restarts", restart );
1855  MyPL.set( "Convergence Tolerance", tol );
1856 
1857  // build status test manager
1858  // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
1859  // = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(10));
1860 
1861  // Create the Block Krylov Schur solver
1862  // This takes as inputs the eigenvalue problem and the solver parameters
1863  Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
1864 
1865  // Solve the eigenvalue problem, and save the return code
1866  Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
1867 
1868  if(solverreturn==Anasazi::Unconverged) {
1869  // cout << "Anasazi::BlockKrylovSchur::solve() did not converge! " << std::endl;
1870  return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
1871  }
1872  else { // solverreturn==Anasazi::Converged
1873  // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
1874  return std::abs(eigProb->getSolution().Evals.begin()->realpart);
1875  }
1876 }
1877 
1886 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt)
1887 {
1888  switch(dt) {
1889  case Diagonal:
1890  return getDiagonalOp(A);
1891  case Lumped:
1892  return getLumpedMatrix(A);
1893  case AbsRowSum:
1894  return getAbsRowSumMatrix(A);
1895  case NotDiag:
1896  default:
1897  TEUCHOS_TEST_FOR_EXCEPT(true);
1898  };
1899 
1900  return Teuchos::null;
1901 }
1902 
1911 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const Teko::DiagonalType & dt)
1912 {
1913  switch(dt) {
1914  case Diagonal:
1915  return getInvDiagonalOp(A);
1916  case Lumped:
1917  return getInvLumpedMatrix(A);
1918  case AbsRowSum:
1919  return getAbsRowSumInvMatrix(A);
1920  case NotDiag:
1921  default:
1922  TEUCHOS_TEST_FOR_EXCEPT(true);
1923  };
1924 
1925  return Teuchos::null;
1926 }
1927 
1934 std::string getDiagonalName(const DiagonalType & dt)
1935 {
1936  switch(dt) {
1937  case Diagonal:
1938  return "Diagonal";
1939  case Lumped:
1940  return "Lumped";
1941  case AbsRowSum:
1942  return "AbsRowSum";
1943  case NotDiag:
1944  return "NotDiag";
1945  case BlkDiag:
1946  return "BlkDiag";
1947  };
1948 
1949  return "<error>";
1950 }
1951 
1960 DiagonalType getDiagonalType(std::string name)
1961 {
1962  if(name=="Diagonal")
1963  return Diagonal;
1964  if(name=="Lumped")
1965  return Lumped;
1966  if(name=="AbsRowSum")
1967  return AbsRowSum;
1968  if(name=="BlkDiag")
1969  return BlkDiag;
1970 
1971  return NotDiag;
1972 }
1973 
1974 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,const LinearOp & Op){
1975 #ifdef Teko_ENABLE_Isorropia
1976  Teuchos::ParameterList probeList;
1977  Prober prober(G,probeList,true);
1978  Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(new Epetra_CrsMatrix(Copy,*G));
1980  prober.probe(Mwrap,*Mat);
1981  return Thyra::epetraLinearOp(Mat);
1982 #else
1983  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,"Probe requires Isorropia");
1984 #endif
1985 }
1986 
1987 double norm_1(const MultiVector & v,std::size_t col)
1988 {
1989  Teuchos::Array<double> n(v->domain()->dim());
1990  Thyra::norms_1<double>(*v,n);
1991 
1992  return n[col];
1993 }
1994 
1995 double norm_2(const MultiVector & v,std::size_t col)
1996 {
1997  Teuchos::Array<double> n(v->domain()->dim());
1998  Thyra::norms_2<double>(*v,n);
1999 
2000  return n[col];
2001 }
2002 
2003 void putScalar(const ModifiableLinearOp & op,double scalar)
2004 {
2005  try {
2006  // get Epetra_Operator
2007  RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
2008 
2009  // cast it to a CrsMatrix
2010  RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
2011 
2012  eCrsOp->PutScalar(scalar);
2013  }
2014  catch (std::exception & e) {
2015  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
2016 
2017  *out << "Teko: putScalar requires an Epetra_CrsMatrix\n";
2018  *out << " Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
2019  *out << " OR\n";
2020  *out << " Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
2021  *out << std::endl;
2022  *out << "*** THROWN EXCEPTION ***\n";
2023  *out << e.what() << std::endl;
2024  *out << "************************\n";
2025 
2026  throw e;
2027  }
2028 }
2029 
2030 void clipLower(MultiVector & v,double lowerBound)
2031 {
2032  using Teuchos::RCP;
2033  using Teuchos::rcp_dynamic_cast;
2034 
2035  // cast so entries are accessible
2036  // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2037  // = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
2038 
2039  for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2040  RCP<Thyra::SpmdVectorBase<double> > spmdVec
2041  = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2042 
2043  Teuchos::ArrayRCP<double> values;
2044  // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2045  spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2046  for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2047  if(values[j]<lowerBound)
2048  values[j] = lowerBound;
2049  }
2050  }
2051 }
2052 
2053 void clipUpper(MultiVector & v,double upperBound)
2054 {
2055  using Teuchos::RCP;
2056  using Teuchos::rcp_dynamic_cast;
2057 
2058  // cast so entries are accessible
2059  // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2060  // = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
2061  for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2062  RCP<Thyra::SpmdVectorBase<double> > spmdVec
2063  = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2064 
2065  Teuchos::ArrayRCP<double> values;
2066  // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2067  spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2068  for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2069  if(values[j]>upperBound)
2070  values[j] = upperBound;
2071  }
2072  }
2073 }
2074 
2075 void replaceValue(MultiVector & v,double currentValue,double newValue)
2076 {
2077  using Teuchos::RCP;
2078  using Teuchos::rcp_dynamic_cast;
2079 
2080  // cast so entries are accessible
2081  // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec
2082  // = rcp_dynamic_cast<Thyra::SpmdMultiVectorBase<double> >(v,true);
2083  for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
2084  RCP<Thyra::SpmdVectorBase<double> > spmdVec
2085  = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
2086 
2087  Teuchos::ArrayRCP<double> values;
2088  // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
2089  spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
2090  for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
2091  if(values[j]==currentValue)
2092  values[j] = newValue;
2093  }
2094  }
2095 }
2096 
2097 void columnAverages(const MultiVector & v,std::vector<double> & averages)
2098 {
2099  averages.resize(v->domain()->dim());
2100 
2101  // sum over each column
2102  Thyra::sums<double>(*v,averages);
2103 
2104  // build averages
2105  Thyra::Ordinal rows = v->range()->dim();
2106  for(std::size_t i=0;i<averages.size();i++)
2107  averages[i] = averages[i]/rows;
2108 }
2109 
2110 double average(const MultiVector & v)
2111 {
2112  Thyra::Ordinal rows = v->range()->dim();
2113  Thyra::Ordinal cols = v->domain()->dim();
2114 
2115  std::vector<double> averages;
2116  columnAverages(v,averages);
2117 
2118  double sum = 0.0;
2119  for(std::size_t i=0;i<averages.size();i++)
2120  sum += averages[i] * rows;
2121 
2122  return sum/(rows*cols);
2123 }
2124 
2125 }
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 .