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

Generated on Thu Sep 18 12:33:03 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1