Thyra Package Browser (Single Doxygen Collection) Version of the Day
ExampleTridiagSerialLinearOp.hpp
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 // 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_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
00030 #define THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
00031 
00032 #include "Thyra_LinearOpDefaultBase.hpp"
00033 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00034 #include "Thyra_DetachedVectorView.hpp"
00035 #include "Teuchos_Assert.hpp"
00036 
00037 
00065 template<class Scalar>
00066 class ExampleTridiagSerialLinearOp : public Thyra::LinearOpDefaultBase<Scalar>
00067 {
00068 public:
00069 
00071   ExampleTridiagSerialLinearOp() {}
00072 
00074   ExampleTridiagSerialLinearOp(
00075     const Thyra::Ordinal dim,
00076     const Teuchos::ArrayView<const Scalar> &lower,
00077     const Teuchos::ArrayView<const Scalar> &diag,
00078     const Teuchos::ArrayView<const Scalar> &upper
00079     )
00080     { this->initialize(dim, lower, diag, upper);  }
00081   
00103   void initialize(
00104     const Thyra::Ordinal dim,
00105     const Teuchos::ArrayView<const Scalar> &lower,
00106     const Teuchos::ArrayView<const Scalar> &diag,
00107     const Teuchos::ArrayView<const Scalar> &upper
00108     )
00109     {
00110       TEST_FOR_EXCEPT( dim < 2 );
00111       space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim);
00112       lower_ = lower;
00113       diag_ = diag;
00114       upper_ = upper;
00115     }
00116 
00117 protected:
00118 
00119   // Overridden from LinearOpBase
00120 
00122   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > range() const
00123     { return space_; }
00124 
00126   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > domain() const
00127     { return space_; }
00128 
00130   bool opSupportedImpl(Thyra::EOpTransp M_trans) const
00131     { return true; }  // This class supports everything!
00132 
00134   void applyImpl(
00135     const Thyra::EOpTransp M_trans,
00136     const Thyra::MultiVectorBase<Scalar> &X_in,
00137     const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00138     const Scalar alpha,
00139     const Scalar beta
00140     ) const;
00141 
00142 private:
00143 
00144   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > space_;
00145   Teuchos::Array<Scalar> lower_;   // size = dim - 1    
00146   Teuchos::Array<Scalar> diag_;    // size = dim
00147   Teuchos::Array<Scalar> upper_;   // size = dim - 1
00148 
00149 };  // end class ExampleTridiagSerialLinearOp
00150 
00151 
00152 template<class Scalar>
00153 void ExampleTridiagSerialLinearOp<Scalar>::applyImpl(
00154   const Thyra::EOpTransp M_trans,
00155   const Thyra::MultiVectorBase<Scalar> &X_in,
00156   const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00157   const Scalar alpha,
00158   const Scalar beta
00159   ) const
00160 {
00161 
00162   typedef Teuchos::ScalarTraits<Scalar> ST;
00163   typedef Thyra::Ordinal Ordinal;
00164 
00165   const Ordinal dim = space_->dim();
00166       
00167   // Loop over the input columns
00168 
00169   const Ordinal m = X_in.domain()->dim();
00170 
00171   for (Ordinal col_j = 0; col_j < m; ++col_j) {
00172 
00173     // Get access the the elements of column j
00174     Thyra::ConstDetachedVectorView<Scalar> x_vec(X_in.col(col_j));
00175     Thyra::DetachedVectorView<Scalar> y_vec(Y_inout->col(col_j));
00176     const Teuchos::ArrayRCP<const Scalar> x = x_vec.sv().values();
00177     const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
00178         
00179     // Perform y = beta*y (being careful to set y=0 if beta=0 in case y is
00180     // uninitialized on input!)
00181     if( beta == ST::zero() ) {
00182       for( Ordinal k = 0; k < dim; ++k ) y[k] = ST::zero();
00183     }
00184     else if( beta != ST::one() ) {
00185       for( Ordinal k = 0; k < dim; ++k ) y[k] *= beta;
00186     }
00187 
00188     // Perform y = alpha*op(M)*x 
00189     Ordinal k = 0;
00190     if( M_trans == Thyra::NOTRANS ) {
00191       y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );  // First row
00192       for( k = 1; k < dim - 1; ++k )   // Middle rows
00193         y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );
00194       y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] ); // Last row
00195     }
00196     else if( M_trans == Thyra::CONJ ) {
00197       y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
00198       for( k = 1; k < dim - 1; ++k )
00199         y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1]
00200           + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
00201       y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
00202     }
00203     else if( M_trans == Thyra::TRANS ) {
00204       y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
00205       for( k = 1; k < dim - 1; ++k )
00206         y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
00207       y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
00208     }
00209     else if( M_trans == Thyra::CONJTRANS ) {
00210       y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
00211       for( k = 1; k < dim - 1; ++k )
00212         y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1]
00213           + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
00214       y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
00215     }
00216     else {
00217       TEST_FOR_EXCEPT(true); // Throw exception if we get here!
00218     }
00219   }
00220 
00221 }
00222 
00223 
00224 #endif  // THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines