Teko Version of the Day
Teko_Utilities.cpp
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00047 #include "Teko_Config.h"
00048 #include "Teko_Utilities.hpp"
00049 
00050 // Thyra includes
00051 #include "Thyra_MultiVectorStdOps.hpp"
00052 #include "Thyra_ZeroLinearOpBase.hpp"
00053 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00054 #include "Thyra_DefaultAddedLinearOp.hpp"
00055 #include "Thyra_EpetraExtDiagScaledMatProdTransformer.hpp"
00056 #include "Thyra_EpetraExtDiagScalingTransformer.hpp"
00057 #include "Thyra_EpetraExtAddTransformer.hpp"
00058 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00059 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00060 #include "Thyra_DefaultZeroLinearOp.hpp"
00061 #include "Thyra_DefaultProductMultiVector.hpp"
00062 #include "Thyra_DefaultProductVectorSpace.hpp"
00063 #include "Thyra_MultiVectorStdOps.hpp"
00064 #include "Thyra_VectorStdOps.hpp"
00065 #include "Thyra_SpmdVectorBase.hpp"
00066 #include "Thyra_get_Epetra_Operator.hpp"
00067 #include "Thyra_EpetraThyraWrappers.hpp"
00068 #include "Thyra_EpetraOperatorWrapper.hpp"
00069 #include "Thyra_EpetraLinearOp.hpp"
00070 
00071 // Teuchos includes
00072 #include "Teuchos_Array.hpp"
00073 
00074 // Epetra includes
00075 #include "Epetra_Operator.h"
00076 #include "Epetra_CrsGraph.h"
00077 #include "Epetra_CrsMatrix.h"
00078 #include "Epetra_Vector.h"
00079 #include "Epetra_Map.h"
00080 
00081 #include "EpetraExt_Transpose_RowMatrix.h"
00082 #include "EpetraExt_MatrixMatrix.h"
00083 
00084 // Anasazi includes
00085 #include "AnasaziBasicEigenproblem.hpp"
00086 #include "AnasaziThyraAdapter.hpp"
00087 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00088 #include "AnasaziBlockKrylovSchur.hpp"
00089 #include "AnasaziStatusTestMaxIters.hpp"
00090 
00091 // Isorropia includes
00092 #ifdef Teko_ENABLE_Isorropia
00093 #include "Isorropia_EpetraProber.hpp"
00094 #endif
00095 
00096 // Teko includes
00097 #include "Epetra/Teko_EpetraHelpers.hpp"
00098 #include "Epetra/Teko_EpetraOperatorWrapper.hpp"
00099 
00100 #include <cmath>
00101 
00102 namespace Teko {
00103 
00104 using Teuchos::rcp;
00105 using Teuchos::rcp_dynamic_cast;
00106 using Teuchos::RCP;
00107 #ifdef Teko_ENABLE_Isorropia
00108 using Isorropia::Epetra::Prober;
00109 #endif
00110 
00111 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
00112 { 
00113    Teuchos::RCP<Teuchos::FancyOStream> os = 
00114          Teuchos::VerboseObjectBase::getDefaultOStream(); 
00115   
00116    //os->setShowProcRank(true);
00117    //os->setOutputToRootOnly(-1);
00118    return os;
00119 }
00120 
00121 // distance function...not parallel...entirely internal to this cpp file
00122 inline double dist(int dim,double * coords,int row,int col)
00123 {
00124    double value = 0.0;
00125    for(int i=0;i<dim;i++)
00126       value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
00127 
00128    // the distance between the two
00129    return std::sqrt(value);
00130 }
00131 
00132 // distance function...not parallel...entirely internal to this cpp file
00133 inline double dist(double * x,double * y,double * z,int stride,int row,int col)
00134 {
00135    double value = 0.0;
00136    if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
00137    if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
00138    if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
00139 
00140    // the distance between the two
00141    return std::sqrt(value);
00142 }
00143 
00162 RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil)
00163 {
00164    // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
00165    RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
00166                                                        stencil.MaxNumEntries()+1),true);
00167 
00168    // allocate an additional value for the diagonal, if neccessary
00169    std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
00170    std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
00171 
00172    // loop over all the rows
00173    for(int j=0;j<gl->NumMyRows();j++) {
00174       int row = gl->GRID(j);
00175       double diagValue = 0.0;
00176       int diagInd = -1;
00177       int rowSz = 0;
00178 
00179       // extract a copy of this row...put it in rowData, rowIndicies
00180       stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
00181  
00182       // loop over elements of row
00183       for(int i=0;i<rowSz;i++) {
00184          int col = rowInd[i];
00185 
00186          // is this a 0 entry masquerading as some thing else?
00187          double value = rowData[i];
00188          if(value==0) continue;
00189 
00190          // for nondiagonal entries
00191          if(row!=col) {
00192             double d = dist(dim,coords,row,col);
00193             rowData[i] = -1.0/d;
00194             diagValue += rowData[i];
00195          }
00196          else 
00197             diagInd = i;
00198       }
00199     
00200       // handle diagonal entry
00201       if(diagInd<0) { // diagonal not in row
00202          rowData[rowSz] = -diagValue;
00203          rowInd[rowSz] = row;
00204          rowSz++;
00205       }
00206       else { // diagonal in row
00207          rowData[diagInd] = -diagValue;
00208          rowInd[diagInd] = row;
00209       }
00210 
00211       // insert row data into graph Laplacian matrix
00212       TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
00213    }
00214 
00215    gl->FillComplete();
00216 
00217    return gl;
00218 }
00219 
00242 RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil)
00243 {
00244    // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
00245    RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
00246                                                        stencil.MaxNumEntries()+1),true);
00247 
00248    // allocate an additional value for the diagonal, if neccessary
00249    std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
00250    std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
00251 
00252    // loop over all the rows
00253    for(int j=0;j<gl->NumMyRows();j++) {
00254       int row = gl->GRID(j);
00255       double diagValue = 0.0;
00256       int diagInd = -1;
00257       int rowSz = 0;
00258 
00259       // extract a copy of this row...put it in rowData, rowIndicies
00260       stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
00261  
00262       // loop over elements of row
00263       for(int i=0;i<rowSz;i++) {
00264          int col = rowInd[i];
00265 
00266          // is this a 0 entry masquerading as some thing else?
00267          double value = rowData[i];
00268          if(value==0) continue;
00269 
00270          // for nondiagonal entries
00271          if(row!=col) {
00272             double d = dist(x,y,z,stride,row,col);
00273             rowData[i] = -1.0/d;
00274             diagValue += rowData[i];
00275          }
00276          else 
00277             diagInd = i;
00278       }
00279     
00280       // handle diagonal entry
00281       if(diagInd<0) { // diagonal not in row
00282          rowData[rowSz] = -diagValue;
00283          rowInd[rowSz] = row;
00284          rowSz++;
00285       }
00286       else { // diagonal in row
00287          rowData[diagInd] = -diagValue;
00288          rowInd[diagInd] = row;
00289       }
00290 
00291       // insert row data into graph Laplacian matrix
00292       TEUCHOS_TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
00293    }
00294 
00295    gl->FillComplete();
00296 
00297    return gl;
00298 }
00299 
00315 void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
00316 {
00317    // Thyra::apply(*A,Thyra::NONCONJ_ELE,*x,&*y,alpha,beta);
00318    Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
00319 }
00320 
00323 void update(double alpha,const MultiVector & x,double beta,MultiVector & y)
00324 {
00325    Teuchos::Array<double> scale;
00326    Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > >  vec;
00327 
00328    // build arrays needed for linear combo
00329    scale.push_back(alpha);
00330    vec.push_back(x.ptr());
00331 
00332    // compute linear combination
00333    Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
00334 }
00335 
00337 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
00338 {
00339    int rows = blockRowCount(blo);
00340 
00341    TEUCHOS_ASSERT(rows==blockColCount(blo));
00342 
00343    RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
00344    RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
00345 
00346    // allocate new operator
00347    BlockedLinearOp upper = createBlockedOp();
00348  
00349    // build new operator 
00350    upper->beginBlockFill(rows,rows);
00351 
00352    for(int i=0;i<rows;i++) {
00353       // put zero operators on the diagonal
00354       // this gurantees the vector space of
00355       // the new operator are fully defined
00356       RCP<const Thyra::LinearOpBase<double> > zed 
00357             = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
00358       upper->setBlock(i,i,zed);
00359 
00360       for(int j=i+1;j<rows;j++) {
00361          // get block i,j
00362          LinearOp uij = blo->getBlock(i,j);
00363 
00364          // stuff it in U
00365          if(uij!=Teuchos::null)
00366             upper->setBlock(i,j,uij);
00367       }
00368    }
00369    if(callEndBlockFill)
00370       upper->endBlockFill();
00371 
00372    return upper;
00373 }
00374 
00376 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
00377 {
00378    int rows = blockRowCount(blo);
00379 
00380    TEUCHOS_ASSERT(rows==blockColCount(blo));
00381 
00382    RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
00383    RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
00384 
00385    // allocate new operator
00386    BlockedLinearOp lower = createBlockedOp();
00387  
00388    // build new operator 
00389    lower->beginBlockFill(rows,rows);
00390 
00391    for(int i=0;i<rows;i++) {
00392       // put zero operators on the diagonal
00393       // this gurantees the vector space of
00394       // the new operator are fully defined
00395       RCP<const Thyra::LinearOpBase<double> > zed 
00396             = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
00397       lower->setBlock(i,i,zed);
00398 
00399       for(int j=0;j<i;j++) {
00400          // get block i,j
00401          LinearOp uij = blo->getBlock(i,j);
00402 
00403          // stuff it in U
00404          if(uij!=Teuchos::null)
00405             lower->setBlock(i,j,uij);
00406       }
00407    }
00408    if(callEndBlockFill)
00409       lower->endBlockFill();
00410 
00411    return lower;
00412 }
00413 
00433 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo)
00434 {
00435    int rows = blockRowCount(blo);
00436 
00437    TEUCHOS_ASSERT(rows==blockColCount(blo)); // assert that matrix is square
00438 
00439    RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
00440    RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
00441 
00442    // allocate new operator
00443    BlockedLinearOp zeroOp = createBlockedOp();
00444  
00445    // build new operator 
00446    zeroOp->beginBlockFill(rows,rows);
00447 
00448    for(int i=0;i<rows;i++) {
00449       // put zero operators on the diagonal
00450       // this gurantees the vector space of
00451       // the new operator are fully defined
00452       RCP<const Thyra::LinearOpBase<double> > zed 
00453             = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
00454       zeroOp->setBlock(i,i,zed);
00455    }
00456 
00457    return zeroOp;
00458 }
00459 
00461 bool isZeroOp(const LinearOp op)
00462 {
00463    // if operator is null...then its zero!
00464    if(op==Teuchos::null) return true;
00465 
00466    // try to cast it to a zero linear operator
00467    const LinearOp test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(op);
00468 
00469    // if it works...then its zero...otherwise its null
00470    return test!=Teuchos::null;
00471 }
00472 
00481 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op)
00482 {
00483    RCP<const Epetra_CrsMatrix> eCrsOp;
00484 
00485    try {
00486       // get Epetra_Operator
00487       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00488 
00489       // cast it to a CrsMatrix
00490       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00491    }
00492    catch (std::exception & e) {
00493       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00494 
00495       *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix\n";
00496       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00497       *out << "           OR\n";
00498       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00499       *out << std::endl;
00500       *out << "*** THROWN EXCEPTION ***\n";
00501       *out << e.what() << std::endl;
00502       *out << "************************\n";
00503       
00504       throw e;
00505    }
00506 
00507    // extract diagonal
00508    const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00509    Epetra_Vector & diag = *ptrDiag;
00510    
00511    // compute absolute value row sum
00512    diag.PutScalar(0.0);
00513    for(int i=0;i<eCrsOp->NumMyRows();i++) {
00514       double * values = 0;
00515       int numEntries;
00516       eCrsOp->ExtractMyRowView(i,numEntries,values);
00517 
00518       // build abs value row sum
00519       for(int j=0;j<numEntries;j++)
00520          diag[i] += std::abs(values[j]);
00521    }
00522 
00523    // build Thyra diagonal operator
00524    return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
00525 }
00526 
00535 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op)
00536 {
00537    RCP<const Epetra_CrsMatrix> eCrsOp;
00538 
00539    try {
00540       // get Epetra_Operator
00541       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00542 
00543       // cast it to a CrsMatrix
00544       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00545    }
00546    catch (std::exception & e) {
00547       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00548 
00549       *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix\n";
00550       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00551       *out << "           OR\n";
00552       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00553       *out << std::endl;
00554       *out << "*** THROWN EXCEPTION ***\n";
00555       *out << e.what() << std::endl;
00556       *out << "************************\n";
00557       
00558       throw e;
00559    }
00560 
00561    // extract diagonal
00562    const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00563    Epetra_Vector & diag = *ptrDiag;
00564    
00565    // compute absolute value row sum
00566    diag.PutScalar(0.0);
00567    for(int i=0;i<eCrsOp->NumMyRows();i++) {
00568       double * values = 0;
00569       int numEntries;
00570       eCrsOp->ExtractMyRowView(i,numEntries,values);
00571 
00572       // build abs value row sum
00573       for(int j=0;j<numEntries;j++)
00574          diag[i] += std::abs(values[j]);
00575    }
00576    diag.Reciprocal(diag); // invert entries
00577 
00578    // build Thyra diagonal operator
00579    return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSumInv( " + op->getObjectLabel() + " )");
00580 }
00581 
00589 ModifiableLinearOp getLumpedMatrix(const LinearOp & op)
00590 {
00591    RCP<Thyra::VectorBase<double> > ones = Thyra::createMember(op->domain());
00592    RCP<Thyra::VectorBase<double> > diag = Thyra::createMember(op->range());
00593 
00594    // set to all ones
00595    Thyra::assign(ones.ptr(),1.0);
00596 
00597    // compute lumped diagonal
00598    // Thyra::apply(*op,Thyra::NONCONJ_ELE,*ones,&*diag);
00599    Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
00600 
00601    return rcp(new Thyra::DefaultDiagonalLinearOp<double>(diag));
00602 }
00603 
00612 ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op)
00613 {
00614    RCP<Thyra::VectorBase<double> > ones = Thyra::createMember(op->domain());
00615    RCP<Thyra::VectorBase<double> > diag = Thyra::createMember(op->range());
00616 
00617    // set to all ones
00618    Thyra::assign(ones.ptr(),1.0);
00619 
00620    // compute lumped diagonal
00621    Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
00622    Thyra::reciprocal(*diag,diag.ptr());
00623 
00624    return rcp(new Thyra::DefaultDiagonalLinearOp<double>(diag));
00625 }
00626 
00638 const ModifiableLinearOp getDiagonalOp(const LinearOp & op)
00639 {
00640    RCP<const Epetra_CrsMatrix> eCrsOp;
00641 
00642    try {
00643       // get Epetra_Operator
00644       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00645 
00646       // cast it to a CrsMatrix
00647       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00648    }
00649    catch (std::exception & e) {
00650       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00651 
00652       *out << "Teko: getDiagonalOp requires an Epetra_CrsMatrix\n";
00653       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00654       *out << "           OR\n";
00655       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00656       *out << std::endl;
00657       *out << "*** THROWN EXCEPTION ***\n";
00658       *out << e.what() << std::endl;
00659       *out << "************************\n";
00660       
00661       throw e;
00662    }
00663 
00664    // extract diagonal
00665    const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00666    TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
00667 
00668    // build Thyra diagonal operator
00669    return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"diag( " + op->getObjectLabel() + " )");
00670 }
00671 
00672 const MultiVector getDiagonal(const LinearOp & op)
00673 {
00674    RCP<const Epetra_CrsMatrix> eCrsOp;
00675 
00676    try {
00677       // get Epetra_Operator
00678       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00679 
00680       // cast it to a CrsMatrix
00681       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00682    }
00683    catch (std::exception & e) {
00684       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00685 
00686       *out << "Teko: getDiagonal requires an Epetra_CrsMatrix\n";
00687       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00688       *out << "           OR\n";
00689       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00690       *out << std::endl;
00691       *out << "*** THROWN EXCEPTION ***\n";
00692       *out << e.what() << std::endl;
00693       *out << "************************\n";
00694       
00695       throw e;
00696    }
00697 
00698    // extract diagonal
00699    const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00700    TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
00701 
00702    return Thyra::create_Vector(diag,Thyra::create_VectorSpace(Teuchos::rcpFromRef(eCrsOp->RowMap())));
00703 }
00704 
00705 const MultiVector getDiagonal(const Teko::LinearOp & A,const DiagonalType & dt)
00706 {
00707    LinearOp diagOp = Teko::getDiagonalOp(A,dt);
00708 
00709    Teuchos::RCP<const Thyra::MultiVectorBase<double> > v = 
00710          Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(diagOp)->getDiag(); 
00711    return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
00712 }
00713 
00725 const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op)
00726 {
00727    RCP<const Epetra_CrsMatrix> eCrsOp;
00728 
00729    try {
00730       // get Epetra_Operator
00731       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00732 
00733       // cast it to a CrsMatrix
00734       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00735    }
00736    catch (std::exception & e) {
00737       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00738 
00739       *out << "Teko: getInvDiagonalOp requires an Epetra_CrsMatrix\n";
00740       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00741       *out << "           OR\n";
00742       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00743       *out << std::endl;
00744       *out << "*** THROWN EXCEPTION ***\n";
00745       *out << e.what() << std::endl;
00746       *out << "************************\n";
00747       
00748       throw e;
00749    }
00750 
00751    // extract diagonal
00752    const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00753    TEUCHOS_TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
00754    diag->Reciprocal(*diag);
00755 
00756    // build Thyra diagonal operator
00757    return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
00758 }
00759 
00772 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr)
00773 {
00774    // build implicit multiply
00775    const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
00776 
00777    // build transformer
00778    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00779        Thyra::epetraExtDiagScaledMatProdTransformer();
00780 
00781    // build operator and multiply
00782    const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
00783    prodTrans->transform(*implicitOp,explicitOp.ptr());
00784    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00785                                      " * " + opm->getObjectLabel() +
00786                                      " * " + opr->getObjectLabel() + " )");
00787 
00788    return explicitOp;
00789 }
00790 
00805 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
00806                                           const ModifiableLinearOp & destOp)
00807 {
00808    // build implicit multiply
00809    const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
00810 
00811    // build transformer
00812    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00813        Thyra::epetraExtDiagScaledMatProdTransformer();
00814 
00815    // build operator destination operator
00816    ModifiableLinearOp explicitOp; 
00817 
00818    // if neccessary build a operator to put the explicit multiply into
00819    if(destOp==Teuchos::null)
00820       explicitOp = prodTrans->createOutputOp();
00821    else
00822       explicitOp = destOp;
00823 
00824    // perform multiplication
00825    prodTrans->transform(*implicitOp,explicitOp.ptr());
00826 
00827    // label it
00828    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00829                                      " * " + opm->getObjectLabel() +
00830                                      " * " + opr->getObjectLabel() + " )");
00831 
00832    return explicitOp;
00833 }
00834 
00845 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr)
00846 {
00847    // build implicit multiply
00848    const LinearOp implicitOp = Thyra::multiply(opl,opr);
00849 
00850    // build a scaling transformer
00851    RCP<Thyra::LinearOpTransformerBase<double> > prodTrans  
00852          = Thyra::epetraExtDiagScalingTransformer();
00853 
00854    // check to see if a scaling transformer works: if not use the 
00855    // DiagScaledMatrixProduct transformer
00856    if(not prodTrans->isCompatible(*implicitOp))
00857        prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
00858 
00859    // build operator and multiply
00860    const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
00861    prodTrans->transform(*implicitOp,explicitOp.ptr());
00862    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00863                                      " * " + opr->getObjectLabel() + " )");
00864 
00865    return explicitOp;
00866 }
00867 
00881 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
00882                                           const ModifiableLinearOp & destOp)
00883 {
00884    // build implicit multiply
00885    const LinearOp implicitOp = Thyra::multiply(opl,opr);
00886 
00887    // build a scaling transformer
00888    RCP<Thyra::LinearOpTransformerBase<double> > prodTrans  
00889          = Thyra::epetraExtDiagScalingTransformer();
00890 
00891    // check to see if a scaling transformer works: if not use the 
00892    // DiagScaledMatrixProduct transformer
00893    if(not prodTrans->isCompatible(*implicitOp))
00894        prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
00895 
00896    // build operator destination operator
00897    ModifiableLinearOp explicitOp; 
00898 
00899    // if neccessary build a operator to put the explicit multiply into
00900    if(destOp==Teuchos::null)
00901       explicitOp = prodTrans->createOutputOp();
00902    else
00903       explicitOp = destOp;
00904 
00905    // perform multiplication
00906    prodTrans->transform(*implicitOp,explicitOp.ptr());
00907 
00908    // label it
00909    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00910                                      " * " + opr->getObjectLabel() + " )");
00911 
00912    return explicitOp;
00913 }
00914 
00925 const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr)
00926 {
00927    // build implicit multiply
00928    const LinearOp implicitOp = Thyra::add(opl,opr);
00929  
00930    // build transformer
00931    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00932        Thyra::epetraExtAddTransformer();
00933 
00934    // build operator and multiply
00935    const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
00936    prodTrans->transform(*implicitOp,explicitOp.ptr());
00937    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00938                                      " + " + opr->getObjectLabel() + " )");
00939 
00940    return explicitOp;
00941 }
00942 
00955 const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
00956                                      const ModifiableLinearOp & destOp)
00957 {
00958    // build implicit multiply
00959    const LinearOp implicitOp = Thyra::add(opl,opr);
00960  
00961    // build transformer
00962    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00963        Thyra::epetraExtAddTransformer();
00964 
00965    // build or reuse destination operator
00966    RCP<Thyra::LinearOpBase<double> > explicitOp;
00967    if(destOp!=Teuchos::null)
00968       explicitOp = destOp;
00969    else
00970       explicitOp = prodTrans->createOutputOp();
00971 
00972    // perform multiplication
00973    prodTrans->transform(*implicitOp,explicitOp.ptr());
00974    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00975                                      " + " + opr->getObjectLabel() + " )");
00976 
00977    return explicitOp;
00978 }
00979 
00984 const ModifiableLinearOp explicitSum(const LinearOp & op,
00985                                      const ModifiableLinearOp & destOp)
00986 {
00987    // convert operators to Epetra_CrsMatrix
00988    const RCP<const Epetra_CrsMatrix> epetraOp =
00989          rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*op), true);
00990 
00991    if(destOp==Teuchos::null) {
00992       Teuchos::RCP<Epetra_Operator> epetraDest = Teuchos::rcp(new Epetra_CrsMatrix(*epetraOp));
00993 
00994       return Thyra::nonconstEpetraLinearOp(epetraDest);
00995    }
00996 
00997    const RCP<Epetra_CrsMatrix> epetraDest =
00998          rcp_dynamic_cast<Epetra_CrsMatrix>(get_Epetra_Operator(*destOp), true);
00999 
01000    EpetraExt::MatrixMatrix::Add(*epetraOp,false,1.0,*epetraDest,1.0);
01001 
01002    return destOp;
01003 }
01004 
01005 const LinearOp explicitTranspose(const LinearOp & op)
01006 {
01007    RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
01008    TEUCHOS_TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
01009                              "Teko::explicitTranspose Not an Epetra_Operator");
01010    RCP<const Epetra_RowMatrix> eRowMatrixOp 
01011          = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(eOp);
01012    TEUCHOS_TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
01013                              "Teko::explicitTranspose Not an Epetra_RowMatrix");
01014 
01015    // we now have a delete transpose operator
01016    EpetraExt::RowMatrix_Transpose tranposeOp;
01017    Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
01018 
01019    // this copy is because of a poor implementation of the EpetraExt::Transform 
01020    // implementation
01021    Teuchos::RCP<Epetra_CrsMatrix> crsMat 
01022          = Teuchos::rcp(new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
01023 
01024    return Thyra::epetraLinearOp(crsMat);
01025 }
01026 
01027 const LinearOp buildDiagonal(const MultiVector & src,const std::string & lbl)
01028 {
01029    RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range()); 
01030    Thyra::copy(*src->col(0),dst.ptr());
01031    
01032    return Thyra::diagonal<double>(dst,lbl);
01033 }
01034 
01035 const LinearOp buildInvDiagonal(const MultiVector & src,const std::string & lbl)
01036 {
01037    const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
01038    RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range()); 
01039    Thyra::reciprocal<double>(*srcV,dst.ptr());
01040 
01041    return Thyra::diagonal<double>(dst,lbl);
01042 }
01043 
01045 BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvv)
01046 {
01047    Teuchos::Array<MultiVector> mvA;
01048    Teuchos::Array<VectorSpace> vsA;
01049 
01050    // build arrays of multi vectors and vector spaces
01051    std::vector<MultiVector>::const_iterator itr;
01052    for(itr=mvv.begin();itr!=mvv.end();++itr) {
01053       mvA.push_back(*itr);
01054       vsA.push_back((*itr)->range());
01055    }
01056 
01057    // construct the product vector space
01058    const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
01059          = Thyra::productVectorSpace<double>(vsA);
01060 
01061    return Thyra::defaultProductMultiVector<double>(vs,mvA);
01062 }
01063 
01074 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
01075                                  const std::vector<int> & indices,
01076                                  const VectorSpace & vs,
01077                                  double onValue, double offValue)
01078 
01079 {
01080    using Teuchos::RCP;
01081  
01082    // create a new vector
01083    RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
01084    Thyra::put_scalar<double>(offValue,v.ptr()); // fill it with "off" values
01085 
01086    // set on values
01087    for(std::size_t i=0;i<indices.size();i++) 
01088       Thyra::set_ele<double>(indices[i],onValue,v.ptr());
01089 
01090    return v;
01091 }
01092 
01117 double computeSpectralRad(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
01118                           bool isHermitian,int numBlocks,int restart,int verbosity)
01119 {
01120    typedef Thyra::LinearOpBase<double> OP;
01121    typedef Thyra::MultiVectorBase<double> MV;
01122 
01123    int startVectors = 1;
01124 
01125    // construct an initial guess
01126    const RCP<MV> ivec = Thyra::createMember(A->domain());
01127    Thyra::randomize(-1.0,1.0,ivec.ptr());
01128    
01129    RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
01130          = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
01131    eigProb->setNEV(1);
01132    eigProb->setHermitian(isHermitian);
01133 
01134    // set the problem up
01135    if(not eigProb->setProblem()) {
01136       // big time failure!
01137       return Teuchos::ScalarTraits<double>::nan();
01138    }
01139 
01140    // we want largert magnitude eigenvalue
01141    std::string which("LM"); // largest magnitude
01142 
01143    // Create the parameter list for the eigensolver
01144    // verbosity+=Anasazi::TimingDetails;
01145    Teuchos::ParameterList MyPL;
01146    MyPL.set( "Verbosity", verbosity );
01147    MyPL.set( "Which", which );
01148    MyPL.set( "Block Size", startVectors );
01149    MyPL.set( "Num Blocks", numBlocks );
01150    MyPL.set( "Maximum Restarts", restart ); 
01151    MyPL.set( "Convergence Tolerance", tol );
01152 
01153    // build status test manager
01154    // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
01155    //       = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(numBlocks*(restart+1)));
01156 
01157    // Create the Block Krylov Schur solver
01158    // This takes as inputs the eigenvalue problem and the solver parameters
01159    Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
01160  
01161    // Solve the eigenvalue problem, and save the return code
01162    Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
01163 
01164    if(solverreturn==Anasazi::Unconverged) {
01165       double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
01166       double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
01167 
01168       return -std::abs(std::sqrt(real*real+comp*comp));
01169 
01170       // cout << "Anasazi::BlockKrylovSchur::solve() did not converge!" << std::endl; 
01171       // return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
01172    }
01173    else { // solverreturn==Anasazi::Converged
01174       double real = eigProb->getSolution().Evals.begin()->realpart;
01175       double comp = eigProb->getSolution().Evals.begin()->imagpart;
01176 
01177       return std::abs(std::sqrt(real*real+comp*comp));
01178 
01179       // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
01180       // return std::abs(eigProb->getSolution().Evals.begin()->realpart);
01181    }
01182 }
01183 
01207 double computeSmallestMagEig(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
01208                              bool isHermitian,int numBlocks,int restart,int verbosity)
01209 {
01210    typedef Thyra::LinearOpBase<double> OP;
01211    typedef Thyra::MultiVectorBase<double> MV;
01212 
01213    int startVectors = 1;
01214 
01215    // construct an initial guess
01216    const RCP<MV> ivec = Thyra::createMember(A->domain());
01217    Thyra::randomize(-1.0,1.0,ivec.ptr());
01218    
01219    RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
01220          = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
01221    eigProb->setNEV(1);
01222    eigProb->setHermitian(isHermitian);
01223 
01224    // set the problem up
01225    if(not eigProb->setProblem()) {
01226       // big time failure!
01227       return Teuchos::ScalarTraits<double>::nan();
01228    }
01229 
01230    // we want largert magnitude eigenvalue
01231    std::string which("SM"); // smallest magnitude
01232 
01233    // Create the parameter list for the eigensolver
01234    Teuchos::ParameterList MyPL;
01235    MyPL.set( "Verbosity", verbosity );
01236    MyPL.set( "Which", which );
01237    MyPL.set( "Block Size", startVectors );
01238    MyPL.set( "Num Blocks", numBlocks );
01239    MyPL.set( "Maximum Restarts", restart ); 
01240    MyPL.set( "Convergence Tolerance", tol );
01241 
01242    // build status test manager
01243    // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
01244    //       = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(10));
01245 
01246    // Create the Block Krylov Schur solver
01247    // This takes as inputs the eigenvalue problem and the solver parameters
01248    Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
01249  
01250    // Solve the eigenvalue problem, and save the return code
01251    Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
01252 
01253    if(solverreturn==Anasazi::Unconverged) {
01254       // cout << "Anasazi::BlockKrylovSchur::solve() did not converge! " << std::endl;
01255       return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
01256    }
01257    else { // solverreturn==Anasazi::Converged
01258       // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
01259       return std::abs(eigProb->getSolution().Evals.begin()->realpart);
01260    }
01261 }
01262 
01271 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt)
01272 {
01273    switch(dt) {
01274    case Diagonal:
01275       return getDiagonalOp(A);  
01276    case Lumped:
01277       return getLumpedMatrix(A);  
01278    case AbsRowSum:
01279       return getAbsRowSumMatrix(A);  
01280    case NotDiag:
01281    default:
01282       TEUCHOS_TEST_FOR_EXCEPT(true);
01283    };
01284 
01285    return Teuchos::null;
01286 }
01287 
01296 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const Teko::DiagonalType & dt)
01297 {
01298    switch(dt) {
01299    case Diagonal:
01300       return getInvDiagonalOp(A);  
01301    case Lumped:
01302       return getInvLumpedMatrix(A);  
01303    case AbsRowSum:
01304       return getAbsRowSumInvMatrix(A);  
01305    case NotDiag:
01306    default:
01307       TEUCHOS_TEST_FOR_EXCEPT(true);
01308    };
01309 
01310    return Teuchos::null;
01311 }
01312 
01319 std::string getDiagonalName(const DiagonalType & dt)
01320 {
01321    switch(dt) {
01322    case Diagonal:
01323       return "Diagonal";
01324    case Lumped:
01325       return "Lumped";
01326    case AbsRowSum:
01327       return "AbsRowSum";
01328    case NotDiag:
01329       return "NotDiag";
01330    case BlkDiag:
01331       return "BlkDiag";
01332    };
01333 
01334    return "<error>";
01335 }
01336 
01345 DiagonalType getDiagonalType(std::string name)
01346 {
01347    if(name=="Diagonal")
01348       return Diagonal;
01349    if(name=="Lumped")
01350      return Lumped;
01351    if(name=="AbsRowSum")
01352       return AbsRowSum;
01353    if(name=="BlkDiag")
01354       return BlkDiag;
01355  
01356    return NotDiag;
01357 }
01358 
01359 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,const LinearOp & Op){
01360 #ifdef Teko_ENABLE_Isorropia
01361   Teuchos::ParameterList probeList;
01362   Prober prober(G,probeList,true);
01363   Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(new Epetra_CrsMatrix(Copy,*G));
01364   Teko::Epetra::EpetraOperatorWrapper Mwrap(Op);
01365   prober.probe(Mwrap,*Mat);
01366   return Thyra::epetraLinearOp(Mat);    
01367 #else
01368   TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,"Probe requires Isorropia");
01369 #endif
01370 }
01371 
01372 double norm_1(const MultiVector & v,std::size_t col)
01373 {
01374    Teuchos::Array<double> n(v->domain()->dim());
01375    Thyra::norms_1<double>(*v,n);
01376 
01377    return n[col];
01378 }
01379 
01380 double norm_2(const MultiVector & v,std::size_t col)
01381 {
01382    Teuchos::Array<double> n(v->domain()->dim());
01383    Thyra::norms_2<double>(*v,n);
01384 
01385    return n[col];
01386 }
01387 
01388 void putScalar(const ModifiableLinearOp & op,double scalar) 
01389 {
01390    try {
01391       // get Epetra_Operator
01392       RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
01393 
01394       // cast it to a CrsMatrix
01395       RCP<Epetra_CrsMatrix> eCrsOp = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
01396 
01397      eCrsOp->PutScalar(scalar);
01398    }
01399    catch (std::exception & e) {
01400       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
01401 
01402       *out << "Teko: putScalar requires an Epetra_CrsMatrix\n";
01403       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
01404       *out << "           OR\n";
01405       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
01406       *out << std::endl;
01407       *out << "*** THROWN EXCEPTION ***\n";
01408       *out << e.what() << std::endl;
01409       *out << "************************\n";
01410       
01411       throw e;
01412    }
01413 }
01414 
01415 void clipLower(MultiVector & v,double lowerBound)
01416 {
01417    using Teuchos::RCP;
01418    using Teuchos::rcp_dynamic_cast;
01419 
01420    // cast so entries are accessible
01421    // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec  
01422    //    = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
01423    
01424    for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
01425       RCP<Thyra::SpmdVectorBase<double> > spmdVec  
01426          = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
01427    
01428       Teuchos::ArrayRCP<double> values;
01429       // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
01430       spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
01431       for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
01432          if(values[j]<lowerBound)
01433             values[j] = lowerBound;
01434       }
01435    }
01436 }
01437 
01438 void clipUpper(MultiVector & v,double upperBound)
01439 {
01440    using Teuchos::RCP;
01441    using Teuchos::rcp_dynamic_cast;
01442 
01443    // cast so entries are accessible
01444    // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec  
01445    //   = rcp_dynamic_cast<Thyra::DefaultSpmdMultiVector<double> >(v);
01446    for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
01447       RCP<Thyra::SpmdVectorBase<double> > spmdVec  
01448          = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
01449    
01450       Teuchos::ArrayRCP<double> values;
01451       // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
01452       spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
01453       for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
01454          if(values[j]>upperBound)
01455             values[j] = upperBound;
01456       }
01457    }
01458 }
01459 
01460 void replaceValue(MultiVector & v,double currentValue,double newValue)
01461 {
01462    using Teuchos::RCP;
01463    using Teuchos::rcp_dynamic_cast;
01464 
01465    // cast so entries are accessible
01466    // RCP<Thyra::SpmdMultiVectorBase<double> > spmdMVec  
01467    //    = rcp_dynamic_cast<Thyra::SpmdMultiVectorBase<double> >(v,true);
01468    for(Thyra::Ordinal i=0;i<v->domain()->dim();i++) {
01469       RCP<Thyra::SpmdVectorBase<double> > spmdVec  
01470          = rcp_dynamic_cast<Thyra::SpmdVectorBase<double> >(v->col(i),true);
01471    
01472       Teuchos::ArrayRCP<double> values;
01473       // spmdMVec->getNonconstLocalData(Teuchos::ptrFromRef(values),Teuchos::ptrFromRef(i));
01474       spmdVec->getNonconstLocalData(Teuchos::ptrFromRef(values));
01475       for(Teuchos::ArrayRCP<double>::size_type j=0;j<values.size();j++) {
01476          if(values[j]==currentValue)
01477             values[j] = newValue;
01478       }
01479    }
01480 }
01481 
01482 void columnAverages(const MultiVector & v,std::vector<double> & averages)
01483 {
01484    averages.resize(v->domain()->dim());
01485 
01486    // sum over each column
01487    Thyra::sums<double>(*v,averages);
01488 
01489    // build averages
01490    Thyra::Ordinal rows = v->range()->dim();
01491    for(std::size_t i=0;i<averages.size();i++)
01492       averages[i] = averages[i]/rows;
01493 }
01494 
01495 double average(const MultiVector & v)
01496 {
01497    Thyra::Ordinal rows = v->range()->dim();
01498    Thyra::Ordinal cols = v->domain()->dim();
01499 
01500    std::vector<double> averages;
01501    columnAverages(v,averages);
01502 
01503    double sum = 0.0;
01504    for(std::size_t i=0;i<averages.size();i++)
01505       sum += averages[i] * rows;
01506 
01507    return sum/(rows*cols);
01508 }
01509 
01510 }
 All Classes Files Functions Variables