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_VectorMutable.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 "AbstractLinAlgPack_VectorMutable.hpp"
43 #include "AbstractLinAlgPack_VectorMutableSubView.hpp"
44 #include "AbstractLinAlgPack_VectorSpace.hpp"
45 #include "RTOp_TOp_assign_scalar.h"
47 #include "RTOp_TOp_axpy.h"
49 #include "RTOpPack_RTOpC.hpp"
50 #include "Teuchos_Assert.hpp"
51 
52 namespace {
53 
54 // vector scalar assignment operator
55 static RTOpPack::RTOpC assign_scalar_op;
56 // vector assignment operator
57 static RTOpPack::RTOpC assign_vec_op;
58 // set element operator
59 static RTOpPack::RTOpC set_ele_op;
60 // set sub-vector operator
61 static RTOpPack::RTOpC set_sub_vector_op;
62 // axpy operator
63 static RTOpPack::RTOpC axpy_op;
64 
65 // Simple class for an object that will initialize the operator objects
66 class init_rtop_server_t {
67 public:
68  init_rtop_server_t() {
69  // Vector scalar assignment operator
70  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_construct(0.0,&assign_scalar_op.op()));
71  // Vector assignment operator
72  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_vectors_construct(&assign_vec_op.op()));
73  // Set sub-vector operator
74  RTOp_SparseSubVector spc_sub_vec;
75  RTOp_sparse_sub_vector_null(&spc_sub_vec);
76  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_set_sub_vector_construct(&spc_sub_vec,&set_sub_vector_op.op()));
77  // axpy operator
78  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_construct(0.0,&axpy_op.op()));
79  }
80 };
81 
82 // When the program starts, this object will be created and the RTOp_Server object will
83 // be initialized before main() gets underway!
84 init_rtop_server_t init_rtop_server;
85 
86 } // end namespace
87 
88 namespace AbstractLinAlgPack {
89 
90 // VectorMutable
91 
93 {
94  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_set_alpha(alpha,&assign_scalar_op.op()));
95  VectorMutable* targ_vecs[1] = { this };
96  AbstractLinAlgPack::apply_op(assign_scalar_op,0,NULL,1,targ_vecs,NULL);
97  return *this;
98 }
99 
101 {
102  if( dynamic_cast<const void*>(&vec) == dynamic_cast<const void*>(this) )
103  return *this; // Assignment to self!
104  const Vector* vecs[1] = { &vec };
105  VectorMutable* targ_vecs[1] = { this };
106  AbstractLinAlgPack::apply_op(assign_vec_op,1,vecs,1,targ_vecs,NULL);
107  return *this;
108 }
109 
111 {
112  return this->operator=(static_cast<const Vector&>(vec));
113 }
114 
115 void VectorMutable::set_ele( index_type i, value_type alpha )
116 {
117  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_set_alpha(alpha,&assign_scalar_op.op()));
118  VectorMutable* targ_vecs[1] = { this };
119  AbstractLinAlgPack::apply_op(
120  assign_scalar_op,0,NULL,1,targ_vecs,NULL
121  ,i,1,0 // first_ele, sub_dim, global_offset
122  );
123 }
124 
127 {
128  namespace rcp = MemMngPack;
129  const index_type dim = this->dim();
130  const Range1D rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
131 #ifdef TEUCHOS_DEBUG
133  rng.ubound() > dim, std::out_of_range
134  ,"VectorMutable::sub_view(rng): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()<<"] "
135  "is not in the range [1,this->dim()] = [1,"<<dim<<"]" );
136 #endif
137  if( rng.lbound() == 1 && rng.ubound() == dim )
138  return Teuchos::rcp( this, false );
139  // We are returning a view that could change this vector so we had better
140  // wipe out the cache
141  //this->has_changed(); // I don't think this line is needed!
142  return Teuchos::rcp(
144  Teuchos::rcp( this, false )
145  ,rng ) );
146 }
147 
149 {
150  this->operator=(0.0);
151 }
152 
153 void VectorMutable::axpy( value_type alpha, const Vector& x )
154 {
155  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_set_alpha(alpha,&axpy_op.op()));
156  const Vector* vecs[1] = { &x };
157  VectorMutable* targ_vecs[1] = { this };
158  AbstractLinAlgPack::apply_op(axpy_op,1,vecs,1,targ_vecs,NULL);
159 }
160 
161 void VectorMutable::get_sub_vector( const Range1D& rng, RTOpPack::MutableSubVector* sub_vec_inout )
162 {
163  //
164  // Here we get a copy of the data for the sub-vector that the
165  // client will modify. We must later commit these changes to the
166  // actual vector when the client calls commitDetachedView(...).
167  // Note, this implementation is very dependent on the behavior of
168  // the default implementation of constant version of
169  // Vector<Scalar>::get_sub_vector(...) and the implementation of
170  // Vector<Scalar>::set_sub_vector(...)!
171  //
172  RTOpPack::SubVector sub_vec;
173  Vector::get_sub_vector( rng, &sub_vec );
174  sub_vec_inout->initialize(
175  sub_vec.globalOffset(),sub_vec.subDim(),const_cast<value_type*>(sub_vec.values()),sub_vec.stride());
176 }
177 
178 void VectorMutable::commit_sub_vector( RTOpPack::MutableSubVector* sub_vec_inout )
179 {
180  RTOpPack::SparseSubVector spc_sub_vec(
181  sub_vec_inout->globalOffset(), sub_vec_inout->subDim(),
182  Teuchos::arcp(sub_vec_inout->values(), 0, sub_vec_inout->stride()*sub_vec_inout->subDim(), false),
183  sub_vec_inout->stride()
184  );
185  VectorMutable::set_sub_vector( spc_sub_vec ); // Commit the changes!
186  RTOpPack::SubVector sub_vec(*sub_vec_inout);
187  Vector::free_sub_vector( &sub_vec ); // Free the memory!
188  sub_vec_inout->set_uninitialized(); // Make null as promised!
189 }
190 
191 void VectorMutable::set_sub_vector( const RTOpPack::SparseSubVector& sub_vec )
192 {
193  RTOp_SparseSubVector spc_sub_vec;
194  if (!is_null(sub_vec.indices())) {
195  RTOp_sparse_sub_vector(
196  sub_vec.globalOffset(), sub_vec.subDim(), sub_vec.subNz(),
197  sub_vec.values().get(), sub_vec.valuesStride(),
198  sub_vec.indices().get(), sub_vec.indicesStride(),
199  sub_vec.localOffset(), sub_vec.isSorted(),
200  &spc_sub_vec
201  );
202  }
203  else {
204  RTOp_SubVector _sub_vec;
206  sub_vec.globalOffset(), sub_vec.subDim(),
207  sub_vec.values().get(), sub_vec.valuesStride(),
208  &_sub_vec
209  );
210  RTOp_sparse_sub_vector_from_dense( &_sub_vec, &spc_sub_vec );
211  }
212  RTOpPack::RTOpC set_sub_vector_op;
213  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_set_sub_vector_construct(&spc_sub_vec,&set_sub_vector_op.op()));
214  VectorMutable* targ_vecs[1] = { this };
215  AbstractLinAlgPack::apply_op(
216  set_sub_vector_op,0,NULL,1,targ_vecs,NULL
217  ,sub_vec.globalOffset()+1,sub_vec.subDim(),sub_vec.globalOffset() // first_ele, sub_dim, global_offset
218  );
219 }
220 
222  value_type alpha
223  ,const GenPermMatrixSlice &P
224  ,BLAS_Cpp::Transp P_trans
225  ,const Vector &x
226  ,value_type beta
227  )
228 {
230  true, std::logic_error
231  ,"VectorMutable::Vp_StMtV(...): Error, this function has not yet been "
232  "given a default implementation and has not been overridden for the "
233  "subclass \'" << typeName(*this) << "\'!"
234  );
235  // ToDo: Implement using reduction or transformation operators that will
236  // work with any type of vector.
237 }
238 
239 // Overridden from Vector
240 
242 VectorMutable::sub_view( const Range1D& rng ) const
243 {
244  namespace rcp = MemMngPack;
245  return const_cast<VectorMutable*>(this)->sub_view(rng);
246 }
247 
248 } // end namespace AbstractLinAlgPack
virtual vec_mut_ptr_t sub_view(const Range1D &rng)
Create a mutable abstract view of a vector object.
virtual void set_sub_vector(const RTOpPack::SparseSubVector &sub_vec)
Set a specific sub-vector.
virtual void get_sub_vector(const Range1D &rng, RTOpPack::SubVector *sub_vec) const
Get a non-mutable explicit view of a sub-vector.
bool is_null(const boost::shared_ptr< T > &p)
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Index ubound() const
virtual void commit_sub_vector(RTOpPack::MutableSubVector *sub_vec)
Free a mutable explicit view of a sub-vector.
void RTOp_sub_vector(RTOp_index_type global_offset, RTOp_index_type sub_dim, const RTOp_value_type values[], ptrdiff_t values_stride, struct RTOp_SubVector *sub_vec)
virtual void set_ele(index_type i, value_type val)
Set a specific element of a vector.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool full_range() const
virtual void axpy(value_type alpha, const Vector &x)
Adds a linear combination of another vector to this vector object.
Concrete subclass for a sub-view of a VectorMutable object.
virtual void free_sub_vector(RTOpPack::SubVector *sub_vec) const
Free an explicit view of a sub-vector.
RTOp_RTOp & op()
Index lbound() const
virtual index_type dim() const
Return the dimension of this vector.
virtual VectorMutable & operator=(value_type alpha)
Assign the elements of this vector to a scalar.
virtual void get_sub_vector(const Range1D &rng, RTOpPack::MutableSubVector *sub_vec)
Get a mutable explicit view of a sub-vector.
Abstract interface for mutable coordinate vectors {abstract}.
virtual void Vp_StMtV(value_type alpha, const GenPermMatrixSlice &P, BLAS_Cpp::Transp P_trans, const Vector &x, value_type beta)
Perform a gather or scatter operation with a vector.
Transp
int RTOp_TOp_set_sub_vector_construct(const struct RTOp_SparseSubVector *sub_vec, struct RTOp_RTOp *op)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Concrete matrix type to represent general permutation (mapping) matrices.
std::string typeName(const T &t)