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 "Epetra_Map.h"
00050 #include "Epetra_LocalMap.h"
00051 #include "Epetra_SerialComm.h"
00052 #include "Epetra_CrsMatrix.h"
00053 #include "Teuchos_Assert.hpp"
00054 #include "EpetraExt_ConfigDefs.h"
00055 #include "EpetraExt_MatrixMatrix.h"
00056 #include "EpetraExt_MMHelpers.h"
00057 #include "EpetraExt_Transpose_RowMatrix.h"
00058 
00059 
00060 #include "EpetraExt_RowMatrixOut.h"
00061 
00062 
00063 namespace Thyra {
00064 
00065 
00066 // Overridden from LinearOpTransformerBase
00067 
00068 
00069 bool EpetraExtAddTransformer::isCompatible(
00070       const LinearOpBase<double> &op_in) const
00071 {
00072    TEUCHOS_TEST_FOR_EXCEPT(true);
00073    return false;
00074 }
00075 
00076 
00077 RCP<LinearOpBase<double> >
00078 EpetraExtAddTransformer::createOutputOp() const
00079 {
00080    return nonconstEpetraLinearOp();
00081 }
00082 
00083 
00084 void EpetraExtAddTransformer::transform(
00085    const LinearOpBase<double> &op_in,
00086    const Ptr<LinearOpBase<double> > &op_inout) const
00087 {
00088    using Thyra::unwrap;
00089    using EpetraExt::MatrixMatrix;
00090    using Teuchos::rcp;
00091    using Teuchos::rcp_dynamic_cast;
00092    using Teuchos::dyn_cast;
00093  
00094    //
00095    // A) Get the component Thyra objects for M = op(B) * D * G
00096    //
00097  
00098    const AddedLinearOpBase<double> & add_op =
00099          dyn_cast<const AddedLinearOpBase<double> >(op_in);
00100 
00101 #ifdef TEUCHOS_DEBUG
00102    TEUCHOS_ASSERT_EQUALITY( add_op.numOps(), 2 );
00103 #endif
00104 
00105  
00106    // get properties of first operator: Transpose, scaler multiply...etc
00107    const RCP<const LinearOpBase<double> > op_A = add_op.getOp(0);
00108    double A_scalar = 0.0;
00109    EOpTransp A_transp = NOTRANS;
00110    RCP<const LinearOpBase<double> > A;
00111    unwrap( op_A, &A_scalar, &A_transp, &A );
00112    TEUCHOS_ASSERT(A_transp==NOTRANS || A_transp==CONJTRANS); // sanity check
00113  
00114    // get properties of third operator: Transpose, scaler multiply...etc
00115    const RCP<const LinearOpBase<double> > op_B = add_op.getOp(1);
00116    double B_scalar = 0.0;
00117    EOpTransp B_transp = NOTRANS;
00118    RCP<const LinearOpBase<double> > B;
00119    unwrap( op_B, &B_scalar, &B_transp, &B );
00120    TEUCHOS_ASSERT(B_transp==NOTRANS || B_transp==CONJTRANS); // sanity check
00121 
00122    //
00123    // B) Extract out the Epetra_CrsMatrix objects and the vector
00124    //
00125    
00126   // convert operators to Epetra_CrsMatrix
00127    const RCP<const Epetra_CrsMatrix> epetra_A =
00128          rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*A), true);
00129    const RCP<const Epetra_CrsMatrix> epetra_B =
00130          rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*B), true);
00131    
00132    const Epetra_Map op_inout_row_map 
00133          = (A_transp==CONJTRANS ? epetra_A->ColMap() : epetra_A->RowMap());
00134    const Epetra_Map op_inout_col_map 
00135          = (A_transp==CONJTRANS ? epetra_A->RowMap() : epetra_A->ColMap());
00136  
00137    //
00138    // C) Do the explicit addition
00139    //
00140   
00141    // allocate space for final addition: 3 steps
00142    //   1. Get destination EpetraLinearOp
00143    //   2. Extract RCP to destination Epetra_CrsMatrix
00144    //   3. If neccessary, allocate new Epetra_CrsMatrix
00145    EpetraLinearOp &thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
00146    RCP<Epetra_CrsMatrix>  epetra_op =
00147          rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.epetra_op());
00148    Epetra_CrsMatrix * epetra_op_raw = epetra_op.get();
00149 
00150    // perform addition
00151    const int add_epetra_B_err 
00152       = EpetraExt::MatrixMatrix::Add(*epetra_A,A_transp==CONJTRANS,A_scalar,*epetra_B,B_transp==CONJTRANS,B_scalar,epetra_op_raw);
00153    if(epetra_op==Teuchos::null)
00154       epetra_op = Teuchos::rcp(epetra_op_raw);
00155    
00156    TEUCHOS_ASSERT_EQUALITY( add_epetra_B_err, 0 );
00157    
00158    epetra_op->FillComplete();
00159 
00160    // set output operator to use newly create epetra_op
00161    thyra_epetra_op_inout.initialize(epetra_op);
00162 }
00163 
00164 
00165 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines