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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <assert.h>
00030 
00031 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00032 #include "AbstractLinAlgPack_VectorSpace.hpp"
00033 #include "AbstractLinAlgPack_VectorMutable.hpp"
00034 #include "AbstractLinAlgPack_AssertOp.hpp"
00035 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00036 #include "AbstractLinAlgPack_SpVectorView.hpp"
00037 #include "RTOp_ROp_dot_prod.h"
00038 #include "RTOp_ROp_max_abs_ele.h"
00039 #include "RTOp_ROp_sum.h"
00040 #include "RTOp_TOp_add_scalar.h"
00041 #include "RTOp_TOp_axpy.h"
00042 #include "RTOp_TOp_ele_wise_divide.h"
00043 #include "RTOp_TOp_ele_wise_prod.h"
00044 //#include "RTOp_TOp_random_vector.h"
00045 #include "RTOpPack_TOpRandomize.hpp"
00046 #include "RTOp_TOp_scale_vector.h"
00047 #include "RTOp_TOp_sign.h"
00048 #include "RTOpPack_RTOpC.hpp"
00049 #include "Teuchos_TestForException.hpp"
00050 
00051 namespace {
00052 
00053 // sum
00054 static RTOpPack::RTOpC                               sum_op;
00055 static Teuchos::RCP<RTOpPack::ReductTarget>  sum_targ;
00056 // dot prod
00057 static RTOpPack::RTOpC                               dot_prod_op;
00058 static Teuchos::RCP<RTOpPack::ReductTarget>  dot_prod_targ;
00059 // number of bounded elements
00060 static RTOpPack::RTOpC                               num_bounded_op;
00061 static Teuchos::RCP<RTOpPack::ReductTarget>  num_bounded_targ;
00062 // add scalar to vector
00063 static RTOpPack::RTOpC                               add_scalar_op;
00064 // scale vector
00065 static RTOpPack::RTOpC                               scale_vector_op;
00066 // axpy
00067 static RTOpPack::RTOpC                               axpy_op;
00068 // random vector
00069 //static RTOpPack::RTOpC                               random_vector_op;
00070 static RTOpPack::TOpRandomize<AbstractLinAlgPack::value_type>  random_vector_op;
00071 // element-wise division
00072 static RTOpPack::RTOpC                               ele_wise_divide_op;
00073 // element-wise product
00074 static RTOpPack::RTOpC                               ele_wise_prod_op;
00075 
00076 // Simple class for an object that will initialize the RTOp_Server.
00077 class init_rtop_server_t {
00078 public:
00079   init_rtop_server_t() {
00080     // Operator and target obj for sum
00081     TEST_FOR_EXCEPT(0!=RTOp_ROp_sum_construct(&sum_op.op()));
00082     sum_targ = sum_op.reduct_obj_create();
00083     // Operator and target obj for dot_prod
00084     TEST_FOR_EXCEPT(0!=RTOp_ROp_dot_prod_construct(&dot_prod_op.op()));
00085     dot_prod_targ = dot_prod_op.reduct_obj_create();
00086     // Operator add_scalar
00087     TEST_FOR_EXCEPT(0!=RTOp_TOp_add_scalar_construct(0.0,&add_scalar_op.op()));
00088     // Operator scale_vector
00089     TEST_FOR_EXCEPT(0!=RTOp_TOp_scale_vector_construct(0.0,&scale_vector_op.op()));
00090     // Operator axpy
00091     TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_construct(0.0,&axpy_op.op()));
00092     // Operator random_vector
00093     //TEST_FOR_EXCEPT(0!=RTOp_TOp_random_vector_construct(0.0,0.0,&random_vector_op.op()));
00094     // Operator ele_wise_divide
00095     TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_divide_construct(0.0,&ele_wise_divide_op.op()));
00096     // Operator ele_wise_prod
00097     TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_prod_construct(0.0,&ele_wise_prod_op.op()));
00098   }
00099 }; 
00100 
00101 // When the program starts, this object will be created and the RTOp_Server object will
00102 // be initialized before main() gets underway!
00103 init_rtop_server_t  init_rtop_server;
00104 
00105 } // end namespace
00106 
00107 AbstractLinAlgPack::value_type
00108 AbstractLinAlgPack::sum( const Vector& v_rhs )
00109 {
00110   sum_op.reduct_obj_reinit(&*sum_targ);
00111   const Vector* vecs[1] = { &v_rhs };
00112   apply_op(sum_op,1,vecs,0,NULL,&*sum_targ);
00113   return RTOp_ROp_sum_val(sum_op(*sum_targ));
00114 }
00115 
00116 AbstractLinAlgPack::value_type
00117 AbstractLinAlgPack::dot( const Vector& v_rhs1, const Vector& v_rhs2 )
00118 {
00119   dot_prod_op.reduct_obj_reinit(&*dot_prod_targ);
00120   const Vector* vecs[2] = { &v_rhs1, &v_rhs2 };
00121   apply_op(dot_prod_op,2,vecs,0,NULL,&*dot_prod_targ);
00122   return RTOp_ROp_dot_prod_val(dot_prod_op(*dot_prod_targ));
00123 }
00124 
00125 AbstractLinAlgPack::value_type
00126 AbstractLinAlgPack::dot( const Vector& v_rhs1, const SpVectorSlice& sv_rhs2 )
00127 {
00128   VopV_assert_compatibility(v_rhs1,sv_rhs2 );
00129   if( sv_rhs2.nz() ) {
00130     VectorSpace::vec_mut_ptr_t
00131       v_rhs2 = v_rhs1.space().create_member();
00132     v_rhs2->set_sub_vector(sub_vec_view(sv_rhs2));
00133     return dot(v_rhs1,*v_rhs2);
00134   }
00135   return 0.0;
00136 }
00137 
00138 void AbstractLinAlgPack::max_abs_ele(
00139   const Vector& v, value_type* max_v_j, index_type* max_j
00140   )
00141 {
00142   TEST_FOR_EXCEPT( !(  max_v_j && max_j  ) );
00143   RTOpPack::RTOpC op;
00144   TEST_FOR_EXCEPT(0!=RTOp_ROp_max_abs_ele_construct(&op.op()));
00145   Teuchos::RCP<RTOpPack::ReductTarget> reduct_obj = op.reduct_obj_create();
00146   const Vector* vecs[1] = { &v };
00147   apply_op(op,1,vecs,0,NULL,&*reduct_obj);
00148   RTOp_value_index_type val = RTOp_ROp_max_abs_ele_val(op(*reduct_obj));
00149   *max_v_j = val.value;
00150   *max_j   = val.index;
00151 }
00152 
00153 void AbstractLinAlgPack::Vp_S( VectorMutable* v_lhs, const value_type& alpha )
00154 {
00155 #ifdef TEUCHOS_DEBUG
00156   TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vp_S(...), Error!");
00157 #endif
00158   TEST_FOR_EXCEPT(0!=RTOp_TOp_add_scalar_set_alpha(alpha,&add_scalar_op.op()));
00159   VectorMutable* targ_vecs[1] = { v_lhs };
00160   apply_op(add_scalar_op,0,NULL,1,targ_vecs,NULL);
00161 }
00162 
00163 void AbstractLinAlgPack::Vt_S( VectorMutable* v_lhs, const value_type& alpha )
00164 {
00165 #ifdef TEUCHOS_DEBUG
00166   TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vt_S(...), Error!");
00167 #endif
00168   if( alpha == 0.0 ) {
00169     *v_lhs = 0.0;
00170   }
00171   else if( alpha != 1.0 ) {
00172     TEST_FOR_EXCEPT(0!=RTOp_TOp_scale_vector_set_alpha( alpha, &scale_vector_op.op() ));
00173     VectorMutable* targ_vecs[1] = { v_lhs };
00174     apply_op(scale_vector_op,0,NULL,1,targ_vecs,NULL);
00175   }
00176 }
00177 
00178 void AbstractLinAlgPack::Vp_StV(
00179   VectorMutable* v_lhs, const value_type& alpha, const Vector& v_rhs)
00180 {
00181 #ifdef TEUCHOS_DEBUG
00182   TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vp_StV(...), Error!");
00183 #endif
00184   TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_set_alpha( alpha, &axpy_op.op() ));
00185   const Vector*  vecs[1]      = { &v_rhs };
00186   VectorMutable* targ_vecs[1] = { v_lhs  };
00187   apply_op(axpy_op,1,vecs,1,targ_vecs,NULL);
00188 }
00189 
00190 void AbstractLinAlgPack::Vp_StV(
00191   VectorMutable* y, const value_type& a, const SpVectorSlice& sx )
00192 {
00193   Vp_V_assert_compatibility(y,sx);
00194   if( sx.nz() ) {
00195     VectorSpace::vec_mut_ptr_t
00196         x = y->space().create_member();
00197     x->set_sub_vector(sub_vec_view(sx));
00198     Vp_StV( y, a, *x );
00199   }
00200 }
00201 
00202 void AbstractLinAlgPack::ele_wise_prod(
00203   const value_type& alpha, const Vector& v_rhs1, const Vector& v_rhs2
00204   , VectorMutable* v_lhs )
00205 {
00206 #ifdef TEUCHOS_DEBUG
00207   TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"ele_wise_prod(...), Error");
00208 #endif
00209   TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_prod_set_alpha(alpha,&ele_wise_prod_op.op()));
00210   const Vector*   vecs[2]      = { &v_rhs1, &v_rhs2 };
00211   VectorMutable*  targ_vecs[1] = { v_lhs };
00212   apply_op(ele_wise_prod_op,2,vecs,1,targ_vecs,NULL);
00213 }
00214 
00215 void AbstractLinAlgPack::ele_wise_divide(
00216   const value_type& alpha, const Vector& v_rhs1, const Vector& v_rhs2
00217   , VectorMutable* v_lhs )
00218 {
00219 #ifdef TEUCHOS_DEBUG
00220   TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"ele_wise_divide(...), Error");
00221 #endif
00222   TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_divide_set_alpha(alpha,&ele_wise_divide_op.op()));
00223   const int num_vecs = 2;
00224   const Vector*   vecs[2]      = { &v_rhs1, &v_rhs2 };
00225   VectorMutable*  targ_vecs[1] = { v_lhs };
00226   apply_op(ele_wise_divide_op,2,vecs,1,targ_vecs,NULL);
00227 }
00228 
00229 void AbstractLinAlgPack::seed_random_vector_generator( unsigned int s )
00230 {
00231   random_vector_op.set_seed(s);
00232 }
00233 
00234 void AbstractLinAlgPack::random_vector( value_type l, value_type u, VectorMutable* v )
00235 {
00236 #ifdef TEUCHOS_DEBUG
00237   TEST_FOR_EXCEPTION(v==NULL,std::logic_error,"Vt_S(...), Error");
00238 #endif
00239   //TEST_FOR_EXCEPT(0!=RTOp_TOp_random_vector_set_bounds( l, u, &random_vector_op.op() ));
00240   random_vector_op.set_bounds(l,u);
00241   VectorMutable* targ_vecs[1] = { v };
00242   apply_op(random_vector_op,0,NULL,1,targ_vecs,NULL);
00243   
00244 }
00245 
00246 void AbstractLinAlgPack::sign(
00247   const Vector      &v
00248   ,VectorMutable    *z
00249   )
00250 {
00251   RTOpPack::RTOpC op;
00252   TEST_FOR_EXCEPT( !( 0==RTOp_TOp_sign_construct(&op.op()) ) );
00253   const Vector*   vecs[1]      = { &v };
00254   VectorMutable*  targ_vecs[1] = { z  };
00255   apply_op(op,1,vecs,1,targ_vecs,NULL);
00256 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines