ROL
ROL_PartitionedVector.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) 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 lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #include "ROL_Vector.hpp"
45 
46 #ifndef ROL_PARTITIONED_VECTOR_H
47 #define ROL_PARTITIONED_VECTOR_H
48 
55 namespace ROL {
56 
57 template<class Real>
58 class PartitionedVector : public Vector<Real> {
59 
60  typedef Vector<Real> V;
61  typedef Teuchos::RCP<V> RCPV;
63 
64 private:
65  const std::vector<RCPV> vecs_;
66  mutable std::vector<RCPV> dual_vecs_;
67  mutable Teuchos::RCP<PV> dual_pvec_;
68 public:
69 
70  typedef typename std::vector<PV>::size_type size_type;
71 
72  PartitionedVector( const std::vector<RCPV> &vecs ) : vecs_(vecs) {
73  for( size_type i=0; i<vecs_.size(); ++i ) {
74  dual_vecs_.push_back((vecs_[i]->dual()).clone());
75  }
76  }
77 
78  void set( const V &x ) {
79  using Teuchos::dyn_cast;
80  const PV &xs = dyn_cast<const PV>(dyn_cast<const V>(x));
81 
82  TEUCHOS_TEST_FOR_EXCEPTION( numVectors() != xs.numVectors(),
83  std::invalid_argument,
84  "Error: Vectors must have the same number of subvectors." );
85 
86  for( size_type i=0; i<vecs_.size(); ++i ) {
87  vecs_[i]->set(*xs.get(i));
88  }
89  }
90 
91  void plus( const V &x ) {
92  using Teuchos::dyn_cast;
93  const PV &xs = dyn_cast<const PV>(dyn_cast<const V>(x));
94 
95  TEUCHOS_TEST_FOR_EXCEPTION( numVectors() != xs.numVectors(),
96  std::invalid_argument,
97  "Error: Vectors must have the same number of subvectors." );
98 
99  for( size_type i=0; i<vecs_.size(); ++i ) {
100  vecs_[i]->plus(*xs.get(i));
101  }
102  }
103 
104  void scale( const Real alpha ) {
105  for( size_type i=0; i<vecs_.size(); ++i ) {
106  vecs_[i]->scale(alpha);
107  }
108  }
109 
110  void axpy( const Real alpha, const V &x ) {
111  using Teuchos::dyn_cast;
112  const PV &xs = dyn_cast<const PV>(x);
113 
114  TEUCHOS_TEST_FOR_EXCEPTION( numVectors() != xs.numVectors(),
115  std::invalid_argument,
116  "Error: Vectors must have the same number of subvectors." );
117 
118  for( size_type i=0; i<vecs_.size(); ++i ) {
119  vecs_[i]->axpy(alpha,*xs.get(i));
120  }
121  }
122 
123  Real dot( const V &x ) const {
124  using Teuchos::dyn_cast;
125  const PV &xs = dyn_cast<const PV>(x);
126 
127  TEUCHOS_TEST_FOR_EXCEPTION( numVectors() != xs.numVectors(),
128  std::invalid_argument,
129  "Error: Vectors must have the same number of subvectors." );
130 
131  Real result = 0;
132  for( size_type i=0; i<vecs_.size(); ++i ) {
133  result += vecs_[i]->dot(*xs.get(i));
134  }
135  return result;
136  }
137 
138  Real norm() const {
139  Real result = 0;
140  for( size_type i=0; i<vecs_.size(); ++i ) {
141  result += std::pow(vecs_[i]->norm(),2);
142  }
143  return std::sqrt(result);
144  }
145 
146  RCPV clone() const {
147  using Teuchos::RCP;
148  using Teuchos::rcp;
149 
150  std::vector<RCPV> clonevec;
151  for( size_type i=0; i<vecs_.size(); ++i ) {
152  clonevec.push_back(vecs_[i]->clone());
153  }
154  return rcp( new PV(clonevec) );
155  }
156 
157  const V& dual(void) const {
158  using Teuchos::rcp;
159 
160  for( size_type i=0; i<vecs_.size(); ++i ) {
161  dual_vecs_[i]->set(vecs_[i]->dual());
162  }
163  dual_pvec_ = rcp( new PV( dual_vecs_ ) );
164  return *dual_pvec_;
165  }
166 
167  RCPV basis( const int i ) const {
168 
169  TEUCHOS_TEST_FOR_EXCEPTION( i >= dimension() || i<0,
170  std::invalid_argument,
171  "Error: Basis index must be between 0 and vector dimension." );
172 
173  using Teuchos::RCP;
174  using Teuchos::rcp;
175  using Teuchos::dyn_cast;
176 
177  RCPV bvec = clone();
178 
179  // Downcast
180  PV &eb = dyn_cast<PV>(*bvec);
181 
182  int begin = 0;
183  int end = 0;
184 
185  // Iterate over subvectors
186  for( size_type j=0; j<vecs_.size(); ++j ) {
187 
188  end += vecs_[j]->dimension();
189 
190  if( begin<= i && i<end ) {
191  eb.set(j, *(vecs_[j]->basis(i-begin)) );
192  }
193  else {
194  eb.zero(j);
195  }
196 
197  begin = end;
198 
199  }
200  return bvec;
201  }
202 
203  int dimension() const {
204  int total_dim = 0;
205  for( size_type j=0; j<vecs_.size(); ++j ) {
206  total_dim += vecs_[j]->dimension();
207  }
208  return total_dim;
209  }
210 
211  void zero() {
212  for( size_type j=0; j<vecs_.size(); ++j ) {
213  vecs_[j]->zero();
214  }
215  }
216 
217  // Apply the same unary function to each subvector
218  void applyUnary( const Elementwise::UnaryFunction<Real> &f ) {
219  for( size_type i=0; i<vecs_.size(); ++i ) {
220  vecs_[i]->applyUnary(f);
221  }
222  }
223 
224  // Apply the same binary function to each pair of subvectors in this vector and x
225  void applyBinary( const Elementwise::BinaryFunction<Real> &f, const V &x ) {
226  const PV &xs = Teuchos::dyn_cast<const PV>(x);
227 
228  for( size_type i=0; i<vecs_.size(); ++i ) {
229  vecs_[i]->applyBinary(f,*xs.get(i));
230  }
231  }
232 
233  Real reduce( const Elementwise::ReductionOp<Real> &r ) const {
234  Real result = r.initialValue();
235 
236  for( size_type i=0; i<vecs_.size(); ++i ) {
237  r.reduce(vecs_[i]->reduce(r),result);
238  }
239  return result;
240  }
241 
242  // Methods that do not exist in the base class
243 
244  Teuchos::RCP<const Vector<Real> > get(size_type i) const {
245  return vecs_[i];
246  }
247 
248  Teuchos::RCP<Vector<Real> > get(size_type i) {
249  return vecs_[i];
250  }
251 
252  void set(size_type i, const V &x) {
253  vecs_[i]->set(x);
254  }
255 
256  void zero(size_type i) {
257  vecs_[i]->zero();
258  }
259 
261  return vecs_.size();
262  }
263 
264 };
265 
266 // Helper methods
267 template<class Real>
268 Teuchos::RCP<Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<Vector<Real> > &a ) {
269  using Teuchos::RCP;
270  using Teuchos::rcp;
271  typedef RCP<Vector<Real> > RCPV;
272  typedef PartitionedVector<Real> PV;
273 
274  RCPV temp[] = {a};
275  return rcp( new PV( std::vector<RCPV>(temp, temp+1) ) );
276 }
277 
278 template<class Real>
279 Teuchos::RCP<const Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<const Vector<Real> > &a ) {
280  using Teuchos::RCP;
281  using Teuchos::rcp;
282  typedef RCP<const Vector<Real> > RCPV;
283  typedef const PartitionedVector<Real> PV;
284 
285  RCPV temp[] = {a};
286  return rcp( new PV( std::vector<RCPV>(temp, temp+1) ) );
287 }
288 
289 template<class Real>
290 Teuchos::RCP<Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<Vector<Real> > &a,
291  const Teuchos::RCP<Vector<Real> > &b ) {
292  using Teuchos::RCP;
293  using Teuchos::rcp;
294  typedef RCP<Vector<Real> > RCPV;
295  typedef PartitionedVector<Real> PV;
296 
297  RCPV temp[] = {a,b};
298  return rcp( new PV( std::vector<RCPV>(temp, temp+2) ) );
299 }
300 
301 template<class Real>
302 Teuchos::RCP<const Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<const Vector<Real> > &a,
303  const Teuchos::RCP<const Vector<Real> > &b ) {
304  using Teuchos::RCP;
305  using Teuchos::rcp;
306  typedef RCP<const Vector<Real> > RCPV;
307  typedef const PartitionedVector<Real> PV;
308 
309  RCPV temp[] = {a,b};
310  return rcp( new PV( std::vector<RCPV>(temp, temp+2) ) );
311 }
312 
313 template<class Real>
314 Teuchos::RCP<Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<Vector<Real> > &a,
315  const Teuchos::RCP<Vector<Real> > &b,
316  const Teuchos::RCP<Vector<Real> > &c ) {
317  using Teuchos::RCP;
318  using Teuchos::rcp;
319  typedef RCP<Vector<Real> > RCPV;
320  typedef PartitionedVector<Real> PV;
321 
322  RCPV temp[] = {a,b,c};
323  return rcp( new PV( std::vector<RCPV>(temp, temp+3) ) );
324 }
325 
326 template<class Real>
327 Teuchos::RCP<const Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<const Vector<Real> > &a,
328  const Teuchos::RCP<const Vector<Real> > &b,
329  const Teuchos::RCP<const Vector<Real> > &c ) {
330  using Teuchos::RCP;
331  using Teuchos::rcp;
332  typedef RCP<const Vector<Real> > RCPV;
333  typedef const PartitionedVector<Real> PV;
334 
335  RCPV temp[] = {a,b,c};
336  return rcp( new PV( std::vector<RCPV>(temp, temp+3) ) );
337 }
338 
339 template<class Real>
340 Teuchos::RCP<Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<Vector<Real> > &a,
341  const Teuchos::RCP<Vector<Real> > &b,
342  const Teuchos::RCP<Vector<Real> > &c,
343  const Teuchos::RCP<Vector<Real> > &d ) {
344  using Teuchos::RCP;
345  using Teuchos::rcp;
346  typedef RCP<Vector<Real> > RCPV;
347  typedef PartitionedVector<Real> PV;
348 
349  RCPV temp[] = {a,b,c,d};
350  return rcp( new PV( std::vector<RCPV>(temp, temp+4) ) );
351 }
352 
353 template<class Real>
354 Teuchos::RCP<const Vector<Real> > CreatePartitionedVector( const Teuchos::RCP<const Vector<Real> > &a,
355  const Teuchos::RCP<const Vector<Real> > &b,
356  const Teuchos::RCP<const Vector<Real> > &c,
357  const Teuchos::RCP<const Vector<Real> > &d ) {
358  using Teuchos::RCP;
359  using Teuchos::rcp;
360  typedef RCP<const Vector<Real> > RCPV;
361  typedef const PartitionedVector<Real> PV;
362 
363  RCPV temp[] = {a,b,c,d};
364  return rcp( new PV( std::vector<RCPV>(temp, temp+4) ) );
365 }
366 
367 } // namespace ROL
368 
369 #endif // ROL_PARTITIONED_VECTOR_H
370 
Real norm() const
Returns where .
Defines the linear algebra of vector space on a generic partitioned vector.
const std::vector< RCPV > vecs_
Teuchos::RCP< const Vector< Real > > get(size_type i) const
RCPV basis(const int i) const
Return i-th basis vector.
void applyBinary(const Elementwise::BinaryFunction< Real > &f, const V &x)
Real dot(const V &x) const
Compute where .
Teuchos::RCP< Vector< Real > > CreatePartitionedVector(const Teuchos::RCP< Vector< Real > > &a)
const V & dual(void) const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
int dimension() const
Return dimension of the vector space.
void zero()
Set to zero vector.
void applyUnary(const Elementwise::UnaryFunction< Real > &f)
PartitionedVector< Real > PV
PartitionedVector(const std::vector< RCPV > &vecs)
void scale(const Real alpha)
Compute where .
std::vector< RCPV > dual_vecs_
void set(const V &x)
Set where .
Real reduce(const Elementwise::ReductionOp< Real > &r) const
std::vector< PV >::size_type size_type
RCPV clone() const
Clone to make a new (uninitialized) vector.
void plus(const V &x)
Compute , where .
void axpy(const Real alpha, const V &x)
Compute where .