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