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_Vector.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 
43 #include "AbstractLinAlgPack_VectorMutable.hpp"
44 #include "AbstractLinAlgPack_VectorSubView.hpp"
45 #include "RTOp_ROp_dot_prod.h"
46 #include "RTOp_ROp_sum.h"
47 #include "RTOp_ROp_norms.h"
48 #include "RTOp_ROp_num_nonzeros.h"
50 #include "RTOpPack_RTOpC.hpp"
51 #include "RTOpPack_print_sub_vector.hpp"
52 #include "Teuchos_dyn_cast.hpp"
53 #include "Teuchos_Assert.hpp"
54 #include "Teuchos_FancyOStream.hpp"
55 
56 #include <limits>
57 #include <ostream>
58 
59 #include <assert.h>
60 
61 
62 // Uncomment to ignore cache of reduction data
63 //#define ALAP_VECTOR_IGNORE_CACHE_DATA
64 
65 
66 namespace {
67 
68 
69 // get element operator
70 static RTOpPack::RTOpC sum_op;
72 // number nonzros
73 static RTOpPack::RTOpC num_nonzeros_op;
74 static Teuchos::RCP<RTOpPack::ReductTarget> num_nonzeros_targ;
75 // Norm 1
76 static RTOpPack::RTOpC norm_1_op;
77 static Teuchos::RCP<RTOpPack::ReductTarget> norm_1_targ;
78 // Norm 2
79 static RTOpPack::RTOpC norm_2_op;
80 static Teuchos::RCP<RTOpPack::ReductTarget> norm_2_targ;
81 // Norm inf
82 static RTOpPack::RTOpC norm_inf_op;
83 static Teuchos::RCP<RTOpPack::ReductTarget> norm_inf_targ;
84 // get sub-vector operator
85 static RTOpPack::RTOpC get_sub_vector_op;
86 
87 
88 // Simple class for an object that will initialize the RTOp_Server and operators.
89 class init_rtop_server_t {
90 public:
91  init_rtop_server_t() {
92  // Operator and target for getting a vector element
93  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_sum_construct(&sum_op.op()));
94  sum_targ = sum_op.reduct_obj_create();
95  // Operator and target for norm 1
96  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_num_nonzeros_construct(&num_nonzeros_op.op()));
97  num_nonzeros_targ = num_nonzeros_op.reduct_obj_create();
98  // Operator and target for norm 1
99  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_norm_1_construct(&norm_1_op.op()));
100  norm_1_targ = norm_1_op.reduct_obj_create();
101  // Operator and target for norm 1
102  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_norm_2_construct(&norm_2_op.op()));
103  norm_2_targ = norm_2_op.reduct_obj_create();
104  // Operator and target for norm 1
105  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_norm_inf_construct(&norm_inf_op.op()));
106  norm_inf_targ = norm_inf_op.reduct_obj_create();
107  // Get sub-vector operator
108  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_get_sub_vector_construct(1,1,&get_sub_vector_op.op()));
109  }
110 };
111 
112 
113 // When the program starts, this object will be created
114 init_rtop_server_t init_rtop_server;
115 
116 
117 } // end namespace
118 
119 
120 namespace AbstractLinAlgPack {
121 
122 
124  :num_nonzeros_(-1), norm_1_(-1.0), norm_2_(-1.0), norm_inf_(-1.0) // uninitalized
125 {}
126 
127 
128 index_type Vector::dim() const
129 {
130  return this->space().dim();
131 }
132 
133 
134 index_type Vector::nz() const
135 {
136 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
137  {
138 #else
139  if( num_nonzeros_ < 0 ) {
140 #endif
141  num_nonzeros_op.reduct_obj_reinit(num_nonzeros_targ.ptr());
142  const Vector *vecs[1] = { this };
143  AbstractLinAlgPack::apply_op(num_nonzeros_op,1,vecs,0,NULL,&*num_nonzeros_targ);
144  num_nonzeros_ = RTOp_ROp_num_nonzeros_val(num_nonzeros_op(*num_nonzeros_targ));
145  }
146  return num_nonzeros_;
147 }
148 
149 
150 std::ostream& Vector::output(
151  std::ostream& out_arg, bool print_dim, bool newline,
152  index_type global_offset
153  ) const
154 {
155  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
156  Teuchos::OSTab tab(out);
157  RTOpPack::SubVector sub_vec;
158  this->get_sub_vector( Range1D(), &sub_vec );
159  RTOpPack::SubVector sub_vec_print( sub_vec.globalOffset() + global_offset, sub_vec.subDim(), sub_vec.values(), sub_vec.stride() );
160  RTOpPack::output(*out,sub_vec_print,print_dim,newline);
161  this->free_sub_vector( &sub_vec );
162  return out_arg;
163 }
164 
165 
167 {
169  vec = this->space().create_member();
170  *vec = *this;
171  return vec;
172 }
173 
174 
175 value_type Vector::get_ele(index_type i) const
176 {
177  sum_op.reduct_obj_reinit(sum_targ.ptr());
178  const Vector *vecs[1] = { this };
179  AbstractLinAlgPack::apply_op(
180  sum_op,1,vecs,0,NULL,&*sum_targ
181  ,i,1,0 // first_ele, sub_dim, global_offset
182  );
183  return RTOp_ROp_sum_val(sum_op(*sum_targ));
184 }
185 
186 
187 value_type Vector::norm_1() const
188 {
189 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
190  {
191 #else
192  if( norm_1_ < 0.0 ) {
193 #endif
194  norm_1_op.reduct_obj_reinit(norm_1_targ.ptr());
195  const Vector *vecs[1] = { this };
196  AbstractLinAlgPack::apply_op(norm_1_op,1,vecs,0,NULL,&*norm_1_targ);
197  norm_1_ = RTOp_ROp_norm_1_val(norm_1_op(*norm_1_targ));
198  }
199  return norm_1_;
200 }
201 
202 
203 value_type Vector::norm_2() const
204 {
205 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
206  {
207 #else
208  if( norm_2_ < 0.0 ) {
209 #endif
210  norm_2_op.reduct_obj_reinit(norm_2_targ.ptr());
211  const Vector *vecs[1] = { this };
212  AbstractLinAlgPack::apply_op(norm_2_op,1,vecs,0,NULL,&*norm_2_targ);
213  norm_2_ = RTOp_ROp_norm_2_val(norm_2_op(*norm_2_targ));
214  }
215  return norm_2_;
216 }
217 
218 
219 value_type Vector::norm_inf() const
220 {
221 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
222  {
223 #else
224  if( norm_inf_ < 0.0 ) {
225 #endif
226  norm_inf_op.reduct_obj_reinit(norm_inf_targ.ptr());
227  const Vector *vecs[1] = { this };
228  AbstractLinAlgPack::apply_op(norm_inf_op,1,vecs,0,NULL,&*norm_inf_targ);
229  norm_inf_ = RTOp_ROp_norm_inf_val(norm_inf_op(*norm_inf_targ));
230  }
231  return norm_inf_;
232 }
233 
234 
235 value_type Vector::inner_product( const Vector& v ) const
236 {
237  return this->space().inner_prod()->inner_prod(*this,v);
238 }
239 
241 Vector::sub_view( const Range1D& rng_in ) const
242 {
243  namespace rcp = MemMngPack;
244  const index_type dim = this->dim();
245  const Range1D rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
246 #ifdef TEUCHOS_DEBUG
248  rng.ubound() > dim, std::out_of_range
249  ,"Vector::sub_view(rng): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()<<"] "
250  "is not in the range [1,this->dim()] = [1,"<<dim<<"]" );
251 #endif
252  if( rng.lbound() == 1 && rng.ubound() == dim )
253  return vec_ptr_t( this, false );
254  return Teuchos::rcp(
255  new VectorSubView(
256  vec_ptr_t( this, false )
257  ,rng ) );
258 }
259 
260 
261 void Vector::get_sub_vector( const Range1D& rng_in, RTOpPack::SubVector* sub_vec_inout ) const
262 {
263  const Range1D rng = rng_in.full_range() ? Range1D(1,this->space().dim()) : rng_in;
264 #ifdef TEUCHOS_DEBUG
266  this->space().dim() < rng.ubound(), std::out_of_range
267  ,"Vector::get_sub_vector(rng,...): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()
268  <<"] is not in range = [1,"<<this->space().dim()<<"]" );
269 #endif
270  // Free sub_vec if needed (note this is dependent on the implemenation of this operator class!)
271  if( sub_vec_inout->values() ) {
272  std::free( (void*)sub_vec_inout->values() );
273  }
274  // Initialize the operator
275  RTOpPack::RTOpC get_sub_vector_op;
276  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_get_sub_vector_construct(rng.lbound(),rng.ubound(),&get_sub_vector_op.op()));
277  // Create the reduction object (another sub_vec)
278  Teuchos::RCP<RTOpPack::ReductTarget> reduct_obj = get_sub_vector_op.reduct_obj_create(); // This is really of type RTOpPack::ConstSubVectorView<Scalar>!
279  // Perform the reduction (get the sub-vector requested)
280  const size_t num_vecs = 1;
281  const Vector* sub_vecs[num_vecs] = { this };
282  AbstractLinAlgPack::apply_op(
283  get_sub_vector_op,num_vecs,sub_vecs,0,NULL,&*reduct_obj
284  ,rng.lbound(),rng.size(),rng.lbound()-1 // first_ele, sub_dim, global_offset
285  );
286  // Set the sub-vector. Note reduct_obj will go out of scope so the sub_vec parameter will
287  // own the memory allocated within get_sub_vector_op.create_reduct_obj_raw(...). This is okay
288  // since the client is required to call release_sub_vector(...) so release memory!
289  RTOp_ReductTarget reduct_obj_raw = get_sub_vector_op(*reduct_obj);
290  RTOp_SubVector sub_vec = RTOp_ROp_get_sub_vector_val(reduct_obj_raw);
291  sub_vec_inout->initialize(sub_vec.global_offset,sub_vec.sub_dim,sub_vec.values,sub_vec.values_stride);\
292  reduct_obj.release(); // Do not allow delete to be called!
293  std::free(reduct_obj_raw); // Now *sub_vec owns the values[] and indices[] arrays!
294 }
295 
296 
297 void Vector::free_sub_vector( RTOpPack::SubVector* sub_vec ) const
298 {
299  // Free sub_vec if needed (note this is dependent on the implemenation of this operator class!)
300  if( sub_vec->values() )
301  std::free( (void*)sub_vec->values() );
302  sub_vec->set_uninitialized();
303 }
304 
305 
307 {
308  num_nonzeros_= -1; // uninitalized;
309  norm_1_ = norm_2_ = norm_inf_ = -1.0;
310 }
311 
312 
313 // protected
314 
315 
317  const size_t num_targ_vecs, VectorMutable** targ_vecs
318  ) const
319 {
320  for( size_t k = 0; k < num_targ_vecs; ++k )
321  targ_vecs[k]->has_changed();
322 }
323 
324 
325 } // end namespace AbstractLinAlgPack
326 
327 
328 // nonmember functions
329 
330 
331 void AbstractLinAlgPack::apply_op(
332  const RTOpPack::RTOp &op
333  ,const size_t num_vecs
334  ,const Vector* vecs[]
335  ,const size_t num_targ_vecs
336  ,VectorMutable* targ_vecs[]
337  ,RTOpPack::ReductTarget *reduct_obj
338  ,const index_type first_ele
339  ,const index_type sub_dim
340  ,const index_type global_offset
341  )
342 {
343  if(num_vecs)
344  vecs[0]->apply_op(op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj,first_ele,sub_dim,global_offset);
345  else if (num_targ_vecs)
346  targ_vecs[0]->apply_op(op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj,first_ele,sub_dim,global_offset);
347 }
virtual vec_ptr_t sub_view(const Range1D &rng) const
Create an abstract view of a vector object .
virtual const VectorSpace & space() const =0
Return the vector space that this vector belongs to.
virtual void get_sub_vector(const Range1D &rng, RTOpPack::SubVector *sub_vec) const
Get a non-mutable explicit view of a sub-vector.
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
virtual value_type norm_1() const
One norm. ||v||_1 = sum( |v(i)|, i = 1,,,this->dim() )
virtual void finalize_apply_op(const size_t num_targ_vecs, VectorMutable **targ_vecs) const
This method usually needs to be called by subclasses at the end of the apply_op() method implementati...
Concrete subclass for a default sub-view implementation for a Vector object.
virtual vec_mut_ptr_t clone() const
Create a clone of this vector objet.
virtual value_type norm_inf() const
Infinity norm. ||v||_inf = max( |v(i)|, i = 1,,,this->dim() )
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Index size() const
Teuchos::RCP< const Vector > vec_ptr_t
Index ubound() const
virtual value_type norm_2() const
Two norm. ||v||_2 = sqrt( sum( v(i)^2, i = 1,,,this->dim() ) )
struct RTOp_SubVector RTOp_ROp_get_sub_vector_val(RTOp_ReductTarget reduct_obj)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void has_changed() const
Must be called by any vector subclass that modifies this vector object!
virtual std::ostream & output(std::ostream &out, bool print_dim=true, bool newline=true, index_type global_offset=0) const
Virtual output function.
bool full_range() const
virtual void free_sub_vector(RTOpPack::SubVector *sub_vec) const
Free an explicit view of a sub-vector.
RTOp_RTOp & op()
Teuchos::RCP< ReductTarget > reduct_obj_create() const
virtual index_type dim() const =0
Return the dimmension of the vector space.
Index lbound() const
virtual index_type dim() const
Return the dimension of this vector.
virtual index_type nz() const
Return the number of nonzero elements in the vector.
virtual value_type get_ele(index_type i) const
Fetch an element in the vector.
Abstract interface for mutable coordinate vectors {abstract}.
virtual vec_mut_ptr_t create_member() const =0
Create a vector member from the vector space.
virtual value_type inner_product(const Vector &v) const
Return the inner product of *this with v.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual value_type inner_prod(const Vector &v1, const Vector &v2) const =0
Compute the inner product of two vectors.