AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_VectorStdOps.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include "AbstractLinAlgPack_VectorStdOps.hpp"
45 #include "AbstractLinAlgPack_VectorSpace.hpp"
46 #include "AbstractLinAlgPack_VectorMutable.hpp"
47 #include "AbstractLinAlgPack_AssertOp.hpp"
48 #include "AbstractLinAlgPack_SpVectorClass.hpp"
49 #include "AbstractLinAlgPack_SpVectorView.hpp"
50 #include "RTOp_ROp_dot_prod.h"
51 #include "RTOp_ROp_max_abs_ele.h"
52 #include "RTOp_ROp_sum.h"
53 #include "RTOp_TOp_add_scalar.h"
54 #include "RTOp_TOp_axpy.h"
56 #include "RTOp_TOp_ele_wise_prod.h"
57 //#include "RTOp_TOp_random_vector.h"
58 #include "RTOpPack_TOpRandomize.hpp"
59 #include "RTOp_TOp_scale_vector.h"
60 #include "RTOp_TOp_sign.h"
61 #include "RTOpPack_RTOpC.hpp"
62 #include "Teuchos_Assert.hpp"
63 
64 namespace {
65 
66 // sum
67 static RTOpPack::RTOpC sum_op;
69 // dot prod
70 static RTOpPack::RTOpC dot_prod_op;
71 static Teuchos::RCP<RTOpPack::ReductTarget> dot_prod_targ;
72 // number of bounded elements
73 static RTOpPack::RTOpC num_bounded_op;
74 static Teuchos::RCP<RTOpPack::ReductTarget> num_bounded_targ;
75 // add scalar to vector
76 static RTOpPack::RTOpC add_scalar_op;
77 // scale vector
78 static RTOpPack::RTOpC scale_vector_op;
79 // axpy
80 static RTOpPack::RTOpC axpy_op;
81 // random vector
82 //static RTOpPack::RTOpC random_vector_op;
84 // element-wise division
85 static RTOpPack::RTOpC ele_wise_divide_op;
86 // element-wise product
87 static RTOpPack::RTOpC ele_wise_prod_op;
88 
89 // Simple class for an object that will initialize the RTOp_Server.
90 class init_rtop_server_t {
91 public:
92  init_rtop_server_t() {
93  // Operator and target obj for sum
94  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_sum_construct(&sum_op.op()));
95  sum_targ = sum_op.reduct_obj_create();
96  // Operator and target obj for dot_prod
97  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_dot_prod_construct(&dot_prod_op.op()));
98  dot_prod_targ = dot_prod_op.reduct_obj_create();
99  // Operator add_scalar
100  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_add_scalar_construct(0.0,&add_scalar_op.op()));
101  // Operator scale_vector
102  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_scale_vector_construct(0.0,&scale_vector_op.op()));
103  // Operator axpy
104  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_construct(0.0,&axpy_op.op()));
105  // Operator random_vector
106  //TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_random_vector_construct(0.0,0.0,&random_vector_op.op()));
107  // Operator ele_wise_divide
108  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_divide_construct(0.0,&ele_wise_divide_op.op()));
109  // Operator ele_wise_prod
110  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_prod_construct(0.0,&ele_wise_prod_op.op()));
111  }
112 };
113 
114 // When the program starts, this object will be created and the RTOp_Server object will
115 // be initialized before main() gets underway!
116 init_rtop_server_t init_rtop_server;
117 
118 } // end namespace
119 
120 AbstractLinAlgPack::value_type
122 {
123  sum_op.reduct_obj_reinit(sum_targ.ptr());
124  const Vector* vecs[1] = { &v_rhs };
125  apply_op(sum_op,1,vecs,0,NULL,&*sum_targ);
126  return RTOp_ROp_sum_val(sum_op(*sum_targ));
127 }
128 
129 AbstractLinAlgPack::value_type
130 AbstractLinAlgPack::dot( const Vector& v_rhs1, const Vector& v_rhs2 )
131 {
132  dot_prod_op.reduct_obj_reinit(dot_prod_targ.ptr());
133  const Vector* vecs[2] = { &v_rhs1, &v_rhs2 };
134  apply_op(dot_prod_op,2,vecs,0,NULL,&*dot_prod_targ);
135  return RTOp_ROp_dot_prod_val(dot_prod_op(*dot_prod_targ));
136 }
137 
138 AbstractLinAlgPack::value_type
139 AbstractLinAlgPack::dot( const Vector& v_rhs1, const SpVectorSlice& sv_rhs2 )
140 {
141  VopV_assert_compatibility(v_rhs1,sv_rhs2 );
142  if( sv_rhs2.nz() ) {
144  v_rhs2 = v_rhs1.space().create_member();
145  v_rhs2->set_sub_vector(sub_vec_view(sv_rhs2));
146  return dot(v_rhs1,*v_rhs2);
147  }
148  return 0.0;
149 }
150 
152  const Vector& v, value_type* max_v_j, index_type* max_j
153  )
154 {
155  TEUCHOS_TEST_FOR_EXCEPT( !( max_v_j && max_j ) );
156  RTOpPack::RTOpC op;
157  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_max_abs_ele_construct(&op.op()));
159  const Vector* vecs[1] = { &v };
160  apply_op(op,1,vecs,0,NULL,&*reduct_obj);
161  RTOp_value_index_type val = RTOp_ROp_max_abs_ele_val(op(*reduct_obj));
162  *max_v_j = val.value;
163  *max_j = val.index;
164 }
165 
166 void AbstractLinAlgPack::Vp_S( VectorMutable* v_lhs, const value_type& alpha )
167 {
168 #ifdef TEUCHOS_DEBUG
169  TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vp_S(...), Error!");
170 #endif
171  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_add_scalar_set_alpha(alpha,&add_scalar_op.op()));
172  VectorMutable* targ_vecs[1] = { v_lhs };
173  apply_op(add_scalar_op,0,NULL,1,targ_vecs,NULL);
174 }
175 
176 void AbstractLinAlgPack::Vt_S( VectorMutable* v_lhs, const value_type& alpha )
177 {
178 #ifdef TEUCHOS_DEBUG
179  TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vt_S(...), Error!");
180 #endif
181  if( alpha == 0.0 ) {
182  *v_lhs = 0.0;
183  }
184  else if( alpha != 1.0 ) {
185  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_scale_vector_set_alpha( alpha, &scale_vector_op.op() ));
186  VectorMutable* targ_vecs[1] = { v_lhs };
187  apply_op(scale_vector_op,0,NULL,1,targ_vecs,NULL);
188  }
189 }
190 
192  VectorMutable* v_lhs, const value_type& alpha, const Vector& v_rhs)
193 {
194 #ifdef TEUCHOS_DEBUG
195  TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"Vp_StV(...), Error!");
196 #endif
197  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_set_alpha( alpha, &axpy_op.op() ));
198  const Vector* vecs[1] = { &v_rhs };
199  VectorMutable* targ_vecs[1] = { v_lhs };
200  apply_op(axpy_op,1,vecs,1,targ_vecs,NULL);
201 }
202 
204  VectorMutable* y, const value_type& a, const SpVectorSlice& sx )
205 {
207  if( sx.nz() ) {
209  x = y->space().create_member();
210  x->set_sub_vector(sub_vec_view(sx));
211  Vp_StV( y, a, *x );
212  }
213 }
214 
216  const value_type& alpha, const Vector& v_rhs1, const Vector& v_rhs2
217  , VectorMutable* v_lhs )
218 {
219 #ifdef TEUCHOS_DEBUG
220  TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"ele_wise_prod(...), Error");
221 #endif
222  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_prod_set_alpha(alpha,&ele_wise_prod_op.op()));
223  const Vector* vecs[2] = { &v_rhs1, &v_rhs2 };
224  VectorMutable* targ_vecs[1] = { v_lhs };
225  apply_op(ele_wise_prod_op,2,vecs,1,targ_vecs,NULL);
226 }
227 
229  const value_type& alpha, const Vector& v_rhs1, const Vector& v_rhs2
230  , VectorMutable* v_lhs )
231 {
232 #ifdef TEUCHOS_DEBUG
233  TEUCHOS_TEST_FOR_EXCEPTION(v_lhs==NULL,std::logic_error,"ele_wise_divide(...), Error");
234 #endif
235  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_ele_wise_divide_set_alpha(alpha,&ele_wise_divide_op.op()));
236  const int num_vecs = 2;
237  const Vector* vecs[2] = { &v_rhs1, &v_rhs2 };
238  VectorMutable* targ_vecs[1] = { v_lhs };
239  apply_op(ele_wise_divide_op,2,vecs,1,targ_vecs,NULL);
240 }
241 
243 {
244  random_vector_op.set_seed(s);
245 }
246 
247 void AbstractLinAlgPack::random_vector( value_type l, value_type u, VectorMutable* v )
248 {
249 #ifdef TEUCHOS_DEBUG
250  TEUCHOS_TEST_FOR_EXCEPTION(v==NULL,std::logic_error,"Vt_S(...), Error");
251 #endif
252  //TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_random_vector_set_bounds( l, u, &random_vector_op.op() ));
253  random_vector_op.set_bounds(l,u);
254  VectorMutable* targ_vecs[1] = { v };
255  apply_op(random_vector_op,0,NULL,1,targ_vecs,NULL);
256 
257 }
258 
260  const Vector &v
261  ,VectorMutable *z
262  )
263 {
264  RTOpPack::RTOpC op;
265  TEUCHOS_TEST_FOR_EXCEPT( !( 0==RTOp_TOp_sign_construct(&op.op()) ) );
266  const Vector* vecs[1] = { &v };
267  VectorMutable* targ_vecs[1] = { z };
268  apply_op(op,1,vecs,1,targ_vecs,NULL);
269 }
virtual const VectorSpace & space() const =0
Return the vector space that this vector belongs to.
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
void sign(const Vector &v, VectorMutable *z)
Compute the sign of each element in an input vector.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
void seed_random_vector_generator(unsigned int)
Seed the random number generator.
void set_bounds(const Scalar &l, const Scalar &u)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
void ele_wise_divide(const value_type &alpha, const Vector &v_rhs1, const Vector &v_rhs2, VectorMutable *v_lhs)
v_lhs(i) = alpha * v_rhs1(i) / v_rhs2(i), i = 1,,,dim.
void Vp_V_assert_compatibility(VectorMutable *v_lhs, const Vector &v_rhs)
v_lhs += op v_rhs
Ptr< T > ptr() const
void reduct_obj_reinit(const Ptr< ReductTarget > &reduct_obj) const
void ele_wise_prod(const value_type &alpha, const Vector &v_rhs1, const Vector &v_rhs2, VectorMutable *v_lhs)
v_lhs(i) += alpha * v_rhs1(i) * v_rhs2(i), i = 1,,,dim.
RTOp_RTOp & op()
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
Teuchos::RCP< ReductTarget > reduct_obj_create() const
Abstract interface for mutable coordinate vectors {abstract}.
void Vp_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs += alpha
void max_abs_ele(const Vector &v, value_type *max_v_j, index_type *max_j)
Compute the maximum element in a vector.
value_type sum(const Vector &v_rhs)
result = sum( v_rhs(i), i = 1,,,dim )
virtual vec_mut_ptr_t create_member() const =0
Create a vector member from the vector space.
void set_seed(const unsigned int seed)
void random_vector(value_type l, value_type u, VectorMutable *v)
Generate a random vector with elements uniformly distrubuted elements.
size_type nz() const
Return the number of non-zero elements.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void VopV_assert_compatibility(const Vector &v_rhs1, const Vector &v_rhs2)
v_rhs1 op v_rhs2