AbstractLinAlgPack_apply_op_helper.cpp

00001 // /////////////////////////////////////////////////////////////////////////////
00002 // apply_op_helper.cpp
00003 
00004 #include <assert.h>
00005 
00006 #include "AbstractLinAlgPack_apply_op_helper.hpp"
00007 #include "AbstractLinAlgPack_VectorSpace.hpp"
00008 #include "AbstractLinAlgPack_VectorMutable.hpp"
00009 #include "Teuchos_Workspace.hpp"
00010 #include "Teuchos_TestForException.hpp"
00011 
00012 void AbstractLinAlgPack::apply_op_validate_input(
00013   const char func_name[]
00014   ,const RTOpPack::RTOp& op
00015   ,const size_t num_vecs, const Vector* vecs[]
00016   ,const size_t num_targ_vecs, VectorMutable* targ_vecs[]
00017   ,RTOpPack::ReductTarget *reduct_obj
00018   ,const index_type first_ele_in, const index_type sub_dim_in, const index_type global_offset_in
00019   )
00020 {
00021   const VectorSpace
00022     &space = (num_vecs ? vecs[0]->space() : targ_vecs[0]->space() );
00023   const index_type
00024     dim = space.dim();
00025   TEST_FOR_EXCEPTION(
00026     global_offset_in < 0, std::logic_error
00027     ,func_name << " : Error!  global_offset_in = "
00028     <<global_offset_in<<" is not valid" );
00029   TEST_FOR_EXCEPTION(
00030     first_ele_in > dim, std::logic_error
00031     ,func_name << " : Error!  first_ele_in = "
00032     <<first_ele_in<<" is not compatible with space.dim() = " << dim );
00033   TEST_FOR_EXCEPTION(
00034     sub_dim_in < 0 || (sub_dim_in > 0 && sub_dim_in > dim-(first_ele_in-1)), std::logic_error
00035     ,func_name << " : Error!  first_ele_in = "
00036     <<first_ele_in<<" and sub_dim_in = "<<sub_dim_in
00037     <<" is not compatible with space.dim() = " << dim );
00038   {for(int k = 0; k < num_vecs; ++k) {
00039     const bool is_compatible = space.is_compatible(vecs[k]->space());
00040     TEST_FOR_EXCEPTION(
00041       !is_compatible, VectorSpace::IncompatibleVectorSpaces
00042       ,func_name << " : Error!  vecs["<<k<<"]->space() of type \'"
00043       << typeName(vecs[k]->space()) << "\' "
00044       << " with dimension vecs["<<k<<"].dim() = " << vecs[k]->dim()
00045       << " is not compatible with space of type \'"
00046       << typeName(space) << " with dimmension space.dim() = " << dim );
00047   }}
00048   {for(int k = 0; k < num_targ_vecs; ++k) {
00049     const bool is_compatible = space.is_compatible(targ_vecs[k]->space());
00050     TEST_FOR_EXCEPTION(
00051       !is_compatible, VectorSpace::IncompatibleVectorSpaces
00052       ,func_name << " : Error!  targ_vecs["<<k<<"]->space() of type \'"
00053       << typeName(targ_vecs[k]->space()) << "\' "
00054       << " with dimension targ_vecs["<<k<<"].dim() = " << targ_vecs[k]->dim()
00055       << " is not compatible with space of type \'"
00056       << typeName(space) << " with dimmension space.dim() = " << dim );
00057   }}
00058 }
00059 
00060 void AbstractLinAlgPack::apply_op_serial(
00061   const RTOpPack::RTOp& op
00062   ,const size_t num_vecs, const Vector* vecs[]
00063   ,const size_t num_targ_vecs, VectorMutable* targ_vecs[]
00064   ,RTOpPack::ReductTarget *reduct_obj
00065   ,const index_type first_ele_in, const index_type sub_dim_in, const index_type global_offset_in
00066   )
00067 {
00068    using Teuchos::Workspace;
00069   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00070   
00071   // Dimension of global sub-vector
00072   const VectorSpace
00073     &space = ( num_vecs ? vecs[0]->space() : targ_vecs[0]->space() );
00074   const index_type
00075     full_dim       = space.dim(),
00076     global_sub_dim = sub_dim_in ? sub_dim_in     : full_dim - (first_ele_in-1);
00077   const Range1D
00078     global_sub_rng = Range1D(first_ele_in,(first_ele_in-1)+global_sub_dim);
00079 
00080   //
00081   // Get explicit views of the vector elements
00082   //
00083 
00084   Workspace<RTOpPack::ConstSubVectorView<value_type> >         local_vecs(wss,num_vecs);
00085   Workspace<RTOpPack::SubVectorView<value_type> >  local_targ_vecs(wss,num_targ_vecs);
00086   int k;
00087   for(k = 0; k < num_vecs; ++k) {
00088     RTOpPack::SubVector _v;
00089     vecs[k]->get_sub_vector( global_sub_rng, &_v );
00090     (local_vecs[k] = _v).setGlobalOffset( _v.globalOffset() + global_offset_in );
00091   }
00092   for(k = 0; k < num_targ_vecs; ++k) {
00093     RTOpPack::MutableSubVector _v;
00094     targ_vecs[k]->get_sub_vector( global_sub_rng, &_v );
00095     (local_targ_vecs[k] = _v).setGlobalOffset( _v.globalOffset() + global_offset_in );
00096   }
00097 
00098   //
00099   // Apply the reduction/transformation operator on all elements all at once!
00100   //
00101 
00102   op.apply_op(
00103     num_vecs,       num_vecs      ? &local_vecs[0]      : NULL
00104     ,num_targ_vecs, num_targ_vecs ? &local_targ_vecs[0] : NULL
00105     ,reduct_obj
00106     );
00107 
00108   //
00109   // Free (and commit) the explicit views of the vector elements
00110   // which should also inform the vectors that they have
00111   // changed.
00112   //
00113 
00114   for(k = 0; k < num_vecs; ++k) {
00115     RTOpPack::ConstSubVectorView<value_type> &v = local_vecs[k];
00116     v.setGlobalOffset( v.globalOffset() - global_offset_in );
00117     RTOpPack::SubVector _v = v;
00118     vecs[k]->free_sub_vector(&_v);
00119   }
00120   for(k = 0; k < num_targ_vecs; ++k) {
00121     RTOpPack::SubVectorView<value_type> &v = local_targ_vecs[k];
00122     v.setGlobalOffset( v.globalOffset() - global_offset_in );
00123     RTOpPack::MutableSubVector _v = v;
00124     targ_vecs[k]->commit_sub_vector(&_v);
00125   }
00126 
00127 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:09:14 2011 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.6.3