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

Generated on Tue Oct 20 12:46:58 2009 for Thyra Operator/Vector Support by doxygen 1.4.7