Thyra_EpetraExtDiagScalingTransformer.cpp

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_EpetraExtDiagScalingTransformer.hpp"
00031 #include "Thyra_MultipliedLinearOpBase.hpp"
00032 #include "Thyra_DiagonalLinearOpBase.hpp"
00033 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
00034 #include "Thyra_EpetraLinearOp.hpp"
00035 #include "Thyra_get_Epetra_Operator.hpp"
00036 #include "Thyra_EpetraThyraWrappers.hpp"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_LocalMap.h"
00039 #include "Epetra_SerialComm.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include "Teuchos_Assert.hpp"
00042 #include "Teuchos_RCP.hpp"
00043 
00044 
00045 namespace Thyra {
00046 
00047 
00048 // Overridden from LinearOpTransformerBase
00049 
00050 
00051 bool EpetraExtDiagScalingTransformer::isCompatible(
00052       const LinearOpBase<double> &op_in) const
00053 {
00054    using Thyra::unwrap;
00055    using Teuchos::rcp;
00056    using Teuchos::rcp_dynamic_cast;
00057    using Teuchos::dyn_cast;
00058 
00059    const MultipliedLinearOpBase<double> &multi_op =
00060          dyn_cast<const MultipliedLinearOpBase<double> >(op_in);
00061 
00062    // this operation must only have two operands
00063    if(multi_op.numOps()!=2)
00064       return false;
00065 
00066    double scalar = 0.0;
00067    EOpTransp transp = NOTRANS;
00068 
00069    // get properties of first operator: Transpose, scaler multiply...etc
00070    RCP<const LinearOpBase<double> > A;
00071    unwrap( multi_op.getOp(0), &scalar, &transp, &A );
00072    if(transp!=NOTRANS) return false;
00073 
00074    // get properties of second operator: Transpose, scaler multiply...etc
00075    RCP<const LinearOpBase<double> > B;
00076    unwrap( multi_op.getOp(1), &scalar, &transp, &B );
00077    if(transp!=NOTRANS) return false;
00078  
00079    // see if exactly one operator is a diagonal linear op
00080    RCP<const DiagonalLinearOpBase<double> > dA 
00081       = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(A);
00082    RCP<const DiagonalLinearOpBase<double> > dB 
00083       = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(B);
00084 
00085    if(dA==Teuchos::null && dB!=Teuchos::null)
00086       return true;
00087 
00088    if(dA!=Teuchos::null && dB==Teuchos::null)
00089       return true;
00090 
00091    return false;
00092 }
00093 
00094 
00095 RCP<LinearOpBase<double> >
00096 EpetraExtDiagScalingTransformer::createOutputOp() const
00097 {
00098    return nonconstEpetraLinearOp();
00099 }
00100 
00101 
00102 void EpetraExtDiagScalingTransformer::transform(
00103    const LinearOpBase<double> &op_in,
00104    const Ptr<LinearOpBase<double> > &op_inout) const
00105 {
00106    using Thyra::unwrap;
00107    using Teuchos::rcp;
00108    using Teuchos::rcp_dynamic_cast;
00109    using Teuchos::dyn_cast;
00110 
00111    //
00112    // A) Get the component Thyra objects for M = op(A) * op(B)
00113    //
00114  
00115    const MultipliedLinearOpBase<double> &multi_op =
00116          dyn_cast<const MultipliedLinearOpBase<double> >(op_in);
00117 
00118    TEUCHOS_ASSERT(multi_op.numOps()==2);
00119  
00120    // get properties of first operator: Transpose, scaler multiply...etc
00121    const RCP<const LinearOpBase<double> > op_A = multi_op.getOp(0);
00122    double A_scalar = 0.0;
00123    EOpTransp A_transp = NOTRANS;
00124    RCP<const LinearOpBase<double> > A;
00125    unwrap( op_A, &A_scalar, &A_transp, &A );
00126    TEUCHOS_ASSERT(A_transp==NOTRANS || A_transp==CONJTRANS); // sanity check
00127  
00128    // get properties of third operator: Transpose, scaler multiply...etc
00129    const RCP<const LinearOpBase<double> > op_B = multi_op.getOp(1);
00130    double B_scalar = 0.0;
00131    EOpTransp B_transp = NOTRANS;
00132    RCP<const LinearOpBase<double> > B;
00133    unwrap( op_B, &B_scalar, &B_transp, &B );
00134    TEUCHOS_ASSERT(B_transp==NOTRANS || B_transp==CONJTRANS); // sanity check
00135 
00136    //
00137    // B) Extract out the CrsMatrix and Diagonal objects
00138    //
00139 
00140    // see if exactly one operator is a diagonal linear op
00141    RCP<const DiagonalLinearOpBase<double> > dA 
00142       = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(A);
00143    RCP<const DiagonalLinearOpBase<double> > dB 
00144       = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(B);
00145 
00146    RCP<const Epetra_CrsMatrix> epetra_A;
00147    RCP<const Epetra_CrsMatrix> epetra_B;
00148    if(dA==Teuchos::null) {
00149       bool exactly_one_op_must_be_diagonal__dB_neq_null = dB!=Teuchos::null;
00150       TEUCHOS_ASSERT(exactly_one_op_must_be_diagonal__dB_neq_null);
00151 
00152       // convert second operator to an Epetra_CrsMatrix
00153       epetra_A = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*A), true);
00154    }
00155    else if(dB==Teuchos::null) {
00156       bool exactly_one_op_must_be_diagonal__dA_neq_null = dA!=Teuchos::null;
00157       TEUCHOS_ASSERT(exactly_one_op_must_be_diagonal__dA_neq_null);
00158 
00159       // convert second operator to an Epetra_CrsMatrix
00160       epetra_B = rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*B), true);
00161    }
00162    else {
00163       bool exactly_one_op_must_be_diagonal=false;
00164       TEUCHOS_ASSERT(exactly_one_op_must_be_diagonal);
00165    }
00166 
00167    TEUCHOS_ASSERT( A_transp == NOTRANS ); // ToDo: Handle the transpose
00168    TEUCHOS_ASSERT( B_transp == NOTRANS ); // ToDo: Handle the transpose
00169 
00170 
00171    // get or allocate an object to use as the destination
00172    EpetraLinearOp & thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
00173    RCP<Epetra_CrsMatrix>  epetra_op =
00174          rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.epetra_op());
00175 
00176    bool rightScale = dB!=Teuchos::null;
00177 
00178    if(rightScale) {
00179       RCP<const Epetra_Vector> v = get_Epetra_Vector(epetra_A->OperatorDomainMap(),dB->getDiag());
00180 
00181       // notice that this does not reuse anything! => could be more efficient
00182       epetra_op = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
00183       epetra_op->RightScale(*v);
00184    }
00185    else {
00186       RCP<const Epetra_Vector> v = get_Epetra_Vector(epetra_B->OperatorRangeMap(),dA->getDiag());
00187 
00188       // notice that this does not reuse anything! => could be more efficient
00189       epetra_op = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_B));
00190       epetra_op->LeftScale(*v);
00191    }
00192 
00193    epetra_op->Scale(A_scalar*B_scalar);
00194 
00195    // set output operator to use newly create epetra_op
00196    thyra_epetra_op_inout.initialize(epetra_op);
00197 }
00198 
00199 
00200 } // namespace Thyra

Generated on Thu Feb 11 16:28:43 2010 for EpetraExt/Thyra Adapters by  doxygen 1.4.7