Thyra_EpetraExtDiagScaledMatProdTransformer.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_EpetraExtDiagScaledMatProdTransformer.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 "EpetraExt_MatrixMatrix.h"
00042 #include "Teuchos_Assert.hpp"
00043 
00044 
00045 namespace Thyra {
00046 
00047 
00048 // Overridden from LinearOpTransformerBase
00049 
00050 
00051 bool EpetraExtDiagScaledMatProdTransformer::isCompatible(
00052       const LinearOpBase<double> &op_in) const
00053 {
00054    TEST_FOR_EXCEPT(true);
00055    return false;
00056 }
00057 
00058 
00059 RCP<LinearOpBase<double> >
00060 EpetraExtDiagScaledMatProdTransformer::createOutputOp() const
00061 {
00062    return nonconstEpetraLinearOp();
00063 }
00064 
00065 
00066 void EpetraExtDiagScaledMatProdTransformer::transform(
00067    const LinearOpBase<double> &op_in,
00068    const Ptr<LinearOpBase<double> > &op_inout) const
00069 {
00070    using Thyra::unwrap;
00071    using EpetraExt::MatrixMatrix;
00072    using Teuchos::rcp;
00073    using Teuchos::rcp_dynamic_cast;
00074    using Teuchos::dyn_cast;
00075  
00076    //
00077    // A) Get the component Thyra objects for M = op(B) * D * G
00078    //
00079  
00080    const MultipliedLinearOpBase<double> &multi_op =
00081          dyn_cast<const MultipliedLinearOpBase<double> >(op_in);
00082 
00083    bool haveDiagScaling = (multi_op.numOps()==3);
00084  
00085    // get properties of first operator: Transpose, scaler multiply...etc
00086    const RCP<const LinearOpBase<double> > op_B = multi_op.getOp(0);
00087    double B_scalar = 0.0;
00088    EOpTransp B_transp = NOTRANS;
00089    RCP<const LinearOpBase<double> > B;
00090    unwrap( op_B, &B_scalar, &B_transp, &B );
00091    TEUCHOS_ASSERT(B_transp==NOTRANS || B_transp==CONJTRANS); // sanity check
00092  
00093    // get diagonal scaling
00094    RCP<const VectorBase<double> > d;
00095    double D_scalar = 1.0;
00096    if(haveDiagScaling) {
00097       const RCP<const LinearOpBase<double> > op_D = multi_op.getOp(1);
00098       EOpTransp D_transp = NOTRANS;
00099       RCP<const LinearOpBase<double> > D;
00100       unwrap( op_D, &D_scalar, &D_transp, &D );
00101       d = rcp_dynamic_cast<const DiagonalLinearOpBase<double> >(D, true)->getDiag();
00102    }
00103  
00104    // get properties of third operator: Transpose, scaler multiply...etc
00105    const RCP<const LinearOpBase<double> > op_G = multi_op.getOp(haveDiagScaling ? 2 : 1);
00106    double G_scalar = 0.0;
00107    EOpTransp G_transp = NOTRANS;
00108    RCP<const LinearOpBase<double> > G;
00109    unwrap( op_G, &G_scalar, &G_transp, &G );
00110    TEUCHOS_ASSERT(G_transp==NOTRANS || G_transp==CONJTRANS); // sanity check
00111 
00112    //
00113    // B) Extract out the Epetra_CrsMatrix objects and the vector
00114    //
00115    
00116    // convert second operator to an Epetra_CrsMatrix
00117    const RCP<const Epetra_CrsMatrix> epetra_B =
00118          rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*B), true);
00119    // TEUCHOS_ASSERT( B_transp == NOTRANS ); // ToDo: Handle the transpose
00120    
00121    // extract dagonal
00122    RCP<const Epetra_Vector> epetra_d;
00123    if(haveDiagScaling) {
00124       epetra_d = (B_transp==CONJTRANS ? get_Epetra_Vector(epetra_B->OperatorRangeMap(), d)
00125                                       : get_Epetra_Vector(epetra_B->OperatorDomainMap(), d));
00126    }
00127  
00128    // convert third operator to an Epetra_CrsMatrix
00129    const RCP<const Epetra_CrsMatrix> epetra_G =
00130      rcp_dynamic_cast<const Epetra_CrsMatrix>(get_Epetra_Operator(*G), true);
00131    
00132    // determine row map for final operator
00133    const Epetra_Map op_inout_row_map 
00134          = (B_transp==CONJTRANS ? epetra_B->ColMap() : epetra_B->RowMap());
00135    const Epetra_Map op_inout_col_map 
00136          = (G_transp==CONJTRANS ? epetra_B->RowMap() : epetra_B->ColMap());
00137  
00138    //
00139    // C) Do the explicit multiplication
00140    //
00141   
00142    // allocate space for final product: 3 steps
00143    //   1. Get destination EpetraLinearOp
00144    //   2. Extract RCP to destination Epetra_CrsMatrix
00145    //   3. If neccessary, allocate new Epetra_CrsMatrix
00146    EpetraLinearOp &thyra_epetra_op_inout = dyn_cast<EpetraLinearOp>(*op_inout);
00147    RCP<Epetra_CrsMatrix>  epetra_op =
00148          rcp_dynamic_cast<Epetra_CrsMatrix>(thyra_epetra_op_inout.epetra_op());
00149    if(is_null(epetra_op)) {
00150       epetra_op = Teuchos::rcp(
00151             new Epetra_CrsMatrix(::Copy, op_inout_row_map, 0));
00152    }
00153     
00154    // if necessary scale B by diagonal vector
00155    RCP<const Epetra_CrsMatrix> epetra_BD;
00156    if(haveDiagScaling) {
00157       // create a temporary to get around const issue
00158       RCP<Epetra_CrsMatrix> epetra_BD_temp = rcp(new Epetra_CrsMatrix(*epetra_B));
00159 
00160       // scale matrix depending on properties of B
00161       if(B_transp==CONJTRANS) 
00162          epetra_BD_temp->LeftScale(*epetra_d);
00163       else
00164          epetra_BD_temp->RightScale(*epetra_d);
00165  
00166       epetra_BD = epetra_BD_temp;
00167    }
00168    else
00169       epetra_BD = epetra_B;
00170  
00171    // perform multiply
00172    TEUCHOS_ASSERT(
00173        MatrixMatrix::Multiply( *epetra_BD,  B_transp==CONJTRANS,
00174                                *epetra_G,   G_transp==CONJTRANS, *epetra_op)==0);
00175 
00176    // scale the whole thing if neccessary
00177    if(B_scalar*G_scalar*D_scalar!=1.0)
00178       epetra_op->Scale(B_scalar*G_scalar*D_scalar);
00179  
00180    // set output operator to use newly create epetra_op
00181    thyra_epetra_op_inout.initialize(epetra_op);
00182 }
00183 
00184 
00185 } // namespace Thyra

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