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