Thyra_LinearOperatorImpl.hpp

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 #ifndef THYRA_LINEAROPERATOR_IMPL_HPP
00030 #define THYRA_LINEAROPERATOR_IMPL_HPP
00031 
00032 #include "Thyra_LinearOperatorDecl.hpp"
00033 #include "Thyra_ConfigDefs.hpp"
00034 #include "Thyra_VectorSpaceImpl.hpp"
00035 #include "Thyra_BlockedLinearOpBase.hpp"
00036 #include "Thyra_DefaultBlockedLinearOp.hpp"
00037 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00038 #include "Thyra_DefaultAddedLinearOp.hpp"
00039 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00040 
00041 namespace Thyra
00042 {
00043 
00044   template <class RangeScalar, class DomainScalar> inline 
00045   const VectorSpace<DomainScalar> ConstLinearOperator<RangeScalar, DomainScalar>
00046   ::domain() const 
00047   {return this->constPtr()->domain();}
00048     
00049   template <class RangeScalar, class DomainScalar> inline 
00050   const VectorSpace<RangeScalar> ConstLinearOperator<RangeScalar, DomainScalar>
00051   ::range() const 
00052   {return this->constPtr()->range();}
00053 
00054   template <class RangeScalar, class DomainScalar> inline 
00055   void ConstLinearOperator<RangeScalar, DomainScalar>
00056   ::apply(const ConstVector<DomainScalar>& in,
00057           Vector<RangeScalar>& out,
00058           const RangeScalar& alpha,
00059           const RangeScalar& beta) const
00060   {
00061     /* the result vector might not be initialized. If it's null,
00062      * create a new vector in the range space */
00063     if (out.ptr().get()==0)
00064       {
00065         out = this->range().createMember();
00066       }
00067     this->constPtr()->apply(NONCONJ_ELE, *(in.constPtr().get()),
00068                        out.ptr().get(), alpha, beta);
00069   }
00070 
00071   template <class RangeScalar, class DomainScalar> inline 
00072   void ConstLinearOperator<RangeScalar, DomainScalar>
00073   ::applyTranspose(const ConstVector<RangeScalar>& in,
00074                    Vector<DomainScalar>& out,
00075                    const DomainScalar& alpha,
00076                    const DomainScalar& beta) const
00077   {
00078     /* the result vector might not be initialized. If it's null,
00079      * create a new vector in the domain space (i.e., the range space
00080      * of the transpose operator */
00081     if (out.ptr().get()==0)
00082       {
00083         out = this->domain().createMember();
00084       }
00085     this->constPtr()->applyTranspose(NONCONJ_ELE, *(in.constPtr().get()),
00086                                 out.ptr().get(), 
00087                                 alpha, beta);
00088   }
00089   
00090   template <class RangeScalar, class DomainScalar> inline 
00091   int ConstLinearOperator<RangeScalar, DomainScalar>::numBlockRows() const
00092   {
00093     return range().numBlocks();
00094   }
00095   
00096   template <class RangeScalar, class DomainScalar> inline 
00097   int ConstLinearOperator<RangeScalar, DomainScalar>::numBlockCols() const
00098   {
00099     return domain().numBlocks();
00100   }
00101 
00102   template <class RangeScalar, class DomainScalar> inline 
00103   ConstLinearOperator<RangeScalar, DomainScalar> 
00104   ConstLinearOperator<RangeScalar, DomainScalar>::getBlock(int blockRow, 
00105                                                            int blockCol) const
00106   {
00107     const Thyra::BlockedLinearOpBase<RangeScalar, DomainScalar>* p = 
00108       dynamic_cast<const Thyra::BlockedLinearOpBase<RangeScalar, DomainScalar>* >(this->constPtr().get());
00109     TEST_FOR_EXCEPTION(p==0 && blockRow != 0, runtime_error,
00110                        "request for block row=" << blockRow << " in a non-block operator");
00111     
00112     TEST_FOR_EXCEPTION(p==0 && blockCol != 0, runtime_error,
00113                        "request for block col=" << blockCol << " in a non-block operator");
00114     
00115     if (p != 0)
00116       {
00117         return p->getBlock(blockRow, blockCol);
00118       }
00119     return *this;
00120   }
00121 
00122   template <class RangeScalar, class DomainScalar> inline 
00123   LinearOperator<RangeScalar, DomainScalar> 
00124   LinearOperator<RangeScalar, DomainScalar>::getBlock(int blockRow, 
00125                                                       int blockCol) 
00126   {
00127     Thyra::BlockedLinearOpBase<RangeScalar, DomainScalar>* p = 
00128       dynamic_cast<Thyra::BlockedLinearOpBase<RangeScalar, DomainScalar>* >(this->ptr().get());
00129     TEST_FOR_EXCEPTION(p==0 && blockRow != 0, runtime_error,
00130                        "request for block row=" << blockRow << " in a non-block operator");
00131     
00132     TEST_FOR_EXCEPTION(p==0 && blockCol != 0, runtime_error,
00133                        "request for block col=" << blockCol << " in a non-block operator");
00134     
00135     if (p != 0)
00136       {
00137         return p->getNonconstBlock(blockRow, blockCol);
00138       }
00139     return *this;
00140   }
00141 
00142   //
00143   // Nonmember functions
00144   //
00145 
00146   template <class Scalar> inline 
00147   ConstLinearOperator<Scalar>
00148   operator*(const Scalar& a, 
00149             const ConstLinearOperator<Scalar>& A)
00150   {
00151     return Thyra::scale<Scalar>(a,A.constPtr());
00152   }
00153 
00154   template <class Scalar> inline 
00155   LinearOperator<Scalar>
00156   operator*(const Scalar& a, 
00157             const LinearOperator<Scalar>& A)
00158   {
00159     return Thyra::nonconstScale<Scalar>(a,A.ptr());
00160   }
00161   
00162   template <class Scalar> inline 
00163   ConstLinearOperator<Scalar>
00164   operator*(const ConstLinearOperator<Scalar>& A,
00165             const ConstLinearOperator<Scalar>& B)
00166   {
00167     return Thyra::multiply<Scalar>(A.constPtr(),B.constPtr());
00168   }
00169   
00170   template <class Scalar> inline 
00171   LinearOperator<Scalar>
00172   operator*(const LinearOperator<Scalar>& A,
00173             const LinearOperator<Scalar>& B)
00174   {
00175     return Thyra::nonconstMultiply<Scalar>(A.ptr(),B.ptr());
00176   }
00177  
00178   template <class Scalar> inline 
00179   ConstLinearOperator<Scalar>
00180   operator+(const ConstLinearOperator<Scalar>& A,
00181             const ConstLinearOperator<Scalar>& B)
00182   {
00183     return Thyra::add<Scalar>(A.constPtr(),B.constPtr());
00184   }
00185   
00186   template <class Scalar> inline 
00187   LinearOperator<Scalar>
00188   operator+(const LinearOperator<Scalar>& A,
00189             const LinearOperator<Scalar>& B)
00190   {
00191     return Thyra::nonconstAdd<Scalar>(A.ptr(),B.ptr());
00192   }
00193   
00194   template <class Scalar> inline 
00195   ConstLinearOperator<Scalar>
00196   block2x2(const ConstLinearOperator<Scalar>& A00,
00197            const ConstLinearOperator<Scalar>& A01,
00198            const ConstLinearOperator<Scalar>& A10,
00199            const ConstLinearOperator<Scalar>& A11)
00200   {
00201     return block2x2(A00.constPtr(),A01.constPtr(),A10.constPtr(),A11.constPtr());
00202   }
00203 
00204   template <class Scalar> inline 
00205   ConstLinearOperator<Scalar>
00206   block2x1(const ConstLinearOperator<Scalar>& A00,
00207            const ConstLinearOperator<Scalar>& A10)
00208   {
00209     return block2x1(A00.constPtr(), A10.constPtr());
00210   }
00211 
00212   template <class Scalar> inline 
00213   ConstLinearOperator<Scalar>
00214   block1x2(const ConstLinearOperator<Scalar>& A00,
00215            const ConstLinearOperator<Scalar>& A01)
00216   {
00217     return block2x1(A00.constPtr(), A01.constPtr());
00218   }
00219 
00220   template <class Scalar> inline 
00221   LinearOperator<Scalar>
00222   block2x2(const LinearOperator<Scalar>& A00,
00223            const LinearOperator<Scalar>& A01,
00224            const LinearOperator<Scalar>& A10,
00225            const LinearOperator<Scalar>& A11)
00226   {
00227     return nonconstBlock2x2(A00.ptr(),A01.ptr(),A10.ptr(),A11.ptr());
00228   }
00229 
00230   template <class Scalar> inline 
00231   LinearOperator<Scalar>
00232   block2x1(const LinearOperator<Scalar>& A00,
00233            const LinearOperator<Scalar>& A10)
00234   {
00235     return nonconstBlock2x1(A00.ptr(),A10.ptr());
00236   }
00237 
00238   template <class Scalar> inline 
00239   LinearOperator<Scalar>
00240   block1x2(const LinearOperator<Scalar>& A00,
00241            const LinearOperator<Scalar>& A01)
00242   {
00243     return nonconstBlock2x1(A00.ptr(),A01.ptr());
00244   }
00245   
00246 } // namespace Thyra
00247 
00248 #endif

Generated on Thu Sep 18 12:32:31 2008 for Thyra Operator/Vector Support by doxygen 1.3.9.1