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 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 
00043 #ifndef RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
00044 #define RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
00045 
00046 
00047 namespace RTOpPack {
00048 
00049 
00050 template<class Scalar>
00051 TOpSetSubVector<Scalar>::TOpSetSubVector()
00052   :RTOpT<Scalar>("TOpSetSubVector")
00053 {}
00054 
00055 
00056 template<class Scalar>
00057 TOpSetSubVector<Scalar>::TOpSetSubVector( const SparseSubVectorT<Scalar> &sub_vec )
00058   :RTOpT<Scalar>("TOpSetSubVector")
00059 {
00060   set_sub_vec(sub_vec);
00061 }
00062 
00063 
00064 template<class Scalar>
00065 void TOpSetSubVector<Scalar>::set_sub_vec( const SparseSubVectorT<Scalar> &sub_vec )
00066 {
00067   sub_vec_ = sub_vec;
00068 }
00069 
00070 
00071 // Overridden from RTOpT
00072 
00073 
00074 template<class Scalar>
00075 bool TOpSetSubVector<Scalar>::coord_invariant_impl() const
00076 {
00077   return false;
00078 }
00079 
00080 
00081 template<class Scalar>
00082 void TOpSetSubVector<Scalar>::apply_op_impl(
00083   const ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
00084   const ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
00085   const Ptr<ReductTarget> &reduct_obj
00086   ) const
00087 {
00088 
00089   typedef typename Teuchos::ArrayRCP<Scalar>::iterator iter_t;
00090   typedef typename Teuchos::ArrayRCP<const Scalar>::iterator const_iter_t;
00091   typedef typename Teuchos::ArrayRCP<const Teuchos_Ordinal>::iterator const_indices_iter_t;
00092   
00093   validate_apply_op( *this, 0, 1, false, sub_vecs, targ_sub_vecs, reduct_obj.getConst() );
00094 
00095   // Get the dense subvector data that we will set
00096   const SubVectorView<Scalar> &z = targ_sub_vecs[0];
00097   const index_type z_global_offset = z.globalOffset();
00098   const index_type z_sub_dim = z.subDim();
00099   iter_t z_val = z.values().begin();
00100   const ptrdiff_t z_val_s = z.stride();
00101 
00102   // Get the sparse subvector data
00103   const index_type v_global_offset = sub_vec_.globalOffset();
00104   const index_type v_sub_dim = sub_vec_.subDim();
00105   const index_type v_sub_nz = sub_vec_.subNz();
00106   const_iter_t v_val = sub_vec_.values().begin();
00107   const ptrdiff_t v_val_s = sub_vec_.valuesStride();
00108   const bool has_v_ind = !is_null(sub_vec_.indices());
00109   const_indices_iter_t v_ind = sub_vec_.indices().begin();
00110   const ptrdiff_t v_ind_s = sub_vec_.indicesStride();
00111   const ptrdiff_t v_l_off = sub_vec_.localOffset();
00112   //const bool v_sorted = sub_vec_.isSorted();
00113  
00114   //
00115   // Set part of the sub-vector v for this chunk z.
00116   //
00117 
00118   if( v_global_offset + v_sub_dim < z_global_offset + 1
00119     || z_global_offset + z_sub_dim < v_global_offset + 1 )
00120   {
00121      // The sub-vector that we are setting does not overlap with this vector
00122      // chunk!
00123     return;
00124   }
00125 
00126   if( v_sub_nz == 0 )
00127     return; // The sparse sub-vector we are reading from is empty?
00128 
00129   // Get the number of elements that overlap
00130   index_type num_overlap;
00131   if( v_global_offset <= z_global_offset ) {
00132     if( v_global_offset + v_sub_dim >= z_global_offset + z_sub_dim )
00133       num_overlap = z_sub_dim;
00134     else
00135       num_overlap = (v_global_offset + v_sub_dim) - z_global_offset;
00136   }
00137   else {
00138     if( z_global_offset + z_sub_dim >= v_global_offset + v_sub_dim )
00139       num_overlap = v_sub_dim;
00140     else
00141       num_overlap = (z_global_offset + z_sub_dim) - v_global_offset;
00142   }
00143  
00144   // Set the part of the sub-vector that overlaps
00145   if (has_v_ind) {
00146     // Sparse elements
00147     // Set the overlapping elements to zero first.
00148     if( v_global_offset >= z_global_offset )
00149       z_val += (v_global_offset - z_global_offset) * z_val_s;
00150     for( index_type k = 0; k < num_overlap; ++k, z_val += z_val_s )
00151       *z_val = 0.0;
00152     // Now set the sparse entries
00153     z_val = targ_sub_vecs[0].values().begin();
00154     for( index_type k = 0; k < v_sub_nz; ++k, v_val += v_val_s, v_ind += v_ind_s ) {
00155       const index_type i = v_global_offset + v_l_off + (*v_ind);
00156       if( z_global_offset < i && i <= z_global_offset + z_sub_dim )
00157         z_val[ z_val_s * (i - z_global_offset - 1) ] = *v_val;
00158     }
00159     // ToDo: Implement a faster version for v sorted and eliminate the
00160     // if statement in the above loop.
00161   }
00162   else {
00163     // Dense elements
00164     if( v_global_offset <= z_global_offset )
00165       v_val += (z_global_offset - v_global_offset) * v_val_s;
00166     else
00167       z_val += (v_global_offset - z_global_offset) * z_val_s;
00168     for( index_type k = 0; k < num_overlap; ++k, v_val += v_val_s, z_val += z_val_s )
00169       *z_val = *v_val;
00170   }
00171 }
00172 
00173 
00174 } // namespace RTOpPack
00175 
00176 
00177 #endif // RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines