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