RTOp Package Browser (Single Doxygen Collection) Version of the Day
RTOpPack_TOpSetSubVector_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // RTOp: Interfaces and Support Software for Vector Reduction Transformation
00005 //       Operations
00006 //                Copyright (2006) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
00031 #define RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
00032 
00033 
00034 namespace RTOpPack {
00035 
00036 
00037 template<class Scalar>
00038 TOpSetSubVector<Scalar>::TOpSetSubVector()
00039   :RTOpT<Scalar>("TOpSetSubVector")
00040 {}
00041 
00042 
00043 template<class Scalar>
00044 TOpSetSubVector<Scalar>::TOpSetSubVector( const SparseSubVectorT<Scalar> &sub_vec )
00045   :RTOpT<Scalar>("TOpSetSubVector")
00046 {
00047   set_sub_vec(sub_vec);
00048 }
00049 
00050 
00051 template<class Scalar>
00052 void TOpSetSubVector<Scalar>::set_sub_vec( const SparseSubVectorT<Scalar> &sub_vec )
00053 {
00054   sub_vec_ = sub_vec;
00055 }
00056 
00057 
00058 // Overridden from RTOpT
00059 
00060 
00061 template<class Scalar>
00062 bool TOpSetSubVector<Scalar>::coord_invariant_impl() const
00063 {
00064   return false;
00065 }
00066 
00067 
00068 template<class Scalar>
00069 void TOpSetSubVector<Scalar>::apply_op_impl(
00070   const ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
00071   const ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
00072   const Ptr<ReductTarget> &reduct_obj
00073   ) const
00074 {
00075 
00076   typedef typename Teuchos::ArrayRCP<Scalar>::iterator iter_t;
00077   typedef typename Teuchos::ArrayRCP<const Scalar>::iterator const_iter_t;
00078   typedef typename Teuchos::ArrayRCP<const Teuchos_Index>::iterator const_indices_iter_t;
00079   
00080   validate_apply_op( *this, 0, 1, false, sub_vecs, targ_sub_vecs, reduct_obj.getConst() );
00081 
00082   // Get the dense subvector data that we will set
00083   const SubVectorView<Scalar> &z = targ_sub_vecs[0];
00084   const index_type z_global_offset = z.globalOffset();
00085   const index_type z_sub_dim = z.subDim();
00086   iter_t z_val = z.values().begin();
00087   const ptrdiff_t z_val_s = z.stride();
00088 
00089   // Get the sparse subvector data
00090   const index_type v_global_offset = sub_vec_.globalOffset();
00091   const index_type v_sub_dim = sub_vec_.subDim();
00092   const index_type v_sub_nz = sub_vec_.subNz();
00093   const_iter_t v_val = sub_vec_.values().begin();
00094   const ptrdiff_t v_val_s = sub_vec_.valuesStride();
00095   const bool has_v_ind = !is_null(sub_vec_.indices());
00096   const_indices_iter_t v_ind = sub_vec_.indices().begin();
00097   const ptrdiff_t v_ind_s = sub_vec_.indicesStride();
00098   const ptrdiff_t v_l_off = sub_vec_.localOffset();
00099   //const bool v_sorted = sub_vec_.isSorted();
00100  
00101   //
00102   // Set part of the sub-vector v for this chunk z.
00103   //
00104 
00105   if( v_global_offset + v_sub_dim < z_global_offset + 1
00106     || z_global_offset + z_sub_dim < v_global_offset + 1 )
00107   {
00108      // The sub-vector that we are setting does not overlap with this vector
00109      // chunk!
00110     return;
00111   }
00112 
00113   if( v_sub_nz == 0 )
00114     return; // The sparse sub-vector we are reading from is empty?
00115 
00116   // Get the number of elements that overlap
00117   index_type num_overlap;
00118   if( v_global_offset <= z_global_offset ) {
00119     if( v_global_offset + v_sub_dim >= z_global_offset + z_sub_dim )
00120       num_overlap = z_sub_dim;
00121     else
00122       num_overlap = (v_global_offset + v_sub_dim) - z_global_offset;
00123   }
00124   else {
00125     if( z_global_offset + z_sub_dim >= v_global_offset + v_sub_dim )
00126       num_overlap = v_sub_dim;
00127     else
00128       num_overlap = (z_global_offset + z_sub_dim) - v_global_offset;
00129   }
00130  
00131   // Set the part of the sub-vector that overlaps
00132   if (has_v_ind) {
00133     // Sparse elements
00134     // Set the overlapping elements to zero first.
00135     if( v_global_offset >= z_global_offset )
00136       z_val += (v_global_offset - z_global_offset) * z_val_s;
00137     for( index_type k = 0; k < num_overlap; ++k, z_val += z_val_s )
00138       *z_val = 0.0;
00139     // Now set the sparse entries
00140     z_val = targ_sub_vecs[0].values().begin();
00141     for( index_type k = 0; k < v_sub_nz; ++k, v_val += v_val_s, v_ind += v_ind_s ) {
00142       const index_type i = v_global_offset + v_l_off + (*v_ind);
00143       if( z_global_offset < i && i <= z_global_offset + z_sub_dim )
00144         z_val[ z_val_s * (i - z_global_offset - 1) ] = *v_val;
00145     }
00146     // ToDo: Implement a faster version for v sorted and eliminate the
00147     // if statement in the above loop.
00148   }
00149   else {
00150     // Dense elements
00151     if( v_global_offset <= z_global_offset )
00152       v_val += (z_global_offset - v_global_offset) * v_val_s;
00153     else
00154       z_val += (v_global_offset - z_global_offset) * z_val_s;
00155     for( index_type k = 0; k < num_overlap; ++k, v_val += v_val_s, z_val += z_val_s )
00156       *z_val = *v_val;
00157   }
00158 }
00159 
00160 
00161 } // namespace RTOpPack
00162 
00163 
00164 #endif // RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines