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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 
00030 #include "Thyra_EpetraExtAddTransformer.hpp"
00031 #include "Thyra_AddedLinearOpBase.hpp"
00032 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
00033 #include "Thyra_EpetraLinearOp.hpp"
00034 #include "Thyra_get_Epetra_Operator.hpp"
00035 #include "Thyra_EpetraThyraWrappers.hpp"
00036 #include "Epetra_Map.h"
00037 #include "Epetra_LocalMap.h"
00038 #include "Epetra_SerialComm.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Teuchos_Assert.hpp"
00041 #include "EpetraExt_ConfigDefs.h"
00042 #include "EpetraExt_MatrixMatrix.h"
00043 #include "EpetraExt_MMHelpers.h"
00044 #include "EpetraExt_Transpose_RowMatrix.h"
00045 
00046 
00047 #include "EpetraExt_RowMatrixOut.h"
00048 
00049 
00050 namespace Thyra {
00051 
00052 
00053 // Overridden from LinearOpTransformerBase
00054 
00055 
00056 bool EpetraExtAddTransformer::isCompatible(
00057       const LinearOpBase<double> &op_in) const
00058 {
00059    TEST_FOR_EXCEPT(true);
00060    return false;
00061 }
00062 
00063 
00064 RCP<LinearOpBase<double> >
00065 EpetraExtAddTransformer::createOutputOp() const
00066 {
00067    return nonconstEpetraLinearOp();
00068 }
00069 
00070 
00071 void EpetraExtAddTransformer::transform(
00072    const LinearOpBase<double> &op_in,
00073    const Ptr<LinearOpBase<double> > &op_inout) const
00074 {
00075    using Thyra::unwrap;
00076    using EpetraExt::MatrixMatrix;
00077    using Teuchos::rcp;
00078    using Teuchos::rcp_dynamic_cast;
00079    using Teuchos::dyn_cast;
00080  
00081    //
00082    // A) Get the component Thyra objects for M = op(B) * D * G
00083    //
00084  
00085    const AddedLinearOpBase<double> & add_op =
00086          dyn_cast<const AddedLinearOpBase<double> >(op_in);
00087 
00088 #ifdef TEUCHOS_DEBUG
00089    TEUCHOS_ASSERT_EQUALITY( add_op.numOps(), 2 );
00090 #endif
00091 
00092  
00093    // get properties of first operator: Transpose, scaler multiply...etc
00094    const RCP<const LinearOpBase<double> > op_A = add_op.getOp(0);
00095    double A_scalar = 0.0;
00096    EOpTransp A_transp = NOTRANS;
00097    RCP<const LinearOpBase<double> > A;
00098    unwrap( op_A, &A_scalar, &A_transp, &A );
00099    TEUCHOS_ASSERT(A_transp==NOTRANS || A_transp==CONJTRANS); // sanity check
00100  
00101    // get properties of third operator: Transpose, scaler multiply...etc
00102    const RCP<const LinearOpBase<double> > op_B = add_op.getOp(1);
00103    double B_scalar = 0.0;
00104    EOpTransp B_transp = NOTRANS;
00105    RCP<const LinearOpBase<double> > B;
00106    unwrap( op_B, &B_scalar, &B_transp, &B );
00107    TEUCHOS_ASSERT(B_transp==NOTRANS || B_transp==CONJTRANS); // sanity check
00108 
00109    //
00110    // B) Extract out the Epetra_CrsMatrix objects and the vector
00111    //
00112    
00113   // convert operators to Epetra_CrsMatrix
00114    const RCP<const Epetra_CrsMatrix> epetra_A =
00115          rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*A), true);
00116    const RCP<const Epetra_CrsMatrix> epetra_B =
00117          rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*B), true);
00118    
00119    const Epetra_Map op_inout_row_map 
00120          = (A_transp==CONJTRANS ? epetra_A->ColMap() : epetra_A->RowMap());
00121    const Epetra_Map op_inout_col_map 
00122          = (A_transp==CONJTRANS ? epetra_A->RowMap() : epetra_A->ColMap());
00123  
00124    //
00125    // C) Do the explicit addition
00126    //
00127   
00128    // allocate space for final addition: 3 steps
00129    //   1. Get destination EpetraLinearOp
00130    //   2. Extract RCP to destination Epetra_CrsMatrix
00131    //   3. If neccessary, allocate new Epetra_CrsMatrix
00132    EpetraLinearOp &thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
00133    RCP<Epetra_CrsMatrix>  epetra_op =
00134          rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.epetra_op());
00135    Epetra_CrsMatrix * epetra_op_raw = epetra_op.get();
00136 
00137    // perform addition
00138    const int add_epetra_B_err 
00139       = EpetraExt::MatrixMatrix::Add(*epetra_A,A_transp==CONJTRANS,A_scalar,*epetra_B,B_transp==CONJTRANS,B_scalar,epetra_op_raw);
00140    if(epetra_op==Teuchos::null)
00141       epetra_op = Teuchos::rcp(epetra_op_raw);
00142    
00143    TEUCHOS_ASSERT_EQUALITY( add_epetra_B_err, 0 );
00144    
00145    epetra_op->FillComplete();
00146 
00147    // set output operator to use newly create epetra_op
00148    thyra_epetra_op_inout.initialize(epetra_op);
00149 }
00150 
00151 
00152 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines