Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_EpetraExtAddTransformer.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 
00043 #include "Thyra_EpetraExtAddTransformer.hpp"
00044 #include "Thyra_AddedLinearOpBase.hpp"
00045 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
00046 #include "Thyra_EpetraLinearOp.hpp"
00047 #include "Thyra_get_Epetra_Operator.hpp"
00048 #include "Thyra_EpetraThyraWrappers.hpp"
00049 #include "Thyra_DiagonalLinearOpBase.hpp"
00050 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00051 #include "Thyra_IdentityLinearOpBase.hpp"
00052 #include "Thyra_VectorStdOps.hpp"
00053 #include "Epetra_Map.h"
00054 #include "Epetra_LocalMap.h"
00055 #include "Epetra_SerialComm.h"
00056 #include "Epetra_Vector.h"
00057 #include "Epetra_CrsMatrix.h"
00058 #include "Teuchos_Assert.hpp"
00059 #include "EpetraExt_ConfigDefs.h"
00060 #include "EpetraExt_MatrixMatrix.h"
00061 #include "EpetraExt_MMHelpers.h"
00062 #include "EpetraExt_Transpose_RowMatrix.h"
00063 
00064 
00065 #include "EpetraExt_RowMatrixOut.h"
00066 
00067 
00068 namespace Thyra {
00069 
00070 
00071 // Overridden from LinearOpTransformerBase
00072 
00073 
00074 bool EpetraExtAddTransformer::isCompatible(
00075       const LinearOpBase<double> &op_in) const
00076 {
00077    TEUCHOS_TEST_FOR_EXCEPT(true);
00078    return false;
00079 }
00080 
00081 
00082 RCP<LinearOpBase<double> >
00083 EpetraExtAddTransformer::createOutputOp() const
00084 {
00085    return nonconstEpetraLinearOp();
00086 }
00087 
00088 
00089 void EpetraExtAddTransformer::transform(
00090    const LinearOpBase<double> &op_in,
00091    const Ptr<LinearOpBase<double> > &op_inout) const
00092 {
00093    using Thyra::unwrap;
00094    using EpetraExt::MatrixMatrix;
00095    using Teuchos::rcp;
00096    using Teuchos::rcp_dynamic_cast;
00097    using Teuchos::dyn_cast;
00098  
00099    //
00100    // A) Get the component Thyra objects
00101    //
00102  
00103    const AddedLinearOpBase<double> & add_op =
00104          dyn_cast<const AddedLinearOpBase<double> >(op_in);
00105 
00106 #ifdef TEUCHOS_DEBUG
00107    TEUCHOS_ASSERT_EQUALITY( add_op.numOps(), 2 );
00108 #endif
00109 
00110  
00111    // get properties of first operator: Transpose, scaler multiply...etc
00112    const RCP<const LinearOpBase<double> > op_A = add_op.getOp(0);
00113    double A_scalar = 0.0;
00114    EOpTransp A_transp = NOTRANS;
00115    RCP<const LinearOpBase<double> > A;
00116    unwrap( op_A, &A_scalar, &A_transp, &A );
00117    TEUCHOS_ASSERT(A_transp==NOTRANS || A_transp==CONJTRANS); // sanity check
00118  
00119    // get properties of third operator: Transpose, scaler multiply...etc
00120    const RCP<const LinearOpBase<double> > op_B = add_op.getOp(1);
00121    double B_scalar = 0.0;
00122    EOpTransp B_transp = NOTRANS;
00123    RCP<const LinearOpBase<double> > B;
00124    unwrap( op_B, &B_scalar, &B_transp, &B );
00125    TEUCHOS_ASSERT(B_transp==NOTRANS || B_transp==CONJTRANS); // sanity check
00126 
00127    //
00128    // B) Extract out the Epetra_CrsMatrix objects and the vector
00129    //
00130 
00131    // first makre sure identity operators are represented as digaonal vectors
00132    if(rcp_dynamic_cast<const Thyra::IdentityLinearOpBase<double> >(A)!=Teuchos::null) {
00133      RCP<Thyra::VectorBase<double> > d = Thyra::createMember(A->domain(), "d");
00134      Thyra::V_S( d.ptr(), 1.0 ); // ToDo: Set ton != 1.0 and generalize
00135      A = Thyra::diagonal(d);
00136    }
00137    if(rcp_dynamic_cast<const Thyra::IdentityLinearOpBase<double> >(B)!=Teuchos::null) {
00138      RCP<Thyra::VectorBase<double> > d = Thyra::createMember(B->domain(), "d");
00139      Thyra::V_S( d.ptr(), 1.0 ); // ToDo: Set ton != 1.0 and generalize
00140      B = Thyra::diagonal(d);
00141    }
00142       
00143 
00144    // see if exactly one operator is a diagonal linear op
00145    RCP<const DiagonalLinearOpBase<double> > dA 
00146       = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(A);
00147    RCP<const DiagonalLinearOpBase<double> > dB 
00148       = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(B);
00149 
00150    // convert operators to Epetra_CrsMatrix
00151    RCP<const Epetra_CrsMatrix> epetra_A;
00152    RCP<const Epetra_CrsMatrix> epetra_B;
00153    if(dA==Teuchos::null)
00154       epetra_A = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*A), true);
00155    if(dB==Teuchos::null)
00156       epetra_B = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*B), true);
00157 
00158    //
00159    // C) Do the explicit addition
00160    //
00161 
00162    if(epetra_A!=Teuchos::null && epetra_B!=Teuchos::null) {
00163   
00164       // allocate space for final addition: 3 steps
00165       //   1. Get destination EpetraLinearOp
00166       //   2. Extract RCP to destination Epetra_CrsMatrix
00167       //   3. If neccessary, allocate new Epetra_CrsMatrix
00168       EpetraLinearOp &thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
00169       RCP<Epetra_CrsMatrix>  epetra_op =
00170             rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.epetra_op());
00171       Epetra_CrsMatrix * epetra_op_raw = epetra_op.get();
00172    
00173       // perform addition
00174       const int add_epetra_B_err 
00175          = EpetraExt::MatrixMatrix::Add(*epetra_A,A_transp==CONJTRANS,A_scalar,*epetra_B,B_transp==CONJTRANS,B_scalar,epetra_op_raw);
00176       if(epetra_op==Teuchos::null)
00177          epetra_op = Teuchos::rcp(epetra_op_raw);
00178       
00179       TEUCHOS_ASSERT_EQUALITY( add_epetra_B_err, 0 );
00180       
00181       epetra_op->FillComplete();
00182    
00183       // set output operator to use newly create epetra_op
00184       thyra_epetra_op_inout.initialize(epetra_op);
00185    }
00186    else if((dA!=Teuchos::null && epetra_B!=Teuchos::null) ||
00187            (dB!=Teuchos::null && epetra_A!=Teuchos::null)) { 
00188 
00189       // get unique addition values
00190       RCP<const Epetra_CrsMatrix> crsMat = (dA!=Teuchos::null) ? epetra_B : epetra_A;
00191       double matScalar = (dA!=Teuchos::null) ? B_scalar : A_scalar; 
00192       RCP<const DiagonalLinearOpBase<double> > diag = (dA!=Teuchos::null) ? dA : dB; 
00193       double diagScalar = (dA!=Teuchos::null) ? A_scalar : B_scalar; 
00194 
00195       TEUCHOS_ASSERT(crsMat!=Teuchos::null);
00196       TEUCHOS_ASSERT(diag!=Teuchos::null);
00197 
00198       // get or allocate an object to use as the destination
00199       EpetraLinearOp & thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
00200       RCP<Epetra_CrsMatrix>  epetra_op =
00201             rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.epetra_op());
00202 
00203       if(epetra_op==Teuchos::null)
00204          epetra_op = Teuchos::rcp(new Epetra_CrsMatrix(*crsMat));
00205       else
00206          *epetra_op = *crsMat;
00207       
00208       // grab vector to add to diagonal
00209       RCP<const Epetra_Vector> v = get_Epetra_Vector(epetra_op->OperatorDomainMap(),diag->getDiag());
00210 
00211       if(matScalar!=1.0)
00212          epetra_op->Scale(matScalar);
00213 
00214       // grab digaonal from matrix, do summation, then replace the values
00215       RCP<Epetra_Vector> diagonal = rcp(new Epetra_Vector(epetra_op->OperatorDomainMap()));
00216       TEUCHOS_TEST_FOR_EXCEPTION(epetra_op->ExtractDiagonalCopy(*diagonal),std::runtime_error,
00217                                  "Thyra::EpetraExtractAddTransformer::transform ExtractDiagonalCopy failed!");;    
00218       diagonal->Update(diagScalar,*v,1.0); // no need to scale matrix, already scaled
00219       TEUCHOS_TEST_FOR_EXCEPTION(epetra_op->ReplaceDiagonalValues(*diagonal),std::runtime_error,
00220                                  "Thyra::EpetraExtractAddTransformer::transform ReplaceDiagonalValues failed!");;    
00221 
00222       // set output operator to use newly create epetra_op
00223       thyra_epetra_op_inout.initialize(epetra_op);
00224    }
00225    else {
00226       TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
00227                                  "Your case of adding Epetra operators is not yet implemented! Contact the Thyra developers.");
00228    }
00229 }
00230 
00231 
00232 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines