# ExampleTridiagSerialLinearOp.hpp

Click here and here for example programs that use this example class.

```// @HEADER
// ***********************************************************************
//
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************

#ifndef THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
#define THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP

#include "Thyra_SpmdLinearOpBase.hpp"

template<class Scalar>
class ExampleTridiagSerialLinearOp : public Thyra::SpmdLinearOpBase<Scalar> {
private:

Thyra::Index         dim_;
std::vector<Scalar>  lower_;   // size = dim - 1
std::vector<Scalar>  diag_;    // size = dim
std::vector<Scalar>  upper_;   // size = dim - 1

public:

using Thyra::SpmdLinearOpBase<Scalar>::euclideanApply;

ExampleTridiagSerialLinearOp() : dim_(0) {}

ExampleTridiagSerialLinearOp( const Thyra::Index dim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
{ this->initialize(dim,lower,diag,upper); }

void initialize(
const Thyra::Index              dim      // >= 2
,const Scalar                   lower[]  // size == dim - 1
,const Scalar                   diag[]   // size == dim
,const Scalar                   upper[]  // size == dim - 1
)
{
TEST_FOR_EXCEPT( dim < 2 );
this->setLocalDimensions(Teuchos::null,dim,dim); // Needed to setup range() and domain()
dim_ = dim;
lower_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) lower_[k] = lower[k];
diag_.resize(dim);     for( int k = 0; k < dim;   ++k ) diag_[k]  = diag[k];
upper_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) upper_[k] = upper[k];
}

// Overridden form Teuchos::Describable */

std::string description() const
{
return (std::string("ExampleTridiagSerialLinearOp<") + Teuchos::ScalarTraits<Scalar>::name() + std::string(">"));
}

protected:

// Overridden from SingleScalarEuclideanLinearOpBase

bool opSupported(Thyra::ETransp M_trans) const { return true; }  // This class supports everything!

// Overridden from SpmdLinearOpBase

void euclideanApply(
const Thyra::ETransp                         M_trans
,const RTOpPack::ConstSubVectorView<Scalar>  &x_in
,const RTOpPack::SubVectorView<Scalar>       *y_out
,const Scalar                                alpha
,const Scalar                                beta
) const
{
typedef Teuchos::ScalarTraits<Scalar> ST;
// Get raw pointers to the values
const Scalar *x = x_in.values();
Scalar       *y = y_out->values();
// Perform y = beta*y (being careful to set y=0 if beta=0 in case y is uninitialized on input!)
Thyra::Index k = 0;
if( beta == ST::zero() ) {
for( k = 0; k < dim_; ++k ) y[k] = ST::zero();
}
else if( beta != ST::one() ) {
for( k = 0; k < dim_; ++k ) y[k] *= beta;
}
// Perform y = alpha*op(M)*x
k = 0;
if( M_trans == Thyra::NOTRANS ) {
y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );                         // First row
for( k = 1; k < dim_ - 1; ++k )
y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );  // Middle rows
y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] );                       // Last row
}
else if( M_trans == Thyra::CONJ ) {
y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
for( k = 1; k < dim_ - 1; ++k )
y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
}
else if( M_trans == Thyra::TRANS ) {
y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
for( k = 1; k < dim_ - 1; ++k )
y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
}
else if( M_trans == Thyra::CONJTRANS ) {
y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
for( k = 1; k < dim_ - 1; ++k )
y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
}
else {
TEST_FOR_EXCEPT(true); // Throw exception if we get here!
}
}

};  // end class ExampleTridiagSerialLinearOp

#endif  // THYRA_EXAMPLE_TRIDIAG_SERIAL_LINEAR_OP_HPP
```

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