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_SpmdLinearOpBase.hpp"
00033 
00068 template<class Scalar>
00069 class ExampleTridiagSerialLinearOp : public Thyra::SpmdLinearOpBase<Scalar> {
00070 private:
00071 
00072   Thyra::Index         dim_;
00073   std::vector<Scalar>  lower_;   // size = dim - 1    
00074   std::vector<Scalar>  diag_;    // size = dim
00075   std::vector<Scalar>  upper_;   // size = dim - 1
00076 
00077 public:
00078 
00080   using Thyra::SpmdLinearOpBase<Scalar>::euclideanApply;
00081 
00083   ExampleTridiagSerialLinearOp() : dim_(0) {}
00084 
00086   ExampleTridiagSerialLinearOp( const Thyra::Index dim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
00087     { this->initialize(dim,lower,diag,upper); }
00088   
00105   void initialize(
00106     const Thyra::Index              dim      // >= 2
00107     ,const Scalar                   lower[]  // size == dim - 1
00108     ,const Scalar                   diag[]   // size == dim
00109     ,const Scalar                   upper[]  // size == dim - 1
00110     )
00111     {
00112       TEST_FOR_EXCEPT( dim < 2 );
00113       this->setLocalDimensions(Teuchos::null,dim,dim); // Needed to setup range() and domain()
00114       dim_ = dim;
00115       lower_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) lower_[k] = lower[k];
00116       diag_.resize(dim);     for( int k = 0; k < dim;   ++k ) diag_[k]  = diag[k];
00117       upper_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) upper_[k] = upper[k];
00118     }
00119 
00120   // Overridden form Teuchos::Describable */
00121 
00122   std::string description() const
00123     {
00124       return (std::string("ExampleTridiagSerialLinearOp<") + Teuchos::ScalarTraits<Scalar>::name() + std::string(">"));
00125     }
00126 
00127 protected:
00128 
00129   // Overridden from SingleScalarEuclideanLinearOpBase
00130 
00131   bool opSupported(Thyra::ETransp M_trans) const { return true; }  // This class supports everything!
00132 
00133   // Overridden from SpmdLinearOpBase
00134 
00135   void euclideanApply(
00136     const Thyra::ETransp                         M_trans
00137     ,const RTOpPack::ConstSubVectorView<Scalar>  &x_in
00138     ,const RTOpPack::SubVectorView<Scalar>       *y_out
00139     ,const Scalar                                alpha
00140     ,const Scalar                                beta
00141     ) const
00142     {
00143       typedef Teuchos::ScalarTraits<Scalar> ST;
00144       // Get raw pointers to the values
00145       const Scalar *x = x_in.values();
00146       Scalar       *y = y_out->values();
00147       // Perform y = beta*y (being careful to set y=0 if beta=0 in case y is uninitialized on input!)
00148       Thyra::Index k = 0;
00149       if( beta == ST::zero() ) {
00150         for( k = 0; k < dim_; ++k ) y[k] = ST::zero();
00151       }
00152       else if( beta != ST::one() ) {
00153         for( k = 0; k < dim_; ++k ) y[k] *= beta;
00154       }
00155       // Perform y = alpha*op(M)*x 
00156       k = 0;
00157       if( M_trans == Thyra::NOTRANS ) {
00158         y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );                         // First row
00159         for( k = 1; k < dim_ - 1; ++k )
00160           y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );  // Middle rows
00161         y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] );                       // Last row
00162       }
00163       else if( M_trans == Thyra::CONJ ) {
00164         y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
00165         for( k = 1; k < dim_ - 1; ++k )
00166           y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
00167         y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
00168       }
00169       else if( M_trans == Thyra::TRANS ) {
00170         y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
00171         for( k = 1; k < dim_ - 1; ++k )
00172           y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
00173         y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
00174       }
00175       else if( M_trans == Thyra::CONJTRANS ) {
00176         y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
00177         for( k = 1; k < dim_ - 1; ++k )
00178           y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
00179         y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
00180       }
00181       else {
00182         TEST_FOR_EXCEPT(true); // Throw exception if we get here!
00183       }
00184     }
00185 
00186 };  // end class ExampleTridiagSerialLinearOp
00187 
00188 #endif  // THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP

Generated on Thu Sep 18 12:33:01 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1