MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_VectorStdOps.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) 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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <assert.h>
00043 
00044 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00045 #include "AbstractLinAlgPack_VectorSpace.hpp"
00046 #include "AbstractLinAlgPack_VectorMutable.hpp"
00047 #include "AbstractLinAlgPack_AssertOp.hpp"
00048 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00049 #include "AbstractLinAlgPack_SpVectorView.hpp"
00050 #include "RTOp_ROp_dot_prod.h"
00051 #include "RTOp_ROp_max_abs_ele.h"
00052 #include "RTOp_ROp_sum.h"
00053 #include "RTOp_TOp_add_scalar.h"
00054 #include "RTOp_TOp_axpy.h"
00055 #include "RTOp_TOp_ele_wise_divide.h"
00056 #include "RTOp_TOp_ele_wise_prod.h"
00057 //#include "RTOp_TOp_random_vector.h"
00058 #include "RTOpPack_TOpRandomize.hpp"
00059 #include "RTOp_TOp_scale_vector.h"
00060 #include "RTOp_TOp_sign.h"
00061 #include "RTOpPack_RTOpC.hpp"
00062 #include "Teuchos_Assert.hpp"
00063 
00064 namespace {
00065 
00066 // sum
00067 static RTOpPack::RTOpC                               sum_op;
00068 static Teuchos::RCP<RTOpPack::ReductTarget>  sum_targ;
00069 // dot prod
00070 static RTOpPack::RTOpC                               dot_prod_op;
00071 static Teuchos::RCP<RTOpPack::ReductTarget>  dot_prod_targ;
00072 // number of bounded elements
00073 static RTOpPack::RTOpC                               num_bounded_op;
00074 static Teuchos::RCP<RTOpPack::ReductTarget>  num_bounded_targ;
00075 // add scalar to vector
00076 static RTOpPack::RTOpC                               add_scalar_op;
00077 // scale vector
00078 static RTOpPack::RTOpC                               scale_vector_op;
00079 // axpy
00080 static RTOpPack::RTOpC                               axpy_op;
00081 // random vector
00082 //static RTOpPack::RTOpC                               random_vector_op;
00083 static RTOpPack::TOpRandomize<AbstractLinAlgPack::value_type>  random_vector_op;
00084 // element-wise division
00085 static RTOpPack::RTOpC                               ele_wise_divide_op;
00086 // element-wise product
00087 static RTOpPack::RTOpC                               ele_wise_prod_op;
00088 
00089 // Simple class for an object that will initialize the RTOp_Server.
00090 class init_rtop_server_t {
00091 public:
00092   init_rtop_server_t() {
00093     // Operator and target obj for sum
00094     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_sum_construct(&sum_op.op()));
00095     sum_targ = sum_op.reduct_obj_create();
00096     // Operator and target obj for dot_prod
00097     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_dot_prod_construct(&dot_prod_op.op()));
00098     dot_prod_targ = dot_prod_op.reduct_obj_create();
00099     // Operator add_scalar
00100     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_add_scalar_construct(0.0,&add_scalar_op.op()));
00101     // Operator scale_vector
00102     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_scale_vector_construct(0.0,&scale_vector_op.op()));
00103     // Operator axpy
00104     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_construct(0.0,&axpy_op.op()));
00105     // Operator random_vector
00106     //TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_random_vector_construct(0.0,0.0,&random_vector_op.op()));
00107     // Operator ele_wise_divide
00108     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_divide_construct(0.0,&ele_wise_divide_op.op()));
00109     // Operator ele_wise_prod
00110     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_prod_construct(0.0,&ele_wise_prod_op.op()));
00111   }
00112 }; 
00113 
00114 // When the program starts, this object will be created and the RTOp_Server object will
00115 // be initialized before main() gets underway!
00116 init_rtop_server_t  init_rtop_server;
00117 
00118 } // end namespace
00119 
00120 AbstractLinAlgPack::value_type
00121 AbstractLinAlgPack::sum( const Vector& v_rhs )
00122 {
00123   sum_op.reduct_obj_reinit(sum_targ.ptr());
00124   const Vector* vecs[1] = { &v_rhs };
00125   apply_op(sum_op,1,vecs,0,NULL,&*sum_targ);
00126   return RTOp_ROp_sum_val(sum_op(*sum_targ));
00127 }
00128 
00129 AbstractLinAlgPack::value_type
00130 AbstractLinAlgPack::dot( const Vector& v_rhs1, const Vector& v_rhs2 )
00131 {
00132   dot_prod_op.reduct_obj_reinit(dot_prod_targ.ptr());
00133   const Vector* vecs[2] = { &v_rhs1, &v_rhs2 };
00134   apply_op(dot_prod_op,2,vecs,0,NULL,&*dot_prod_targ);
00135   return RTOp_ROp_dot_prod_val(dot_prod_op(*dot_prod_targ));
00136 }
00137 
00138 AbstractLinAlgPack::value_type
00139 AbstractLinAlgPack::dot( const Vector& v_rhs1, const SpVectorSlice& sv_rhs2 )
00140 {
00141   VopV_assert_compatibility(v_rhs1,sv_rhs2 );
00142   if( sv_rhs2.nz() ) {
00143     VectorSpace::vec_mut_ptr_t
00144       v_rhs2 = v_rhs1.space().create_member();
00145     v_rhs2->set_sub_vector(sub_vec_view(sv_rhs2));
00146     return dot(v_rhs1,*v_rhs2);
00147   }
00148   return 0.0;
00149 }
00150 
00151 void AbstractLinAlgPack::max_abs_ele(
00152   const Vector& v, value_type* max_v_j, index_type* max_j
00153   )
00154 {
00155   TEUCHOS_TEST_FOR_EXCEPT( !(  max_v_j && max_j  ) );
00156   RTOpPack::RTOpC op;
00157   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_max_abs_ele_construct(&op.op()));
00158   Teuchos::RCP<RTOpPack::ReductTarget> reduct_obj = op.reduct_obj_create();
00159   const Vector* vecs[1] = { &v };
00160   apply_op(op,1,vecs,0,NULL,&*reduct_obj);
00161   RTOp_value_index_type val = RTOp_ROp_max_abs_ele_val(op(*reduct_obj));
00162   *max_v_j = val.value;
00163   *max_j   = val.index;
00164 }
00165 
00166 void AbstractLinAlgPack::Vp_S( VectorMutable* v_lhs, const value_type& alpha )
00167 {
00168 #ifdef TEUCHOS_DEBUG
00169   TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vp_S(...), Error!");
00170 #endif
00171   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_add_scalar_set_alpha(alpha,&add_scalar_op.op()));
00172   VectorMutable* targ_vecs[1] = { v_lhs };
00173   apply_op(add_scalar_op,0,NULL,1,targ_vecs,NULL);
00174 }
00175 
00176 void AbstractLinAlgPack::Vt_S( VectorMutable* v_lhs, const value_type& alpha )
00177 {
00178 #ifdef TEUCHOS_DEBUG
00179   TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vt_S(...), Error!");
00180 #endif
00181   if( alpha == 0.0 ) {
00182     *v_lhs = 0.0;
00183   }
00184   else if( alpha != 1.0 ) {
00185     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_scale_vector_set_alpha( alpha, &scale_vector_op.op() ));
00186     VectorMutable* targ_vecs[1] = { v_lhs };
00187     apply_op(scale_vector_op,0,NULL,1,targ_vecs,NULL);
00188   }
00189 }
00190 
00191 void AbstractLinAlgPack::Vp_StV(
00192   VectorMutable* v_lhs, const value_type& alpha, const Vector& v_rhs)
00193 {
00194 #ifdef TEUCHOS_DEBUG
00195   TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vp_StV(...), Error!");
00196 #endif
00197   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_set_alpha( alpha, &axpy_op.op() ));
00198   const Vector*  vecs[1]      = { &v_rhs };
00199   VectorMutable* targ_vecs[1] = { v_lhs  };
00200   apply_op(axpy_op,1,vecs,1,targ_vecs,NULL);
00201 }
00202 
00203 void AbstractLinAlgPack::Vp_StV(
00204   VectorMutable* y, const value_type& a, const SpVectorSlice& sx )
00205 {
00206   Vp_V_assert_compatibility(y,sx);
00207   if( sx.nz() ) {
00208     VectorSpace::vec_mut_ptr_t
00209         x = y->space().create_member();
00210     x->set_sub_vector(sub_vec_view(sx));
00211     Vp_StV( y, a, *x );
00212   }
00213 }
00214 
00215 void AbstractLinAlgPack::ele_wise_prod(
00216   const value_type& alpha, const Vector& v_rhs1, const Vector& v_rhs2
00217   , VectorMutable* v_lhs )
00218 {
00219 #ifdef TEUCHOS_DEBUG
00220   TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"ele_wise_prod(...), Error");
00221 #endif
00222   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_prod_set_alpha(alpha,&ele_wise_prod_op.op()));
00223   const Vector*   vecs[2]      = { &v_rhs1, &v_rhs2 };
00224   VectorMutable*  targ_vecs[1] = { v_lhs };
00225   apply_op(ele_wise_prod_op,2,vecs,1,targ_vecs,NULL);
00226 }
00227 
00228 void AbstractLinAlgPack::ele_wise_divide(
00229   const value_type& alpha, const Vector& v_rhs1, const Vector& v_rhs2
00230   , VectorMutable* v_lhs )
00231 {
00232 #ifdef TEUCHOS_DEBUG
00233   TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"ele_wise_divide(...), Error");
00234 #endif
00235   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_divide_set_alpha(alpha,&ele_wise_divide_op.op()));
00236   const int num_vecs = 2;
00237   const Vector*   vecs[2]      = { &v_rhs1, &v_rhs2 };
00238   VectorMutable*  targ_vecs[1] = { v_lhs };
00239   apply_op(ele_wise_divide_op,2,vecs,1,targ_vecs,NULL);
00240 }
00241 
00242 void AbstractLinAlgPack::seed_random_vector_generator( unsigned int s )
00243 {
00244   random_vector_op.set_seed(s);
00245 }
00246 
00247 void AbstractLinAlgPack::random_vector( value_type l, value_type u, VectorMutable* v )
00248 {
00249 #ifdef TEUCHOS_DEBUG
00250   TEUCHOS_TEST_FOR_EXCEPTION(v==NULL,std::logic_error,"Vt_S(...), Error");
00251 #endif
00252   //TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_random_vector_set_bounds( l, u, &random_vector_op.op() ));
00253   random_vector_op.set_bounds(l,u);
00254   VectorMutable* targ_vecs[1] = { v };
00255   apply_op(random_vector_op,0,NULL,1,targ_vecs,NULL);
00256   
00257 }
00258 
00259 void AbstractLinAlgPack::sign(
00260   const Vector      &v
00261   ,VectorMutable    *z
00262   )
00263 {
00264   RTOpPack::RTOpC op;
00265   TEUCHOS_TEST_FOR_EXCEPT( !( 0==RTOp_TOp_sign_construct(&op.op()) ) );
00266   const Vector*   vecs[1]      = { &v };
00267   VectorMutable*  targ_vecs[1] = { z  };
00268   apply_op(op,1,vecs,1,targ_vecs,NULL);
00269 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines