Thyra Version of the Day
ExampleTridiagSerialLinearOp.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 // 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 #ifndef THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
00043 #define THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
00044 
00045 #include "Thyra_LinearOpDefaultBase.hpp"
00046 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00047 #include "Thyra_DetachedVectorView.hpp"
00048 #include "Teuchos_Assert.hpp"
00049 
00050 
00078 template<class Scalar>
00079 class ExampleTridiagSerialLinearOp : public Thyra::LinearOpDefaultBase<Scalar>
00080 {
00081 public:
00082 
00084   ExampleTridiagSerialLinearOp() {}
00085 
00087   ExampleTridiagSerialLinearOp(
00088     const Thyra::Ordinal dim,
00089     const Teuchos::ArrayView<const Scalar> &lower,
00090     const Teuchos::ArrayView<const Scalar> &diag,
00091     const Teuchos::ArrayView<const Scalar> &upper
00092     )
00093     { this->initialize(dim, lower, diag, upper);  }
00094   
00116   void initialize(
00117     const Thyra::Ordinal dim,
00118     const Teuchos::ArrayView<const Scalar> &lower,
00119     const Teuchos::ArrayView<const Scalar> &diag,
00120     const Teuchos::ArrayView<const Scalar> &upper
00121     )
00122     {
00123       TEUCHOS_TEST_FOR_EXCEPT( dim < 2 );
00124       space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim);
00125       lower_ = lower;
00126       diag_ = diag;
00127       upper_ = upper;
00128     }
00129 
00130 protected:
00131 
00132   // Overridden from LinearOpBase
00133 
00135   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > range() const
00136     { return space_; }
00137 
00139   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > domain() const
00140     { return space_; }
00141 
00143   bool opSupportedImpl(Thyra::EOpTransp M_trans) const
00144     { return true; }  // This class supports everything!
00145 
00147   void applyImpl(
00148     const Thyra::EOpTransp M_trans,
00149     const Thyra::MultiVectorBase<Scalar> &X_in,
00150     const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00151     const Scalar alpha,
00152     const Scalar beta
00153     ) const;
00154 
00155 private:
00156 
00157   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > space_;
00158   Teuchos::Array<Scalar> lower_;   // size = dim - 1    
00159   Teuchos::Array<Scalar> diag_;    // size = dim
00160   Teuchos::Array<Scalar> upper_;   // size = dim - 1
00161 
00162 };  // end class ExampleTridiagSerialLinearOp
00163 
00164 
00165 template<class Scalar>
00166 void ExampleTridiagSerialLinearOp<Scalar>::applyImpl(
00167   const Thyra::EOpTransp M_trans,
00168   const Thyra::MultiVectorBase<Scalar> &X_in,
00169   const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
00170   const Scalar alpha,
00171   const Scalar beta
00172   ) const
00173 {
00174 
00175   typedef Teuchos::ScalarTraits<Scalar> ST;
00176   typedef Thyra::Ordinal Ordinal;
00177 
00178   const Ordinal dim = space_->dim();
00179       
00180   // Loop over the input columns
00181 
00182   const Ordinal m = X_in.domain()->dim();
00183 
00184   for (Ordinal col_j = 0; col_j < m; ++col_j) {
00185 
00186     // Get access the the elements of column j
00187     Thyra::ConstDetachedVectorView<Scalar> x_vec(X_in.col(col_j));
00188     Thyra::DetachedVectorView<Scalar> y_vec(Y_inout->col(col_j));
00189     const Teuchos::ArrayRCP<const Scalar> x = x_vec.sv().values();
00190     const Teuchos::ArrayRCP<Scalar> y = y_vec.sv().values();
00191         
00192     // Perform y = beta*y (being careful to set y=0 if beta=0 in case y is
00193     // uninitialized on input!)
00194     if( beta == ST::zero() ) {
00195       for( Ordinal k = 0; k < dim; ++k ) y[k] = ST::zero();
00196     }
00197     else if( beta != ST::one() ) {
00198       for( Ordinal k = 0; k < dim; ++k ) y[k] *= beta;
00199     }
00200 
00201     // Perform y = alpha*op(M)*x 
00202     Ordinal k = 0;
00203     if( M_trans == Thyra::NOTRANS ) {
00204       y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );  // First row
00205       for( k = 1; k < dim - 1; ++k )   // Middle rows
00206         y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );
00207       y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] ); // Last row
00208     }
00209     else if( M_trans == Thyra::CONJ ) {
00210       y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
00211       for( k = 1; k < dim - 1; ++k )
00212         y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1]
00213           + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
00214       y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
00215     }
00216     else if( M_trans == Thyra::TRANS ) {
00217       y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
00218       for( k = 1; k < dim - 1; ++k )
00219         y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
00220       y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
00221     }
00222     else if( M_trans == Thyra::CONJTRANS ) {
00223       y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
00224       for( k = 1; k < dim - 1; ++k )
00225         y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1]
00226           + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
00227       y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
00228     }
00229     else {
00230       TEUCHOS_TEST_FOR_EXCEPT(true); // Throw exception if we get here!
00231     }
00232   }
00233 
00234 }
00235 
00236 
00237 #endif  // THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines