Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultColumnwiseMultiVector_def.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_DEFAULT_COLUMNWISE_MULTI_VECTOR_DEF_HPP
00030 #define THYRA_DEFAULT_COLUMNWISE_MULTI_VECTOR_DEF_HPP
00031 
00032 #include "Thyra_DefaultColumnwiseMultiVector_decl.hpp"
00033 #include "Thyra_MultiVectorDefaultBase.hpp"
00034 #include "Thyra_VectorSpaceBase.hpp"
00035 #include "Thyra_VectorBase.hpp"
00036 #include "Thyra_MultiVectorBase.hpp"
00037 #include "Thyra_VectorStdOps.hpp"
00038 #include "Thyra_VectorSpaceFactoryBase.hpp"
00039 #include "Thyra_AssertOp.hpp"
00040 #include "Teuchos_TestForException.hpp"
00041 #include "Teuchos_as.hpp"
00042 
00043 namespace Thyra {
00044 
00045 
00046 // Constructors/Initializers
00047 
00048 
00049 template<class Scalar>
00050 DefaultColumnwiseMultiVector<Scalar>::DefaultColumnwiseMultiVector()
00051 {}
00052 
00053 
00054 template<class Scalar>
00055 DefaultColumnwiseMultiVector<Scalar>::DefaultColumnwiseMultiVector(
00056   const RCP<VectorBase<Scalar> > &col_vec
00057   )
00058 {
00059   this->initialize(col_vec);
00060 }
00061 
00062 
00063 template<class Scalar>
00064 DefaultColumnwiseMultiVector<Scalar>::DefaultColumnwiseMultiVector(
00065   const RCP<const VectorSpaceBase<Scalar> > &range_in,
00066   const RCP<const VectorSpaceBase<Scalar> > &domain_in,
00067   const ArrayView<const RCP<VectorBase<Scalar> > > &col_vecs_in
00068   )
00069 {
00070   this->initialize(range_in, domain_in, col_vecs_in);
00071 }
00072 
00073 
00074 template<class Scalar>
00075 void DefaultColumnwiseMultiVector<Scalar>::initialize(
00076   const RCP<VectorBase<Scalar> > &col_vec
00077   )
00078 {
00079 #ifdef TEUCHOS_DEBUG
00080   const std::string err_msg =
00081     "DefaultColumnwiseMultiVector<Scalar>::initialize(...): Error!";
00082   TEST_FOR_EXCEPT_MSG( is_null(col_vec), err_msg ); 
00083   TEST_FOR_EXCEPT_MSG( is_null(col_vec->space()), err_msg ); 
00084 #endif
00085   range_  = col_vec->space();
00086   domain_ = range_->smallVecSpcFcty()->createVecSpc(1);
00087   col_vecs_.resize(1);
00088   col_vecs_[0] = col_vec;
00089 }
00090 
00091   
00092 template<class Scalar>
00093 void DefaultColumnwiseMultiVector<Scalar>::initialize(
00094   const RCP<const VectorSpaceBase<Scalar> > &range_in,
00095   const RCP<const VectorSpaceBase<Scalar> > &domain_in,
00096   const ArrayView<const RCP<VectorBase<Scalar> > > &col_vecs
00097   )
00098 {
00099 #ifdef TEUCHOS_DEBUG
00100   const std::string err_msg =
00101     "DefaultColumnwiseMultiVector<Scalar>::initialize(...): Error!";
00102   TEST_FOR_EXCEPT_MSG( is_null(range_in), err_msg ); 
00103   TEST_FOR_EXCEPT_MSG( is_null(domain_in), err_msg ); 
00104   TEST_FOR_EXCEPT_MSG( range_in->dim()  == 0, err_msg ); 
00105   TEST_FOR_EXCEPT_MSG( domain_in->dim() == 0, err_msg );
00106   // ToDo: Check the compatibility of the vectors in col_vecs!
00107 #endif
00108   const int domainDim = domain_in->dim();
00109   range_ = range_in;
00110   domain_ = domain_in;
00111   col_vecs_.clear();
00112   col_vecs_.reserve(domainDim);
00113   if (col_vecs.size()) {
00114     for( Ordinal j = 0; j < domainDim; ++j )
00115       col_vecs_.push_back(col_vecs[j]);
00116   }
00117   else {
00118     for( Ordinal j = 0; j < domainDim; ++j )
00119       col_vecs_.push_back(createMember(range_));
00120   }
00121 }
00122 
00123 
00124 template<class Scalar>
00125 void DefaultColumnwiseMultiVector<Scalar>::uninitialize()
00126 {
00127   col_vecs_.resize(0);
00128   range_ = Teuchos::null;
00129   domain_ = Teuchos::null;
00130 }
00131 
00132 
00133 // Overridden from OpBase
00134 
00135 
00136 template<class Scalar>
00137 RCP<const VectorSpaceBase<Scalar> >
00138 DefaultColumnwiseMultiVector<Scalar>::range() const
00139 {
00140   return range_;
00141 }
00142 
00143 
00144 template<class Scalar>
00145 RCP<const VectorSpaceBase<Scalar> >
00146 DefaultColumnwiseMultiVector<Scalar>::domain() const
00147 {
00148   return domain_;
00149 }
00150 
00151 
00152 // Overridden from LinearOpBase
00153 
00154 
00155 template<class Scalar>
00156 bool DefaultColumnwiseMultiVector<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00157 {
00158   typedef Teuchos::ScalarTraits<Scalar> ST;
00159   return ( ST::isComplex ? ( M_trans==NOTRANS || M_trans==CONJTRANS ) : true );
00160 }
00161 
00162 
00163 template<class Scalar>
00164 void DefaultColumnwiseMultiVector<Scalar>::applyImpl(
00165   const EOpTransp M_trans,
00166   const MultiVectorBase<Scalar> &X,
00167   const Ptr<MultiVectorBase<Scalar> > &Y,
00168   const Scalar alpha,
00169   const Scalar beta
00170   ) const
00171 {
00172 #ifdef TEUCHOS_DEBUG
00173   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00174     "MultiVectorBase<Scalar>::apply()", *this, M_trans, X, &*Y);
00175 #endif
00176   const Ordinal nc = this->domain()->dim();
00177   const Ordinal m = X.domain()->dim();
00178   for (Ordinal col_j = 0; col_j < m; ++col_j) {
00179     const RCP<const VectorBase<Scalar> > x_j = X.col(col_j);
00180     const RCP<VectorBase<Scalar> > y_j = Y->col(col_j);
00181     // y_j *= beta
00182     Vt_S(y_j.ptr(), beta);
00183     // y_j += alpha*op(M)*x_j
00184     if(M_trans == NOTRANS) {
00185       //
00186       // y_j += alpha*M*x_j = alpha*M.col(0)*x_j(0) + ... + alpha*M.col(nc-1)*x_j(nc-1)
00187       //
00188       // Extract an explicit view of x_j
00189       RTOpPack::ConstSubVectorView<Scalar> x_sub_vec;               
00190       x_j->acquireDetachedView(Range1D(), &x_sub_vec);
00191       // Loop through and add the multiple of each column
00192       for (Ordinal j = 0; j < nc; ++j )
00193         Vp_StV( y_j.ptr(), Scalar(alpha*x_sub_vec(j)), *this->col(j) );
00194       // Release the view of x
00195       x_j->releaseDetachedView(&x_sub_vec);
00196     }
00197     else {
00198       //
00199       //                        [ alpha*dot(M.col(0),x_j)    ]
00200       // y_j += alpha*M^T*x_j = [ alpha*dot(M.col(1),x_j)    ]
00201       //                        [ ...                        ]
00202       //                        [ alpha*dot(M.col(nc-1),x_j) ]
00203       //
00204       // Extract an explicit view of y_j
00205       RTOpPack::SubVectorView<Scalar> y_sub_vec;               
00206       y_j->acquireDetachedView(Range1D(), &y_sub_vec);
00207       // Loop through and add to each element in y_j
00208       for (Ordinal j = 0; j < nc; ++j )
00209         y_sub_vec(j) += alpha*dot(*this->col(j), *x_j);
00210       // Commit explicit view of y
00211       y_j->commitDetachedView(&y_sub_vec);
00212     }
00213   }
00214 }
00215 
00216 
00217 // Overridden from MultiVectorBase
00218 
00219 
00220 template<class Scalar>
00221 RCP<VectorBase<Scalar> >
00222 DefaultColumnwiseMultiVector<Scalar>::nonconstColImpl(Ordinal j)
00223 {
00224   using Teuchos::as;
00225   const int num_cols = col_vecs_.size();
00226   TEST_FOR_EXCEPTION(
00227     !(  0 <= j  && j < num_cols ), std::logic_error
00228     ,"Error, j = " << j << " does not fall in the range [0,"<<(num_cols-1)<< "]!"
00229     );
00230   return col_vecs_[j];
00231 }
00232 
00233 
00234 template<class Scalar>
00235 RCP<MultiVectorBase<Scalar> >
00236 DefaultColumnwiseMultiVector<Scalar>::nonconstContigSubViewImpl(
00237   const Range1D& col_rng_in
00238   )
00239 {
00240   const Ordinal numCols = domain_->dim();
00241   const Range1D col_rng = Teuchos::full_range(col_rng_in,0,numCols-1);
00242 #ifdef TEUCHOS_DEBUG
00243   TEST_FOR_EXCEPTION(
00244     !( col_rng.ubound() < numCols ), std::logic_error
00245     ,"DefaultColumnwiseMultiVector<Scalar>::subView(col_rng):"
00246     "Error, the input range col_rng = ["<<col_rng.lbound()
00247     <<","<<col_rng.ubound()<<"] "
00248     "is not in the range [0,"<<(numCols-1)<<"]!"
00249     );
00250 #endif
00251   return Teuchos::rcp(
00252     new DefaultColumnwiseMultiVector<Scalar>(
00253       range_,
00254       domain_->smallVecSpcFcty()->createVecSpc(col_rng.size()),
00255       col_vecs_(col_rng.lbound(),col_rng.size())
00256       )
00257     );
00258 }
00259 
00260   
00261 } // end namespace Thyra
00262 
00263 
00264 #endif // THYRA_DEFAULT_COLUMNWISE_MULTI_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines