Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_MultiVectorDefaultBase_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_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
00030 #define THYRA_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
00031 
00032 
00033 #include "Thyra_MultiVectorDefaultBase_decl.hpp"
00034 #include "Thyra_LinearOpDefaultBase.hpp"
00035 #include "Thyra_MultiVectorStdOps.hpp"
00036 #include "Thyra_VectorSpaceFactoryBase.hpp"
00037 #include "Thyra_VectorSpaceBase.hpp"
00038 #include "Thyra_VectorBase.hpp"
00039 #include "Thyra_AssertOp.hpp"
00040 #include "Thyra_DefaultColumnwiseMultiVector.hpp"
00041 #include "Teuchos_Workspace.hpp"
00042 #include "Teuchos_TestForException.hpp"
00043 #include "Teuchos_as.hpp"
00044 
00045 
00046 namespace Thyra {
00047 
00048 
00049 // Overridden public member functions from MultiVectorBase
00050 
00051 
00052 template<class Scalar>
00053 RCP<MultiVectorBase<Scalar> >
00054 MultiVectorDefaultBase<Scalar>::clone_mv() const
00055 {
00056   const VectorSpaceBase<Scalar>
00057     &l_domain = *this->domain(),
00058     &l_range = *this->range();
00059   RCP<MultiVectorBase<Scalar> >
00060     copy = createMembers(l_range,l_domain.dim());
00061   assign( copy.ptr(), *this );
00062   return copy;
00063 }
00064 
00065 
00066 // protected
00067 
00068 
00069 // Overridden protected member functions from MultiVectorBase
00070 
00071 
00072 template<class Scalar>
00073 RCP<const MultiVectorBase<Scalar> >
00074 MultiVectorDefaultBase<Scalar>::contigSubViewImpl( const Range1D& colRng_in ) const
00075 {
00076   using Teuchos::Workspace;
00077   using Teuchos::as;
00078   Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
00079   const VectorSpaceBase<Scalar> &l_domain = *this->domain();
00080   const VectorSpaceBase<Scalar> &l_range = *this->range();
00081   const Ordinal dimDomain = l_domain.dim();
00082   const Range1D colRng = Teuchos::full_range(colRng_in,0,dimDomain-1);
00083   if( colRng.lbound() == 0 && as<Ordinal>(colRng.ubound()) == dimDomain-1 )
00084     return Teuchos::rcp(this,false); // Takes all of the columns!
00085   if( colRng.size() ) {
00086     // We have to create a view of a subset of the columns
00087     Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,colRng.size());
00088     for( Ordinal j = colRng.lbound(); j <= colRng.ubound(); ++j )
00089       col_vecs[j-colRng.lbound()] = Teuchos::rcp_const_cast<VectorBase<Scalar> >(this->col(j));
00090     return Teuchos::rcp(
00091       new DefaultColumnwiseMultiVector<Scalar>(
00092         this->range(),l_range.smallVecSpcFcty()->createVecSpc(colRng.size()),col_vecs
00093         )
00094       );
00095   }
00096   return Teuchos::null; // There was an empty set in colRng_in!
00097 }
00098 
00099 
00100 template<class Scalar>
00101 RCP<MultiVectorBase<Scalar> >
00102 MultiVectorDefaultBase<Scalar>::nonconstContigSubViewImpl( const Range1D& colRng_in )
00103 {
00104   using Teuchos::Workspace;
00105   using Teuchos::as;
00106   Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
00107   const VectorSpaceBase<Scalar> &l_domain = *this->domain();
00108   const VectorSpaceBase<Scalar> &l_range = *this->range();
00109   const Ordinal dimDomain = l_domain.dim();
00110   const Range1D colRng = Teuchos::full_range(colRng_in,0,dimDomain-1);
00111   if( colRng.lbound() == 0 && as<Ordinal>(colRng.ubound()) == dimDomain-1 )
00112     return Teuchos::rcp(this,false); // Takes all of the columns!
00113   if( colRng.size() ) {
00114     // We have to create a view of a subset of the columns
00115     Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,colRng.size());
00116     for( Ordinal j = colRng.lbound(); j <= colRng.ubound(); ++j )
00117       col_vecs[j-colRng.lbound()] = this->col(j);
00118     return Teuchos::rcp(
00119       new DefaultColumnwiseMultiVector<Scalar>(
00120         this->range(),l_range.smallVecSpcFcty()->createVecSpc(colRng.size()),col_vecs
00121         )
00122       );
00123   }
00124   return Teuchos::null; // There was an empty set in colRng_in!
00125 }
00126 
00127 
00128 template<class Scalar>
00129 RCP<const MultiVectorBase<Scalar> >
00130 MultiVectorDefaultBase<Scalar>::nonContigSubViewImpl(
00131   const ArrayView<const int> &cols
00132   ) const
00133 {
00134   using Teuchos::Workspace;
00135   Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
00136   const VectorSpaceBase<Scalar> &l_range = *this->range();
00137   const int numCols = cols.size();
00138 #ifdef TEUCHOS_DEBUG
00139   const VectorSpaceBase<Scalar> &l_domain = *this->domain();
00140   const Ordinal dimDomain = l_domain.dim();
00141   const char msg_err[] = "MultiVectorDefaultBase<Scalar>::subView(numCols,cols[]): Error!";
00142   TEST_FOR_EXCEPTION( numCols < 1 || dimDomain < numCols, std::invalid_argument, msg_err );
00143 #endif
00144   // We have to create a view of a subset of the columns
00145   Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,numCols);
00146   for( int k = 0; k < numCols; ++k ) {
00147     const int col_k = cols[k];
00148 #ifdef TEUCHOS_DEBUG
00149     TEST_FOR_EXCEPTION(
00150       !( 0 <= col_k && col_k < dimDomain ), std::invalid_argument
00151       ,msg_err << " col["<<k<<"] = " << col_k << " is not in the range [0,"<<(dimDomain-1)<<"]!"
00152       );
00153 #endif
00154     col_vecs[k] = Teuchos::rcp_const_cast<VectorBase<Scalar> >(this->col(col_k));
00155   }
00156   return Teuchos::rcp(
00157     new DefaultColumnwiseMultiVector<Scalar>(
00158       this->range(), l_range.smallVecSpcFcty()->createVecSpc(numCols), col_vecs
00159       )
00160     );
00161 }
00162 
00163 
00164 template<class Scalar>
00165 RCP<MultiVectorBase<Scalar> >
00166 MultiVectorDefaultBase<Scalar>::nonconstNonContigSubViewImpl(
00167   const ArrayView<const int> &cols
00168   )
00169 {
00170   using Teuchos::Workspace;
00171   Teuchos::WorkspaceStore *wss = Teuchos::get_default_workspace_store().get();
00172   const VectorSpaceBase<Scalar> &l_range = *this->range();
00173   const int numCols = cols.size();
00174 #ifdef TEUCHOS_DEBUG
00175   const VectorSpaceBase<Scalar> &l_domain = *this->domain();
00176   const Ordinal dimDomain = l_domain.dim();
00177   const char msg_err[] = "MultiVectorDefaultBase<Scalar>::subView(numCols,cols[]): Error!";
00178   TEST_FOR_EXCEPTION( numCols < 1 || dimDomain < numCols, std::invalid_argument, msg_err );
00179 #endif
00180   // We have to create a view of a subset of the columns
00181   Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,numCols);
00182   for( int k = 0; k < numCols; ++k ) {
00183     const int col_k = cols[k];
00184 #ifdef TEUCHOS_DEBUG
00185     TEST_FOR_EXCEPTION(
00186       !( 0 <= col_k && col_k < dimDomain ), std::invalid_argument
00187       ,msg_err << " col["<<k<<"] = " << col_k << " is not in the range [0,"<<(dimDomain-1)<<"]!"
00188       );
00189 #endif
00190     col_vecs[k] = this->col(col_k);
00191   }
00192   return Teuchos::rcp(
00193     new DefaultColumnwiseMultiVector<Scalar>(
00194       this->range(), l_range.smallVecSpcFcty()->createVecSpc(numCols), col_vecs
00195       )
00196     );
00197 }
00198 
00199 
00200 template<class Scalar>
00201 void MultiVectorDefaultBase<Scalar>::mvMultiReductApplyOpImpl(
00202   const RTOpPack::RTOpT<Scalar> &prim_op,
00203   const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs,
00204   const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs,
00205   const ArrayView<const Ptr<RTOpPack::ReductTarget> > &reduct_objs,
00206   const Ordinal prim_global_offset_in
00207   ) const
00208 {
00209 
00210   using Teuchos::Workspace;
00211   using Teuchos::as;
00212   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00213 
00214   const int num_multi_vecs = multi_vecs.size();
00215   const int num_targ_multi_vecs = targ_multi_vecs.size();
00216 
00217   // ToDo: Validate the input!
00218 
00219   const VectorSpaceBase<Scalar> &l_domain = *this->domain();
00220 
00221   // Get the primary and secondary dimensions.
00222 
00223   const Ordinal sec_dim = l_domain.dim();
00224 
00225   //
00226   // Apply the reduction/transformation operator and transform the
00227   // target vectors and reduce each of the reduction objects.
00228   //
00229 
00230   Workspace<RCP<const VectorBase<Scalar> > > vecs_s(wss, num_multi_vecs);
00231   Workspace<Ptr<const VectorBase<Scalar> > > vecs(wss, num_multi_vecs);
00232   Workspace<RCP<VectorBase<Scalar> > > targ_vecs_s(wss, num_targ_multi_vecs);
00233   Workspace<Ptr<VectorBase<Scalar> > > targ_vecs(wss, num_targ_multi_vecs);
00234 
00235   for(Ordinal j = 0; j < sec_dim; ++j) {
00236     // Fill the arrays of vector arguments
00237     {for(Ordinal k = 0; k < as<Ordinal>(num_multi_vecs); ++k) {
00238         vecs_s[k] = multi_vecs[k]->col(j);
00239         vecs[k] = vecs_s[k].ptr();
00240       }}
00241     {for(Ordinal k = 0; k < as<Ordinal>(num_targ_multi_vecs); ++k) {
00242         targ_vecs_s[k] = targ_multi_vecs[k]->col(j);
00243         targ_vecs[k] = targ_vecs_s[k].ptr();
00244       }}
00245     // Apply the reduction/transformation operator
00246     Thyra::applyOp(
00247       prim_op,
00248       vecs().getConst(),
00249       targ_vecs().getConst(),
00250       reduct_objs.size() ? reduct_objs[j] : Ptr<RTOpPack::ReductTarget>(),
00251       prim_global_offset_in);
00252   }
00253   // At this point all of the designated targ vectors in the target multi-vectors have
00254   // been transformed and all the reduction objects in reduct_obj[] have accumulated
00255   // the reductions.
00256 }
00257 
00258 
00259 template<class Scalar>
00260 void MultiVectorDefaultBase<Scalar>::mvSingleReductApplyOpImpl(
00261   const RTOpPack::RTOpT<Scalar> &prim_op,
00262   const RTOpPack::RTOpT<Scalar> &sec_op,
00263   const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs,
00264   const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs,
00265   const Ptr<RTOpPack::ReductTarget> &reduct_obj,
00266   const Ordinal prim_global_offset_in
00267   ) const
00268 {
00269 
00270   using Teuchos::Workspace;
00271   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00272 
00273   // ToDo: Validate the input!
00274 
00275   const VectorSpaceBase<Scalar> &l_domain = *this->domain();
00276 
00277   // Get the primary and secondary dimensions.
00278   const Ordinal sec_dim = l_domain.dim();
00279 
00280   // Create a temporary buffer for the reduction objects of the primary reduction
00281   // so that we can call the companion version of this method.
00282   const int reduct_objs_size = (!is_null(reduct_obj) ? sec_dim : 0);
00283   Workspace<RCP<RTOpPack::ReductTarget> > rcp_reduct_objs(wss, reduct_objs_size);
00284   Workspace<Ptr<RTOpPack::ReductTarget> > reduct_objs(wss, reduct_objs_size);
00285   if (!is_null(reduct_obj)) {
00286     for(Ordinal k = 0; k < sec_dim; ++k) {
00287       rcp_reduct_objs[k] = prim_op.reduct_obj_create();
00288       reduct_objs[k] = rcp_reduct_objs[k].ptr();
00289     }
00290   }
00291  
00292   // Call the companion version that accepts an array of reduction objects
00293   this->applyOp(
00294     prim_op, multi_vecs, targ_multi_vecs, reduct_objs,
00295     prim_global_offset_in);
00296 
00297   // Reduce all the reduction objects using the secondary reduction operator
00298   // into one reduction object and free the intermediate reduction objects.
00299   if (!is_null(reduct_obj)) {
00300     for (Ordinal k = 0; k < sec_dim; ++k) {
00301       sec_op.reduce_reduct_objs( *reduct_objs[k], reduct_obj );
00302     }
00303   }
00304 }
00305 
00306 
00307 template<class Scalar>
00308 void MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(
00309   const Range1D &rowRng_in,
00310   const Range1D &colRng_in,
00311   RTOpPack::ConstSubMultiVectorView<Scalar> *sub_mv
00312   ) const
00313 {
00314   const Ordinal
00315     rangeDim = this->range()->dim(),
00316     domainDim = this->domain()->dim();
00317   const Range1D
00318     rowRng = rowRng_in.full_range() ? Range1D(0,rangeDim-1) : rowRng_in,
00319     colRng = colRng_in.full_range() ? Range1D(0,domainDim-1) : colRng_in;
00320 #ifdef TEUCHOS_DEBUG
00321   TEST_FOR_EXCEPTION(
00322     !(rowRng.ubound() < rangeDim), std::out_of_range
00323     ,"MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(...): Error, rowRng = ["
00324     <<rowRng.lbound()<<","<<rowRng.ubound()<<"] is not in the range = [0,"
00325     <<(rangeDim-1)<<"]!"
00326     );
00327   TEST_FOR_EXCEPTION(
00328     !(colRng.ubound() < domainDim), std::out_of_range
00329     ,"MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(...): Error, colRng = ["
00330     <<colRng.lbound()<<","<<colRng.ubound()<<"] is not in the range = [0,"
00331     <<(domainDim-1)<<"]!"
00332     );
00333 #endif
00334   // Allocate storage for the multi-vector (stored column-major)
00335   const ArrayRCP<Scalar> values = Teuchos::arcp<Scalar>(rowRng.size() * colRng.size());
00336   // Extract multi-vector values column by column
00337   RTOpPack::ConstSubVectorView<Scalar> sv; // uninitialized by default
00338   for( int k = colRng.lbound(); k <= colRng.ubound(); ++k ) {
00339     RCP<const VectorBase<Scalar> > col_k = this->col(k);
00340     col_k->acquireDetachedView( rowRng, &sv );
00341     for( int i = 0; i < rowRng.size(); ++i )
00342       values[ i + k*rowRng.size() ] = sv[i];
00343     col_k->releaseDetachedView( &sv );
00344   }
00345   // Initialize the multi-vector view object
00346   sub_mv->initialize(
00347     rowRng.lbound(), // globalOffset
00348     rowRng.size(), // subDim
00349     colRng.lbound(), // colOffset
00350     colRng.size(), // numSubCols
00351     values, // values
00352     rowRng.size() // leadingDim
00353     );
00354 }
00355 
00356 
00357 template<class Scalar>
00358 void MultiVectorDefaultBase<Scalar>::releaseDetachedMultiVectorViewImpl(
00359   RTOpPack::ConstSubMultiVectorView<Scalar>* sub_mv
00360   ) const
00361 {
00362   // Here we just need to free the view and that is it!
00363   sub_mv->uninitialize();
00364 }
00365 
00366 
00367 template<class Scalar>
00368 void MultiVectorDefaultBase<Scalar>::acquireNonconstDetachedMultiVectorViewImpl(
00369   const Range1D &rowRng,
00370   const Range1D &colRng,
00371   RTOpPack::SubMultiVectorView<Scalar> *sub_mv
00372   )
00373 {
00374   using Teuchos::as;
00375   // Use the non-const implementation since it does exactly the
00376   // correct thing in this case also!
00377   MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(
00378     rowRng, colRng,
00379     as<RTOpPack::ConstSubMultiVectorView<Scalar>*>(sub_mv)
00380     // This cast will work as long as SubMultiVectorView
00381     // maintains no extra state over ConstSubMultiVectorView (which it
00382     // currently does not) but this is something that I should
00383     // technically check for some how.
00384     );
00385 }
00386 
00387 
00388 template<class Scalar>
00389 void MultiVectorDefaultBase<Scalar>::commitNonconstDetachedMultiVectorViewImpl(
00390   RTOpPack::SubMultiVectorView<Scalar>* sub_mv
00391   )
00392 {
00393 #ifdef TEUCHOS_DEBUG
00394   TEST_FOR_EXCEPTION(
00395     sub_mv==NULL, std::logic_error,
00396     "MultiVectorDefaultBase<Scalar>::commitNonconstDetachedMultiVectorViewImpl(...): Error!"
00397     );
00398 #endif
00399   // Set back the multi-vector values column by column
00400   const Range1D rowRng(sub_mv->globalOffset(),sub_mv->globalOffset()+sub_mv->subDim()-1);
00401   RTOpPack::SubVectorView<Scalar> msv; // uninitialized by default
00402   for( int k = sub_mv->colOffset(); k < sub_mv->numSubCols(); ++k ) {
00403     RCP<VectorBase<Scalar> > col_k = this->col(k);
00404     col_k->acquireDetachedView( rowRng, &msv );
00405     for( int i = 0; i < rowRng.size(); ++i )
00406       msv[i] = sub_mv->values()[ i + k*rowRng.size() ];
00407     col_k->commitDetachedView( &msv );
00408   }
00409   // Zero out the view
00410   sub_mv->uninitialize();
00411 }
00412 
00413 
00414 } // end namespace Thyra
00415 
00416 
00417 #endif // THYRA_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines