AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_apply_op_helper.cpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00006 //                  Copyright (2003) 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 
00044 // /////////////////////////////////////////////////////////////////////////////
00045 // apply_op_helper.cpp
00046 
00047 #include <assert.h>
00048 
00049 #include "AbstractLinAlgPack_apply_op_helper.hpp"
00050 #include "AbstractLinAlgPack_VectorSpace.hpp"
00051 #include "AbstractLinAlgPack_VectorMutable.hpp"
00052 #include "Teuchos_Workspace.hpp"
00053 #include "Teuchos_Assert.hpp"
00054 
00055 void AbstractLinAlgPack::apply_op_validate_input(
00056   const char func_name[]
00057   ,const RTOpPack::RTOp& op
00058   ,const size_t num_vecs, const Vector* vecs[]
00059   ,const size_t num_targ_vecs, VectorMutable* targ_vecs[]
00060   ,RTOpPack::ReductTarget *reduct_obj
00061   ,const index_type first_ele_in, const index_type sub_dim_in, const index_type global_offset_in
00062   )
00063 {
00064   const VectorSpace
00065     &space = (num_vecs ? vecs[0]->space() : targ_vecs[0]->space() );
00066   const index_type
00067     dim = space.dim();
00068   TEUCHOS_TEST_FOR_EXCEPTION(
00069     global_offset_in < 0, std::logic_error
00070     ,func_name << " : Error!  global_offset_in = "
00071     <<global_offset_in<<" is not valid" );
00072   TEUCHOS_TEST_FOR_EXCEPTION(
00073     first_ele_in > dim, std::logic_error
00074     ,func_name << " : Error!  first_ele_in = "
00075     <<first_ele_in<<" is not compatible with space.dim() = " << dim );
00076   TEUCHOS_TEST_FOR_EXCEPTION(
00077     sub_dim_in < 0 || (sub_dim_in > 0 && sub_dim_in > dim-(first_ele_in-1)), std::logic_error
00078     ,func_name << " : Error!  first_ele_in = "
00079     <<first_ele_in<<" and sub_dim_in = "<<sub_dim_in
00080     <<" is not compatible with space.dim() = " << dim );
00081   {for(int k = 0; k < num_vecs; ++k) {
00082     const bool is_compatible = space.is_compatible(vecs[k]->space());
00083     TEUCHOS_TEST_FOR_EXCEPTION(
00084       !is_compatible, VectorSpace::IncompatibleVectorSpaces
00085       ,func_name << " : Error!  vecs["<<k<<"]->space() of type \'"
00086       << typeName(vecs[k]->space()) << "\' "
00087       << " with dimension vecs["<<k<<"].dim() = " << vecs[k]->dim()
00088       << " is not compatible with space of type \'"
00089       << typeName(space) << " with dimmension space.dim() = " << dim );
00090   }}
00091   {for(int k = 0; k < num_targ_vecs; ++k) {
00092     const bool is_compatible = space.is_compatible(targ_vecs[k]->space());
00093     TEUCHOS_TEST_FOR_EXCEPTION(
00094       !is_compatible, VectorSpace::IncompatibleVectorSpaces
00095       ,func_name << " : Error!  targ_vecs["<<k<<"]->space() of type \'"
00096       << typeName(targ_vecs[k]->space()) << "\' "
00097       << " with dimension targ_vecs["<<k<<"].dim() = " << targ_vecs[k]->dim()
00098       << " is not compatible with space of type \'"
00099       << typeName(space) << " with dimmension space.dim() = " << dim );
00100   }}
00101 }
00102 
00103 void AbstractLinAlgPack::apply_op_serial(
00104   const RTOpPack::RTOp& op
00105   ,const size_t num_vecs, const Vector* vecs[]
00106   ,const size_t num_targ_vecs, VectorMutable* targ_vecs[]
00107   ,RTOpPack::ReductTarget *reduct_obj
00108   ,const index_type first_ele_in, const index_type sub_dim_in, const index_type global_offset_in
00109   )
00110 {
00111    using Teuchos::Workspace;
00112   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00113   
00114   // Dimension of global sub-vector
00115   const VectorSpace
00116     &space = ( num_vecs ? vecs[0]->space() : targ_vecs[0]->space() );
00117   const index_type
00118     full_dim       = space.dim(),
00119     global_sub_dim = sub_dim_in ? sub_dim_in     : full_dim - (first_ele_in-1);
00120   const Range1D
00121     global_sub_rng = Range1D(first_ele_in,(first_ele_in-1)+global_sub_dim);
00122 
00123   //
00124   // Get explicit views of the vector elements
00125   //
00126 
00127   Workspace<RTOpPack::ConstSubVectorView<value_type> >         local_vecs(wss,num_vecs);
00128   Workspace<RTOpPack::SubVectorView<value_type> >  local_targ_vecs(wss,num_targ_vecs);
00129   int k;
00130   for(k = 0; k < num_vecs; ++k) {
00131     RTOpPack::SubVector _v;
00132     vecs[k]->get_sub_vector( global_sub_rng, &_v );
00133     (local_vecs[k] = _v).setGlobalOffset( _v.globalOffset() + global_offset_in );
00134   }
00135   for(k = 0; k < num_targ_vecs; ++k) {
00136     RTOpPack::MutableSubVector _v;
00137     targ_vecs[k]->get_sub_vector( global_sub_rng, &_v );
00138     (local_targ_vecs[k] = _v).setGlobalOffset( _v.globalOffset() + global_offset_in );
00139   }
00140 
00141   //
00142   // Apply the reduction/transformation operator on all elements all at once!
00143   //
00144 
00145   op.apply_op(
00146     num_vecs,       num_vecs      ? &local_vecs[0]      : NULL
00147     ,num_targ_vecs, num_targ_vecs ? &local_targ_vecs[0] : NULL
00148     ,reduct_obj
00149     );
00150 
00151   //
00152   // Free (and commit) the explicit views of the vector elements
00153   // which should also inform the vectors that they have
00154   // changed.
00155   //
00156 
00157   for(k = 0; k < num_vecs; ++k) {
00158     RTOpPack::ConstSubVectorView<value_type> &v = local_vecs[k];
00159     v.setGlobalOffset( v.globalOffset() - global_offset_in );
00160     RTOpPack::SubVector _v = v;
00161     vecs[k]->free_sub_vector(&_v);
00162   }
00163   for(k = 0; k < num_targ_vecs; ++k) {
00164     RTOpPack::SubVectorView<value_type> &v = local_targ_vecs[k];
00165     v.setGlobalOffset( v.globalOffset() - global_offset_in );
00166     RTOpPack::MutableSubVector _v = v;
00167     targ_vecs[k]->commit_sub_vector(&_v);
00168   }
00169 
00170 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends