Thyra_apply_op_helper.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 // 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_APPLY_OP_HELPER_HPP
00030 #define THYRA_APPLY_OP_HELPER_HPP
00031 
00032 #include "Thyra_apply_op_helper_decl.hpp"
00033 #include "Thyra_VectorBase.hpp"
00034 #include "Thyra_VectorSpaceBase.hpp"
00035 #include "Thyra_AssertOp.hpp"
00036 #include "Teuchos_Workspace.hpp"
00037 #include "Teuchos_TestForException.hpp"
00038 
00039 template<class Scalar>
00040 void Thyra::apply_op_validate_input(
00041   const char                      func_name[]
00042   ,const VectorSpaceBase<Scalar>  &space
00043   ,const RTOpPack::RTOpT<Scalar>  &op
00044   ,const int                      num_vecs
00045   ,const VectorBase<Scalar>*const vecs[]
00046   ,const int                      num_targ_vecs
00047   ,VectorBase<Scalar>*const       targ_vecs[]
00048   ,RTOpPack::ReductTarget         *reduct_obj
00049   ,const Index                    first_ele_offset_in
00050   ,const Index                    sub_dim_in
00051   ,const Index                    global_offset_in
00052   )
00053 {
00054   int k;
00055   const Index
00056     dim = space.dim();
00057   TEST_FOR_EXCEPTION(
00058     global_offset_in < 0, std::logic_error
00059     ,func_name << " : Error!  global_offset_in = "
00060     <<global_offset_in<<" is not valid" );
00061   TEST_FOR_EXCEPTION(
00062     first_ele_offset_in+1 > dim, std::logic_error
00063     ,func_name << " : Error!  first_ele_offset_in = "
00064     <<first_ele_offset_in<<" is not compatible with space.dim() = " << dim );
00065   TEST_FOR_EXCEPTION(
00066     (sub_dim_in > 0 && sub_dim_in > dim-first_ele_offset_in), std::logic_error
00067     ,func_name << " : Error!  first_ele_offset_in = "
00068     <<first_ele_offset_in<<" and sub_dim_in = "<<sub_dim_in
00069     <<" is not compatible with space.dim() = " << dim );
00070   for(k = 0; k < num_vecs; ++k)
00071     THYRA_ASSERT_VEC_SPACES(func_name,space,*vecs[k]->space());
00072   for(k = 0; k < num_targ_vecs; ++k)
00073     THYRA_ASSERT_VEC_SPACES(func_name,space,*targ_vecs[k]->space());
00074 }
00075 
00076 template<class Scalar>
00077 void Thyra::apply_op_validate_input(
00078   const char                            func_name[]
00079   ,const VectorSpaceBase<Scalar>        &domain
00080   ,const VectorSpaceBase<Scalar>        &range
00081   ,const RTOpPack::RTOpT<Scalar>        &primary_op
00082   ,const int                            num_multi_vecs
00083   ,const MultiVectorBase<Scalar>*const  multi_vecs[]
00084   ,const int                            num_targ_multi_vecs
00085   ,MultiVectorBase<Scalar>*const        targ_multi_vecs[]
00086   ,RTOpPack::ReductTarget*const         reduct_objs[]
00087   ,const Index                          primary_first_ele_offset_in
00088   ,const Index                          primary_sub_dim_in
00089   ,const Index                          primary_global_offset_in
00090   ,const Index                          secondary_first_ele_offset_in
00091   ,const Index                          secondary_sub_dim_in
00092   )
00093 {
00094   int k;
00095   // Validate primary range arguments
00096   const Index
00097     range_dim = range.dim();
00098   TEST_FOR_EXCEPTION(
00099     primary_global_offset_in < 0, std::logic_error
00100     ,func_name << " : Error!  primary_global_offset_in = "
00101     <<primary_global_offset_in<<" is not valid" );
00102   TEST_FOR_EXCEPTION(
00103     primary_first_ele_offset_in < 0 || range_dim < primary_first_ele_offset_in+1, std::logic_error
00104     ,func_name << " : Error!  primary_first_ele_offset_in = "
00105     <<primary_first_ele_offset_in<<" is not compatible with range.dim() = " << range_dim );
00106   TEST_FOR_EXCEPTION(
00107     (primary_sub_dim_in > 0 && primary_sub_dim_in > range_dim-primary_first_ele_offset_in)
00108     , std::logic_error
00109     ,func_name << " : Error!  primary_first_ele_offset_in = "
00110     <<primary_first_ele_offset_in<<" and primary_sub_dim_in = "<<primary_sub_dim_in
00111     <<" are not compatible with range.dim() = " << range_dim );
00112   // Validate secondary domain arguments
00113   const Index
00114     domain_dim = domain.dim();
00115   TEST_FOR_EXCEPTION(
00116     secondary_first_ele_offset_in < 0 || domain_dim < secondary_first_ele_offset_in+1, std::logic_error
00117     ,func_name << " : Error!  secondary_first_ele_offset_in = "
00118     <<secondary_first_ele_offset_in<<" is not compatible with domain.dim() = " << domain_dim );
00119   TEST_FOR_EXCEPTION(
00120     (secondary_sub_dim_in > 0 && secondary_sub_dim_in > domain_dim-secondary_first_ele_offset_in)
00121     , std::logic_error
00122     ,func_name << " : Error!  secondary_first_ele_offset_in = "
00123     <<secondary_first_ele_offset_in<<" and secondary_sub_dim_in = "<<secondary_sub_dim_in
00124     <<" are not compatible with domain.dim() = " << domain_dim );
00125   // Validate spaces
00126   for(k = 0; k < num_multi_vecs; ++k) {
00127     THYRA_ASSERT_VEC_SPACES(func_name,domain,*multi_vecs[k]->domain());
00128     THYRA_ASSERT_VEC_SPACES(func_name,range,*multi_vecs[k]->range());
00129   }
00130   for(k = 0; k < num_targ_multi_vecs; ++k) {
00131     THYRA_ASSERT_VEC_SPACES(func_name,domain,*targ_multi_vecs[k]->domain());
00132     THYRA_ASSERT_VEC_SPACES(func_name,range,*targ_multi_vecs[k]->range());
00133   }
00134 }
00135 
00136 template<class Scalar>
00137 void Thyra::apply_op_serial(
00138   const VectorSpaceBase<Scalar>   &space
00139   ,const RTOpPack::RTOpT<Scalar>  &op
00140   ,const int                      num_vecs
00141   ,const VectorBase<Scalar>*const vecs[]
00142   ,const int                      num_targ_vecs
00143   ,VectorBase<Scalar>*const       targ_vecs[]
00144   ,RTOpPack::ReductTarget         *reduct_obj
00145   ,const Index                    first_ele_offset_in
00146   ,const Index                    sub_dim_in
00147   ,const Index                    global_offset_in
00148   )
00149 {
00150    using Teuchos::Workspace;
00151   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00152   
00153   // Dimension of global sub-vector
00154   const Index
00155     full_dim       = space.dim(),
00156     global_sub_dim = sub_dim_in >= 0 ? sub_dim_in : full_dim - first_ele_offset_in;
00157   const Range1D
00158     global_sub_rng = Range1D(first_ele_offset_in,first_ele_offset_in+global_sub_dim-1);
00159 
00160   //
00161   // Get explicit views of the vector elements
00162   //
00163 
00164   Workspace<RTOpPack::ConstSubVectorView<Scalar> >         local_vecs(wss,num_vecs);
00165   Workspace<RTOpPack::SubVectorView<Scalar> >  local_targ_vecs(wss,num_targ_vecs);
00166   int k;
00167   for(k = 0; k < num_vecs; ++k) {
00168     RTOpPack::ConstSubVectorView<Scalar> &v = local_vecs[k];
00169     vecs[k]->acquireDetachedView( global_sub_rng, &v );
00170     v.setGlobalOffset( global_offset_in );
00171   }
00172   for(k = 0; k < num_targ_vecs; ++k) {
00173     RTOpPack::SubVectorView<Scalar> &v = local_targ_vecs[k];
00174     targ_vecs[k]->acquireDetachedView( global_sub_rng, &v );
00175     v.setGlobalOffset( global_offset_in );
00176   }
00177 
00178   //
00179   // Apply the reduction/transformation operator on all elements all at once!
00180   //
00181 
00182   op.apply_op(
00183     num_vecs,       num_vecs      ? &local_vecs[0]      : NULL
00184     ,num_targ_vecs, num_targ_vecs ? &local_targ_vecs[0] : NULL
00185     ,reduct_obj
00186     );
00187 
00188   //
00189   // Free (and commit) the explicit views of the vector elements which
00190   // should also inform the vectors that they have changed.
00191   //
00192 
00193   for(k = 0; k < num_vecs; ++k) {
00194     RTOpPack::ConstSubVectorView<Scalar> &v = local_vecs[k];
00195     v.setGlobalOffset(global_sub_rng.lbound());
00196     vecs[k]->releaseDetachedView(&v);
00197   }
00198   for(k = 0; k < num_targ_vecs; ++k) {
00199     RTOpPack::SubVectorView<Scalar> &v = local_targ_vecs[k];
00200     v.setGlobalOffset(global_sub_rng.lbound());
00201     targ_vecs[k]->commitDetachedView(&v);
00202   }
00203 
00204 }
00205 
00206 template<class Scalar>
00207 void Thyra::apply_op_serial(
00208   const VectorSpaceBase<Scalar>         &domain
00209   ,const VectorSpaceBase<Scalar>        &range
00210   ,const RTOpPack::RTOpT<Scalar>        &pri_op
00211   ,const int                            num_multi_vecs
00212   ,const MultiVectorBase<Scalar>*const  multi_vecs[]
00213   ,const int                            num_targ_multi_vecs
00214   ,MultiVectorBase<Scalar>*const        targ_multi_vecs[]
00215   ,RTOpPack::ReductTarget*const         reduct_objs[]
00216   ,const Index                          pri_first_ele_offset_in
00217   ,const Index                          pri_sub_dim_in
00218   ,const Index                          pri_global_offset_in
00219   ,const Index                          sec_first_ele_offset_in
00220   ,const Index                          sec_sub_dim_in
00221   )
00222 {
00223 
00224    using Teuchos::Workspace;
00225   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00226   
00227   // Primary range global sub-vector
00228   const Index
00229     range_dim          = range.dim(),
00230     pri_global_sub_dim = pri_sub_dim_in >= 0 ? pri_sub_dim_in : range_dim - pri_first_ele_offset_in;
00231   const Range1D
00232     pri_global_sub_rng = Range1D(pri_first_ele_offset_in,pri_first_ele_offset_in+pri_global_sub_dim-1);
00233   // Secondary domain
00234   const Index
00235     domain_dim         = domain.dim(),
00236     sec_global_sub_dim = sec_sub_dim_in >= 0 ? sec_sub_dim_in : domain_dim - sec_first_ele_offset_in;
00237   const Range1D
00238     sec_global_sub_rng = Range1D(sec_first_ele_offset_in,sec_first_ele_offset_in+sec_global_sub_dim-1);
00239 
00240   //
00241   // Get explicit views of the multi-vector elements
00242   //
00243 
00244   Workspace<RTOpPack::ConstSubMultiVectorView<Scalar> >         local_multi_vecs(wss,num_multi_vecs);
00245   Workspace<RTOpPack::SubMultiVectorView<Scalar> >  local_targ_multi_vecs(wss,num_targ_multi_vecs);
00246   int k;
00247   for(k = 0; k < num_multi_vecs; ++k) {
00248     RTOpPack::ConstSubMultiVectorView<Scalar> &mv = local_multi_vecs[k];
00249     multi_vecs[k]->acquireDetachedView( pri_global_sub_rng, sec_global_sub_rng, &mv );
00250     mv.setGlobalOffset( pri_global_offset_in );
00251   }
00252   for(k = 0; k < num_targ_multi_vecs; ++k) {
00253     RTOpPack::SubMultiVectorView<Scalar> &mv = local_targ_multi_vecs[k];
00254     targ_multi_vecs[k]->acquireDetachedView( pri_global_sub_rng, sec_global_sub_rng, &mv );
00255     mv.setGlobalOffset( pri_global_offset_in );
00256   }
00257 
00258   //
00259   // Apply the reduction/transformation operator one column at a time
00260   //
00261 
00262   Workspace<RTOpPack::ConstSubVectorView<Scalar> >         local_vecs(wss,num_multi_vecs);
00263   Workspace<RTOpPack::SubVectorView<Scalar> >  local_targ_vecs(wss,num_targ_multi_vecs);
00264 
00265   for(int j = 0; j < sec_global_sub_dim; ++j ) {
00266     for(k = 0; k < num_multi_vecs; ++k)       local_vecs[k]      = local_multi_vecs[k].col(j);
00267     for(k = 0; k < num_targ_multi_vecs; ++k)  local_targ_vecs[k] = local_targ_multi_vecs[k].col(j);
00268     pri_op.apply_op(
00269       num_multi_vecs,       num_multi_vecs      ? &local_vecs[0]      : NULL
00270       ,num_targ_multi_vecs, num_targ_multi_vecs ? &local_targ_vecs[0] : NULL
00271       ,reduct_objs ? reduct_objs[j] : NULL
00272       );
00273   }
00274 
00275   //
00276   // Free (and commit) the explicit views of the multi-vector elements
00277   // which should also inform the vectors that they have changed.
00278   //
00279 
00280   for(k = 0; k < num_multi_vecs; ++k) {
00281     RTOpPack::ConstSubMultiVectorView<Scalar> &mv = local_multi_vecs[k];
00282     mv.setGlobalOffset(pri_global_sub_rng.lbound());
00283     multi_vecs[k]->releaseDetachedView(&mv);
00284   }
00285   for(k = 0; k < num_targ_multi_vecs; ++k) {
00286     RTOpPack::SubMultiVectorView<Scalar> &mv = local_targ_multi_vecs[k];
00287     mv.setGlobalOffset(pri_global_sub_rng.lbound());
00288     targ_multi_vecs[k]->commitDetachedView(&mv);
00289   }
00290 
00291 }
00292 
00293 #endif // THYRA_APPLY_OP_HELPER_HPP

Generated on Thu Sep 18 12:32:31 2008 for Thyra Operator/Vector Support by doxygen 1.3.9.1