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

Generated on Tue Oct 20 12:47:27 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7