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_Vector.h"
00077 #include "Epetra_Map.h"
00078 
00079 #include "EpetraExt_Transpose_RowMatrix.h"
00080 
00081 // Anasazi includes
00082 #include "AnasaziBasicEigenproblem.hpp"
00083 #include "AnasaziThyraAdapter.hpp"
00084 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00085 #include "AnasaziBlockKrylovSchur.hpp"
00086 #include "AnasaziStatusTestMaxIters.hpp"
00087 
00088 // Isorropia includes
00089 #ifdef Teko_ENABLE_Isorropia
00090 #include "Isorropia_EpetraProber.hpp"
00091 #endif
00092 
00093 // Teko includes
00094 #include "Epetra/Teko_EpetraHelpers.hpp"
00095 #include "Epetra/Teko_EpetraOperatorWrapper.hpp"
00096 
00097 #include <cmath>
00098 
00099 namespace Teko {
00100 
00101 using Teuchos::rcp;
00102 using Teuchos::rcp_dynamic_cast;
00103 using Teuchos::RCP;
00104 #ifdef Teko_ENABLE_Isorropia
00105 using Isorropia::Epetra::Prober;
00106 #endif
00107 
00108 const Teuchos::RCP<Teuchos::FancyOStream> getOutputStream()
00109 { 
00110    Teuchos::RCP<Teuchos::FancyOStream> os = 
00111          Teuchos::VerboseObjectBase::getDefaultOStream(); 
00112   
00113    //os->setShowProcRank(true);
00114    //os->setOutputToRootOnly(-1);
00115    return os;
00116 }
00117 
00118 // distance function...not parallel...entirely internal to this cpp file
00119 inline double dist(int dim,double * coords,int row,int col)
00120 {
00121    double value = 0.0;
00122    for(int i=0;i<dim;i++)
00123       value += std::pow(coords[dim*row+i]-coords[dim*col+i],2.0);
00124 
00125    // the distance between the two
00126    return std::sqrt(value);
00127 }
00128 
00129 // distance function...not parallel...entirely internal to this cpp file
00130 inline double dist(double * x,double * y,double * z,int stride,int row,int col)
00131 {
00132    double value = 0.0;
00133    if(x!=0) value += std::pow(x[stride*row]-x[stride*col],2.0);
00134    if(y!=0) value += std::pow(y[stride*row]-y[stride*col],2.0);
00135    if(z!=0) value += std::pow(z[stride*row]-z[stride*col],2.0);
00136 
00137    // the distance between the two
00138    return std::sqrt(value);
00139 }
00140 
00159 RCP<Epetra_CrsMatrix> buildGraphLaplacian(int dim,double * coords,const Epetra_CrsMatrix & stencil)
00160 {
00161    // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
00162    RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
00163                                                        stencil.MaxNumEntries()+1),true);
00164 
00165    // allocate an additional value for the diagonal, if neccessary
00166    std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
00167    std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
00168 
00169    // loop over all the rows
00170    for(int j=0;j<gl->NumMyRows();j++) {
00171       int row = gl->GRID(j);
00172       double diagValue = 0.0;
00173       int diagInd = -1;
00174       int rowSz = 0;
00175 
00176       // extract a copy of this row...put it in rowData, rowIndicies
00177       stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
00178  
00179       // loop over elements of row
00180       for(int i=0;i<rowSz;i++) {
00181          int col = rowInd[i];
00182 
00183          // is this a 0 entry masquerading as some thing else?
00184          double value = rowData[i];
00185          if(value==0) continue;
00186 
00187          // for nondiagonal entries
00188          if(row!=col) {
00189             double d = dist(dim,coords,row,col);
00190             rowData[i] = -1.0/d;
00191             diagValue += rowData[i];
00192          }
00193          else 
00194             diagInd = i;
00195       }
00196     
00197       // handle diagonal entry
00198       if(diagInd<0) { // diagonal not in row
00199          rowData[rowSz] = -diagValue;
00200          rowInd[rowSz] = row;
00201          rowSz++;
00202       }
00203       else { // diagonal in row
00204          rowData[diagInd] = -diagValue;
00205          rowInd[diagInd] = row;
00206       }
00207 
00208       // insert row data into graph Laplacian matrix
00209       TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
00210    }
00211 
00212    gl->FillComplete();
00213 
00214    return gl;
00215 }
00216 
00239 RCP<Epetra_CrsMatrix> buildGraphLaplacian(double * x,double * y,double * z,int stride,const Epetra_CrsMatrix & stencil)
00240 {
00241    // allocate a new matrix with storage for the laplacian...in case of diagonals add one extra storage
00242    RCP<Epetra_CrsMatrix> gl = rcp(new Epetra_CrsMatrix(Copy,stencil.RowMap(),stencil.ColMap(),
00243                                                        stencil.MaxNumEntries()+1),true);
00244 
00245    // allocate an additional value for the diagonal, if neccessary
00246    std::vector<double> rowData(stencil.GlobalMaxNumEntries()+1);
00247    std::vector<int> rowInd(stencil.GlobalMaxNumEntries()+1);
00248 
00249    // loop over all the rows
00250    for(int j=0;j<gl->NumMyRows();j++) {
00251       int row = gl->GRID(j);
00252       double diagValue = 0.0;
00253       int diagInd = -1;
00254       int rowSz = 0;
00255 
00256       // extract a copy of this row...put it in rowData, rowIndicies
00257       stencil.ExtractGlobalRowCopy(row,stencil.MaxNumEntries(),rowSz,&rowData[0],&rowInd[0]);
00258  
00259       // loop over elements of row
00260       for(int i=0;i<rowSz;i++) {
00261          int col = rowInd[i];
00262 
00263          // is this a 0 entry masquerading as some thing else?
00264          double value = rowData[i];
00265          if(value==0) continue;
00266 
00267          // for nondiagonal entries
00268          if(row!=col) {
00269             double d = dist(x,y,z,stride,row,col);
00270             rowData[i] = -1.0/d;
00271             diagValue += rowData[i];
00272          }
00273          else 
00274             diagInd = i;
00275       }
00276     
00277       // handle diagonal entry
00278       if(diagInd<0) { // diagonal not in row
00279          rowData[rowSz] = -diagValue;
00280          rowInd[rowSz] = row;
00281          rowSz++;
00282       }
00283       else { // diagonal in row
00284          rowData[diagInd] = -diagValue;
00285          rowInd[diagInd] = row;
00286       }
00287 
00288       // insert row data into graph Laplacian matrix
00289       TEST_FOR_EXCEPT(gl->InsertGlobalValues(row,rowSz,&rowData[0],&rowInd[0]));
00290    }
00291 
00292    gl->FillComplete();
00293 
00294    return gl;
00295 }
00296 
00312 void applyOp(const LinearOp & A,const MultiVector & x,MultiVector & y,double alpha,double beta)
00313 {
00314    // Thyra::apply(*A,Thyra::NONCONJ_ELE,*x,&*y,alpha,beta);
00315    Thyra::apply(*A,Thyra::NOTRANS,*x,y.ptr(),alpha,beta);
00316 }
00317 
00320 void update(double alpha,const MultiVector & x,double beta,MultiVector & y)
00321 {
00322    Teuchos::Array<double> scale;
00323    Teuchos::Array<Teuchos::Ptr<const Thyra::MultiVectorBase<double> > >  vec;
00324 
00325    // build arrays needed for linear combo
00326    scale.push_back(alpha);
00327    vec.push_back(x.ptr());
00328 
00329    // compute linear combination
00330    Thyra::linear_combination<double>(scale,vec,beta,y.ptr());
00331 }
00332 
00334 BlockedLinearOp getUpperTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
00335 {
00336    int rows = blockRowCount(blo);
00337 
00338    TEUCHOS_ASSERT(rows==blockColCount(blo));
00339 
00340    RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
00341    RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
00342 
00343    // allocate new operator
00344    BlockedLinearOp upper = createBlockedOp();
00345  
00346    // build new operator 
00347    upper->beginBlockFill(rows,rows);
00348 
00349    for(int i=0;i<rows;i++) {
00350       // put zero operators on the diagonal
00351       // this gurantees the vector space of
00352       // the new operator are fully defined
00353       RCP<const Thyra::LinearOpBase<double> > zed 
00354             = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
00355       upper->setBlock(i,i,zed);
00356 
00357       for(int j=i+1;j<rows;j++) {
00358          // get block i,j
00359          LinearOp uij = blo->getBlock(i,j);
00360 
00361          // stuff it in U
00362          if(uij!=Teuchos::null)
00363             upper->setBlock(i,j,uij);
00364       }
00365    }
00366    if(callEndBlockFill)
00367       upper->endBlockFill();
00368 
00369    return upper;
00370 }
00371 
00373 BlockedLinearOp getLowerTriBlocks(const BlockedLinearOp & blo,bool callEndBlockFill)
00374 {
00375    int rows = blockRowCount(blo);
00376 
00377    TEUCHOS_ASSERT(rows==blockColCount(blo));
00378 
00379    RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
00380    RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
00381 
00382    // allocate new operator
00383    BlockedLinearOp lower = createBlockedOp();
00384  
00385    // build new operator 
00386    lower->beginBlockFill(rows,rows);
00387 
00388    for(int i=0;i<rows;i++) {
00389       // put zero operators on the diagonal
00390       // this gurantees the vector space of
00391       // the new operator are fully defined
00392       RCP<const Thyra::LinearOpBase<double> > zed 
00393             = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
00394       lower->setBlock(i,i,zed);
00395 
00396       for(int j=0;j<i;j++) {
00397          // get block i,j
00398          LinearOp uij = blo->getBlock(i,j);
00399 
00400          // stuff it in U
00401          if(uij!=Teuchos::null)
00402             lower->setBlock(i,j,uij);
00403       }
00404    }
00405    if(callEndBlockFill)
00406       lower->endBlockFill();
00407 
00408    return lower;
00409 }
00410 
00430 BlockedLinearOp zeroBlockedOp(const BlockedLinearOp & blo)
00431 {
00432    int rows = blockRowCount(blo);
00433 
00434    TEUCHOS_ASSERT(rows==blockColCount(blo)); // assert that matrix is square
00435 
00436    RCP<const Thyra::ProductVectorSpaceBase<double> > range = blo->productRange();
00437    RCP<const Thyra::ProductVectorSpaceBase<double> > domain = blo->productDomain();
00438 
00439    // allocate new operator
00440    BlockedLinearOp zeroOp = createBlockedOp();
00441  
00442    // build new operator 
00443    zeroOp->beginBlockFill(rows,rows);
00444 
00445    for(int i=0;i<rows;i++) {
00446       // put zero operators on the diagonal
00447       // this gurantees the vector space of
00448       // the new operator are fully defined
00449       RCP<const Thyra::LinearOpBase<double> > zed 
00450             = Thyra::zero<double>(range->getBlock(i),domain->getBlock(i));
00451       zeroOp->setBlock(i,i,zed);
00452    }
00453 
00454    return zeroOp;
00455 }
00456 
00458 bool isZeroOp(const LinearOp op)
00459 {
00460    // if operator is null...then its zero!
00461    if(op==Teuchos::null) return true;
00462 
00463    // try to cast it to a zero linear operator
00464    const LinearOp test = rcp_dynamic_cast<const Thyra::ZeroLinearOpBase<double> >(op);
00465 
00466    // if it works...then its zero...otherwise its null
00467    return test!=Teuchos::null;
00468 }
00469 
00478 ModifiableLinearOp getAbsRowSumMatrix(const LinearOp & op)
00479 {
00480    RCP<const Epetra_CrsMatrix> eCrsOp;
00481 
00482    try {
00483       // get Epetra_Operator
00484       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00485 
00486       // cast it to a CrsMatrix
00487       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00488    }
00489    catch (std::exception & e) {
00490       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00491 
00492       *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix\n";
00493       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00494       *out << "           OR\n";
00495       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00496       *out << std::endl;
00497       *out << "*** THROWN EXCEPTION ***\n";
00498       *out << e.what() << std::endl;
00499       *out << "************************\n";
00500       
00501       throw e;
00502    }
00503 
00504    // extract diagonal
00505    const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00506    Epetra_Vector & diag = *ptrDiag;
00507    
00508    // compute absolute value row sum
00509    diag.PutScalar(0.0);
00510    for(int i=0;i<eCrsOp->NumMyRows();i++) {
00511       double * values = 0;
00512       int numEntries;
00513       eCrsOp->ExtractMyRowView(i,numEntries,values);
00514 
00515       // build abs value row sum
00516       for(int j=0;j<numEntries;j++)
00517          diag[i] += std::abs(values[j]);
00518    }
00519 
00520    // build Thyra diagonal operator
00521    return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSum( " + op->getObjectLabel() + " )");
00522 }
00523 
00532 ModifiableLinearOp getAbsRowSumInvMatrix(const LinearOp & op)
00533 {
00534    RCP<const Epetra_CrsMatrix> eCrsOp;
00535 
00536    try {
00537       // get Epetra_Operator
00538       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00539 
00540       // cast it to a CrsMatrix
00541       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00542    }
00543    catch (std::exception & e) {
00544       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00545 
00546       *out << "Teko: getAbsRowSumMatrix requires an Epetra_CrsMatrix\n";
00547       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00548       *out << "           OR\n";
00549       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00550       *out << std::endl;
00551       *out << "*** THROWN EXCEPTION ***\n";
00552       *out << e.what() << std::endl;
00553       *out << "************************\n";
00554       
00555       throw e;
00556    }
00557 
00558    // extract diagonal
00559    const RCP<Epetra_Vector> ptrDiag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00560    Epetra_Vector & diag = *ptrDiag;
00561    
00562    // compute absolute value row sum
00563    diag.PutScalar(0.0);
00564    for(int i=0;i<eCrsOp->NumMyRows();i++) {
00565       double * values = 0;
00566       int numEntries;
00567       eCrsOp->ExtractMyRowView(i,numEntries,values);
00568 
00569       // build abs value row sum
00570       for(int j=0;j<numEntries;j++)
00571          diag[i] += std::abs(values[j]);
00572    }
00573    diag.Reciprocal(diag); // invert entries
00574 
00575    // build Thyra diagonal operator
00576    return Teko::Epetra::thyraDiagOp(ptrDiag,eCrsOp->RowMap(),"absRowSumInv( " + op->getObjectLabel() + " )");
00577 }
00578 
00586 ModifiableLinearOp getLumpedMatrix(const LinearOp & op)
00587 {
00588    RCP<Thyra::VectorBase<double> > ones = Thyra::createMember(op->domain());
00589    RCP<Thyra::VectorBase<double> > diag = Thyra::createMember(op->range());
00590 
00591    // set to all ones
00592    Thyra::assign(ones.ptr(),1.0);
00593 
00594    // compute lumped diagonal
00595    // Thyra::apply(*op,Thyra::NONCONJ_ELE,*ones,&*diag);
00596    Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
00597 
00598    return rcp(new Thyra::DefaultDiagonalLinearOp<double>(diag));
00599 }
00600 
00609 ModifiableLinearOp getInvLumpedMatrix(const LinearOp & op)
00610 {
00611    RCP<Thyra::VectorBase<double> > ones = Thyra::createMember(op->domain());
00612    RCP<Thyra::VectorBase<double> > diag = Thyra::createMember(op->range());
00613 
00614    // set to all ones
00615    Thyra::assign(ones.ptr(),1.0);
00616 
00617    // compute lumped diagonal
00618    // Thyra::apply(*op,Thyra::NONCONJ_ELE,*ones,&*diag);
00619    Thyra::apply(*op,Thyra::NOTRANS,*ones,diag.ptr());
00620    Thyra::reciprocal(diag.ptr(),*diag);
00621 
00622    return rcp(new Thyra::DefaultDiagonalLinearOp<double>(diag));
00623 }
00624 
00636 const ModifiableLinearOp getDiagonalOp(const LinearOp & op)
00637 {
00638    RCP<const Epetra_CrsMatrix> eCrsOp;
00639 
00640    try {
00641       // get Epetra_Operator
00642       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00643 
00644       // cast it to a CrsMatrix
00645       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00646    }
00647    catch (std::exception & e) {
00648       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00649 
00650       *out << "Teko: getDiagonalOp requires an Epetra_CrsMatrix\n";
00651       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00652       *out << "           OR\n";
00653       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00654       *out << std::endl;
00655       *out << "*** THROWN EXCEPTION ***\n";
00656       *out << e.what() << std::endl;
00657       *out << "************************\n";
00658       
00659       throw e;
00660    }
00661 
00662    // extract diagonal
00663    const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00664    TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
00665 
00666    // build Thyra diagonal operator
00667    return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"diag( " + op->getObjectLabel() + " )");
00668 }
00669 
00670 const MultiVector getDiagonal(const LinearOp & op)
00671 {
00672    RCP<const Epetra_CrsMatrix> eCrsOp;
00673 
00674    try {
00675       // get Epetra_Operator
00676       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00677 
00678       // cast it to a CrsMatrix
00679       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00680    }
00681    catch (std::exception & e) {
00682       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00683 
00684       *out << "Teko: getDiagonal requires an Epetra_CrsMatrix\n";
00685       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00686       *out << "           OR\n";
00687       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00688       *out << std::endl;
00689       *out << "*** THROWN EXCEPTION ***\n";
00690       *out << e.what() << std::endl;
00691       *out << "************************\n";
00692       
00693       throw e;
00694    }
00695 
00696    // extract diagonal
00697    const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00698    TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
00699 
00700    return Thyra::create_Vector(diag,Thyra::create_VectorSpace(Teuchos::rcpFromRef(eCrsOp->RowMap())));
00701 }
00702 
00703 const MultiVector getDiagonal(const Teko::LinearOp & A,const DiagonalType & dt)
00704 {
00705    LinearOp diagOp = Teko::getDiagonalOp(A,dt);
00706 
00707    Teuchos::RCP<const Thyra::MultiVectorBase<double> > v = 
00708          Teuchos::rcp_dynamic_cast<const Thyra::DiagonalLinearOpBase<double> >(A)->getDiag(); 
00709    return Teuchos::rcp_const_cast<Thyra::MultiVectorBase<double> >(v);
00710 }
00711 
00723 const ModifiableLinearOp getInvDiagonalOp(const LinearOp & op)
00724 {
00725    RCP<const Epetra_CrsMatrix> eCrsOp;
00726 
00727    try {
00728       // get Epetra_Operator
00729       RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00730 
00731       // cast it to a CrsMatrix
00732       eCrsOp = rcp_dynamic_cast<const Epetra_CrsMatrix>(eOp,true);
00733    }
00734    catch (std::exception & e) {
00735       RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
00736 
00737       *out << "Teko: getInvDiagonalOp requires an Epetra_CrsMatrix\n";
00738       *out << "    Could not extract an Epetra_Operator from a \"" << op->description() << std::endl;
00739       *out << "           OR\n";
00740       *out << "    Could not cast an Epetra_Operator to a Epetra_CrsMatrix\n";
00741       *out << std::endl;
00742       *out << "*** THROWN EXCEPTION ***\n";
00743       *out << e.what() << std::endl;
00744       *out << "************************\n";
00745       
00746       throw e;
00747    }
00748 
00749    // extract diagonal
00750    const RCP<Epetra_Vector> diag = rcp(new Epetra_Vector(eCrsOp->RowMap()));
00751    TEST_FOR_EXCEPT(eCrsOp->ExtractDiagonalCopy(*diag));
00752    diag->Reciprocal(*diag);
00753 
00754    // build Thyra diagonal operator
00755    return Teko::Epetra::thyraDiagOp(diag,eCrsOp->RowMap(),"inv(diag( " + op->getObjectLabel() + " ))");
00756 }
00757 
00770 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr)
00771 {
00772    // build implicit multiply
00773    const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
00774 
00775    // build transformer
00776    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00777        Thyra::epetraExtDiagScaledMatProdTransformer();
00778 
00779    // build operator and multiply
00780    const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
00781    prodTrans->transform(*implicitOp,explicitOp.ptr());
00782    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00783                                      " * " + opm->getObjectLabel() +
00784                                      " * " + opr->getObjectLabel() + " )");
00785 
00786    return explicitOp;
00787 }
00788 
00803 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opm,const LinearOp & opr,
00804                                           const ModifiableLinearOp & destOp)
00805 {
00806    // build implicit multiply
00807    const LinearOp implicitOp = Thyra::multiply(opl,opm,opr);
00808 
00809    // build transformer
00810    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00811        Thyra::epetraExtDiagScaledMatProdTransformer();
00812 
00813    // build operator destination operator
00814    ModifiableLinearOp explicitOp; 
00815 
00816    // if neccessary build a operator to put the explicit multiply into
00817    if(destOp==Teuchos::null)
00818       explicitOp = prodTrans->createOutputOp();
00819    else
00820       explicitOp = destOp;
00821 
00822    // perform multiplication
00823    prodTrans->transform(*implicitOp,explicitOp.ptr());
00824 
00825    // label it
00826    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00827                                      " * " + opm->getObjectLabel() +
00828                                      " * " + opr->getObjectLabel() + " )");
00829 
00830    return explicitOp;
00831 }
00832 
00843 const LinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr)
00844 {
00845    // build implicit multiply
00846    const LinearOp implicitOp = Thyra::multiply(opl,opr);
00847 
00848    // build a scaling transformer
00849    RCP<Thyra::LinearOpTransformerBase<double> > prodTrans  
00850          = Thyra::epetraExtDiagScalingTransformer();
00851 
00852    // check to see if a scaling transformer works: if not use the 
00853    // DiagScaledMatrixProduct transformer
00854    if(not prodTrans->isCompatible(*implicitOp))
00855        prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
00856 
00857    // build operator and multiply
00858    const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
00859    prodTrans->transform(*implicitOp,explicitOp.ptr());
00860    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00861                                      " * " + opr->getObjectLabel() + " )");
00862 
00863    return explicitOp;
00864 }
00865 
00879 const ModifiableLinearOp explicitMultiply(const LinearOp & opl,const LinearOp & opr,
00880                                           const ModifiableLinearOp & destOp)
00881 {
00882    // build implicit multiply
00883    const LinearOp implicitOp = Thyra::multiply(opl,opr);
00884 
00885    // build a scaling transformer
00886    RCP<Thyra::LinearOpTransformerBase<double> > prodTrans  
00887          = Thyra::epetraExtDiagScalingTransformer();
00888 
00889    // check to see if a scaling transformer works: if not use the 
00890    // DiagScaledMatrixProduct transformer
00891    if(not prodTrans->isCompatible(*implicitOp))
00892        prodTrans = Thyra::epetraExtDiagScaledMatProdTransformer();
00893 
00894    // build operator destination operator
00895    ModifiableLinearOp explicitOp; 
00896 
00897    // if neccessary build a operator to put the explicit multiply into
00898    if(destOp==Teuchos::null)
00899       explicitOp = prodTrans->createOutputOp();
00900    else
00901       explicitOp = destOp;
00902 
00903    // perform multiplication
00904    prodTrans->transform(*implicitOp,explicitOp.ptr());
00905 
00906    // label it
00907    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00908                                      " * " + opr->getObjectLabel() + " )");
00909 
00910    return explicitOp;
00911 }
00912 
00923 const LinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr)
00924 {
00925    // build implicit multiply
00926    const LinearOp implicitOp = Thyra::add(opl,opr);
00927  
00928    // build transformer
00929    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00930        Thyra::epetraExtAddTransformer();
00931 
00932    // build operator and multiply
00933    const RCP<Thyra::LinearOpBase<double> > explicitOp = prodTrans->createOutputOp();
00934    prodTrans->transform(*implicitOp,explicitOp.ptr());
00935    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00936                                      " + " + opr->getObjectLabel() + " )");
00937 
00938    return explicitOp;
00939 }
00940 
00953 const ModifiableLinearOp explicitAdd(const LinearOp & opl,const LinearOp & opr,
00954                                      const ModifiableLinearOp & destOp)
00955 {
00956    // build implicit multiply
00957    const LinearOp implicitOp = Thyra::add(opl,opr);
00958  
00959    // build transformer
00960    const RCP<Thyra::LinearOpTransformerBase<double> > prodTrans =
00961        Thyra::epetraExtAddTransformer();
00962 
00963    // build or reuse destination operator
00964    RCP<Thyra::LinearOpBase<double> > explicitOp;
00965    if(destOp!=Teuchos::null)
00966       explicitOp = destOp;
00967    else
00968       explicitOp = prodTrans->createOutputOp();
00969 
00970    // perform multiplication
00971    prodTrans->transform(*implicitOp,explicitOp.ptr());
00972    explicitOp->setObjectLabel("explicit( " + opl->getObjectLabel() +
00973                                      " + " + opr->getObjectLabel() + " )");
00974 
00975    return explicitOp;
00976 }
00977 
00978 const LinearOp explicitTranspose(const LinearOp & op)
00979 {
00980    RCP<const Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*op);
00981    TEST_FOR_EXCEPTION(eOp==Teuchos::null,std::logic_error,
00982                              "Teko::explicitTranspose Not an Epetra_Operator");
00983    RCP<const Epetra_RowMatrix> eRowMatrixOp 
00984          = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(eOp);
00985    TEST_FOR_EXCEPTION(eRowMatrixOp==Teuchos::null,std::logic_error,
00986                              "Teko::explicitTranspose Not an Epetra_RowMatrix");
00987 
00988    // we now have a delete transpose operator
00989    EpetraExt::RowMatrix_Transpose tranposeOp;
00990    Epetra_RowMatrix & eMat = tranposeOp(const_cast<Epetra_RowMatrix &>(*eRowMatrixOp));
00991 
00992    // this copy is because of a poor implementation of the EpetraExt::Transform 
00993    // implementation
00994    Teuchos::RCP<Epetra_CrsMatrix> crsMat 
00995          = Teuchos::rcp(new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix &>(eMat)));
00996 
00997    return Thyra::epetraLinearOp(crsMat);
00998 }
00999 
01000 const LinearOp buildDiagonal(const MultiVector & src,const std::string & lbl)
01001 {
01002    RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(src->range()); 
01003    Thyra::copy(*src->col(0),dst.ptr());
01004    
01005    return Thyra::diagonal<double>(dst,lbl);
01006 }
01007 
01008 const LinearOp buildInvDiagonal(const MultiVector & src,const std::string & lbl)
01009 {
01010    const RCP<const Thyra::VectorBase<double> > srcV = src->col(0);
01011    RCP<Thyra::VectorBase<double> > dst = Thyra::createMember(srcV->range()); 
01012    Thyra::reciprocal<double>(dst.ptr(),*srcV);
01013 
01014    return Thyra::diagonal<double>(dst,lbl);
01015 }
01016 
01018 BlockedMultiVector buildBlockedMultiVector(const std::vector<MultiVector> & mvv)
01019 {
01020    Teuchos::Array<MultiVector> mvA;
01021    Teuchos::Array<VectorSpace> vsA;
01022 
01023    // build arrays of multi vectors and vector spaces
01024    std::vector<MultiVector>::const_iterator itr;
01025    for(itr=mvv.begin();itr!=mvv.end();++itr) {
01026       mvA.push_back(*itr);
01027       vsA.push_back((*itr)->range());
01028    }
01029 
01030    // construct the product vector space
01031    const RCP<const Thyra::DefaultProductVectorSpace<double> > vs
01032          = Thyra::productVectorSpace<double>(vsA);
01033 
01034    return Thyra::defaultProductMultiVector<double>(vs,mvA);
01035 }
01036 
01047 Teuchos::RCP<Thyra::VectorBase<double> > indicatorVector(
01048                                  const std::vector<int> & indices,
01049                                  const VectorSpace & vs,
01050                                  double onValue, double offValue)
01051 
01052 {
01053    using Teuchos::RCP;
01054  
01055    // create a new vector
01056    RCP<Thyra::VectorBase<double> > v = Thyra::createMember<double>(vs);
01057    Thyra::put_scalar<double>(offValue,v.ptr()); // fill it with "off" values
01058 
01059    // set on values
01060    for(std::size_t i=0;i<indices.size();i++) 
01061       Thyra::set_ele<double>(indices[i],onValue,v.ptr());
01062 
01063    return v;
01064 }
01065 
01090 double computeSpectralRad(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
01091                           bool isHermitian,int numBlocks,int restart,int verbosity)
01092 {
01093    typedef Thyra::LinearOpBase<double> OP;
01094    typedef Thyra::MultiVectorBase<double> MV;
01095 
01096    int startVectors = 1;
01097 
01098    // construct an initial guess
01099    const RCP<MV> ivec = Thyra::createMember(A->domain());
01100    Thyra::randomize(-1.0,1.0,&*ivec);
01101    
01102    RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
01103          = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
01104    eigProb->setNEV(1);
01105    eigProb->setHermitian(isHermitian);
01106 
01107    // set the problem up
01108    if(not eigProb->setProblem()) {
01109       // big time failure!
01110       return Teuchos::ScalarTraits<double>::nan();
01111    }
01112 
01113    // we want largert magnitude eigenvalue
01114    std::string which("LM"); // largest magnitude
01115 
01116    // Create the parameter list for the eigensolver
01117    // verbosity+=Anasazi::TimingDetails;
01118    Teuchos::ParameterList MyPL;
01119    MyPL.set( "Verbosity", verbosity );
01120    MyPL.set( "Which", which );
01121    MyPL.set( "Block Size", startVectors );
01122    MyPL.set( "Num Blocks", numBlocks );
01123    MyPL.set( "Maximum Restarts", restart ); 
01124    MyPL.set( "Convergence Tolerance", tol );
01125 
01126    // build status test manager
01127    // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
01128    //       = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(numBlocks*(restart+1)));
01129 
01130    // Create the Block Krylov Schur solver
01131    // This takes as inputs the eigenvalue problem and the solver parameters
01132    Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
01133  
01134    // Solve the eigenvalue problem, and save the return code
01135    Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
01136 
01137    if(solverreturn==Anasazi::Unconverged) {
01138       double real = MyBlockKrylovSchur.getRitzValues().begin()->realpart;
01139       double comp = MyBlockKrylovSchur.getRitzValues().begin()->imagpart;
01140 
01141       return -std::abs(std::sqrt(real*real+comp*comp));
01142 
01143       // cout << "Anasazi::BlockKrylovSchur::solve() did not converge!" << std::endl; 
01144       // return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
01145    }
01146    else { // solverreturn==Anasazi::Converged
01147       double real = eigProb->getSolution().Evals.begin()->realpart;
01148       double comp = eigProb->getSolution().Evals.begin()->imagpart;
01149 
01150       return std::abs(std::sqrt(real*real+comp*comp));
01151 
01152       // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
01153       // return std::abs(eigProb->getSolution().Evals.begin()->realpart);
01154    }
01155 }
01156 
01180 double computeSmallestMagEig(const RCP<const Thyra::LinearOpBase<double> > & A, double tol,
01181                              bool isHermitian,int numBlocks,int restart,int verbosity)
01182 {
01183    typedef Thyra::LinearOpBase<double> OP;
01184    typedef Thyra::MultiVectorBase<double> MV;
01185 
01186    int startVectors = 1;
01187 
01188    // construct an initial guess
01189    const RCP<MV> ivec = Thyra::createMember(A->domain());
01190    Thyra::randomize(-1.0,1.0,&*ivec);
01191    
01192    RCP<Anasazi::BasicEigenproblem<double,MV,OP> > eigProb
01193          = rcp(new Anasazi::BasicEigenproblem<double,MV,OP>(A,ivec));
01194    eigProb->setNEV(1);
01195    eigProb->setHermitian(isHermitian);
01196 
01197    // set the problem up
01198    if(not eigProb->setProblem()) {
01199       // big time failure!
01200       return Teuchos::ScalarTraits<double>::nan();
01201    }
01202 
01203    // we want largert magnitude eigenvalue
01204    std::string which("SM"); // smallest magnitude
01205 
01206    // Create the parameter list for the eigensolver
01207    Teuchos::ParameterList MyPL;
01208    MyPL.set( "Verbosity", verbosity );
01209    MyPL.set( "Which", which );
01210    MyPL.set( "Block Size", startVectors );
01211    MyPL.set( "Num Blocks", numBlocks );
01212    MyPL.set( "Maximum Restarts", restart ); 
01213    MyPL.set( "Convergence Tolerance", tol );
01214 
01215    // build status test manager
01216    // RCP<Anasazi::StatusTestMaxIters<double,MV,OP> > statTest
01217    //       = rcp(new Anasazi::StatusTestMaxIters<double,MV,OP>(10));
01218 
01219    // Create the Block Krylov Schur solver
01220    // This takes as inputs the eigenvalue problem and the solver parameters
01221    Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MyBlockKrylovSchur(eigProb, MyPL );
01222  
01223    // Solve the eigenvalue problem, and save the return code
01224    Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
01225 
01226    if(solverreturn==Anasazi::Unconverged) {
01227       // cout << "Anasazi::BlockKrylovSchur::solve() did not converge! " << std::endl;
01228       return -std::abs(MyBlockKrylovSchur.getRitzValues().begin()->realpart);
01229    }
01230    else { // solverreturn==Anasazi::Converged
01231       // cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
01232       return std::abs(eigProb->getSolution().Evals.begin()->realpart);
01233    }
01234 }
01235 
01244 ModifiableLinearOp getDiagonalOp(const Teko::LinearOp & A,const DiagonalType & dt)
01245 {
01246    switch(dt) {
01247    case Diagonal:
01248       return getDiagonalOp(A);  
01249    case Lumped:
01250       return getLumpedMatrix(A);  
01251    case AbsRowSum:
01252       return getAbsRowSumMatrix(A);  
01253    case NotDiag:
01254    default:
01255       TEST_FOR_EXCEPT(true);
01256    };
01257 
01258    return Teuchos::null;
01259 }
01260 
01269 ModifiableLinearOp getInvDiagonalOp(const Teko::LinearOp & A,const Teko::DiagonalType & dt)
01270 {
01271    switch(dt) {
01272    case Diagonal:
01273       return getInvDiagonalOp(A);  
01274    case Lumped:
01275       return getInvLumpedMatrix(A);  
01276    case AbsRowSum:
01277       return getAbsRowSumInvMatrix(A);  
01278    case NotDiag:
01279    default:
01280       TEST_FOR_EXCEPT(true);
01281    };
01282 
01283    return Teuchos::null;
01284 }
01285 
01292 std::string getDiagonalName(const DiagonalType & dt)
01293 {
01294    switch(dt) {
01295    case Diagonal:
01296       return "Diagonal";
01297    case Lumped:
01298       return "Lumped";
01299    case AbsRowSum:
01300       return "AbsRowSum";
01301    case NotDiag:
01302       return "NotDiag";
01303    case BlkDiag:
01304       return "BlkDiag";
01305    };
01306 
01307    return "<error>";
01308 }
01309 
01318 DiagonalType getDiagonalType(std::string name)
01319 {
01320    if(name=="Diagonal")
01321       return Diagonal;
01322    if(name=="Lumped")
01323      return Lumped;
01324    if(name=="AbsRowSum")
01325       return AbsRowSum;
01326    if(name=="BlkDiag")
01327       return BlkDiag;
01328  
01329    return NotDiag;
01330 }
01331 
01332 LinearOp probe(Teuchos::RCP<const Epetra_CrsGraph> &G,const LinearOp & Op){
01333 #ifdef Teko_ENABLE_Isorropia
01334   Teuchos::ParameterList probeList;
01335   Prober prober(G,probeList,true);
01336   Teuchos::RCP<Epetra_CrsMatrix> Mat=rcp(new Epetra_CrsMatrix(Copy,*G));
01337   Teko::Epetra::EpetraOperatorWrapper Mwrap(Op);
01338   prober.probe(Mwrap,*Mat);
01339   return Thyra::epetraLinearOp(Mat);    
01340 #else
01341   TEST_FOR_EXCEPTION(true,std::runtime_error,"Probe requires Isorropia");
01342 #endif
01343 }
01344 
01345 double norm_1(const MultiVector & v,std::size_t col)
01346 {
01347    Teuchos::Array<double> n(v->domain()->dim());
01348    Thyra::norms_1<double>(*v,n);
01349 
01350    return n[col];
01351 }
01352 
01353 double norm_2(const MultiVector & v,std::size_t col)
01354 {
01355    Teuchos::Array<double> n(v->domain()->dim());
01356    Thyra::norms_2<double>(*v,n);
01357 
01358    return n[col];
01359 }
01360 
01361 }
 All Classes Files Functions Variables