ROL
ROL_BlockOperator.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 #ifndef ROL_BLOCKOPERATOR_H
45 #define ROL_BLOCKOPERATOR_H
46 
47 #include "ROL_LinearOperator.hpp"
49 
58 namespace ROL {
59 
60 template<class Real>
61 class BlockOperator : public LinearOperator<Real> {
62 
63  typedef Vector<Real> V; // ROL Vector
64  typedef PartitionedVector<Real> PV; // Partitioned Vector
65  typedef LinearOperator<Real> OP; // Linear Operator
66 
67  typedef std::vector<Teuchos::RCP<OP> > OpVec; // Vector (column-stacked matrix) of pointers to operators
68  typedef typename OpVec::size_type uint; // index type
69 
70 private:
71 
72  Teuchos::RCP<OpVec> blocks_;
73 
74 public:
75 
76  BlockOperator( const Teuchos::RCP<OpVec> &blocks ) : blocks_(blocks) {}
77 
78  virtual void apply( V &Hv, const V &v, Real &tol ) const {
79 
80  // Downcast to Partitioned Vectors
81  PV &Hv_part = Teuchos::dyn_cast<PV>(Hv);
82  const PV &v_part = Teuchos::dyn_cast<const PV>(v);
83 
84  uint nvec1 = v_part.numVectors();
85  uint nvec2 = Hv_part.numVectors();
86  uint nblks = blocks_->size();
87 
88  TEUCHOS_TEST_FOR_EXCEPTION( (nvec1 != nvec2), std::invalid_argument,
89  ">>> ERROR (ROL_BlockOperator, apply): "
90  "Mismatch between input and output number of subvectors.");
91 
92  TEUCHOS_TEST_FOR_EXCEPTION( (nblks != nvec1*nvec2 ) , std::invalid_argument,
93  ">>> ERROR (ROL_BlockOperator, apply): "
94  "Block operator dimension mismatch.");
95 
96  for( uint i=0; i<nvec1; ++i ) {
97 
98  Teuchos::RCP<V> Hvi = Hv_part.get(i);
99  Teuchos::RCP<V> u = Hvi->clone();
100 
101  u->zero();
102 
103  for( uint j=0; j<nvec2; ++j ) {
104  uint k = j+nvec1*i;
105  (*blocks_)[k]->apply(*u,*v_part.get(j),tol);
106  Hvi->plus(*u);
107  }
108  }
109  }
110 
111 }; // class BlockOperator
112 
113 
114 // Simplified Construction for a 2-by-2 block operator
115 template<class Real>
116 class BlockOperator2 : public LinearOperator<Real> {
117 
118  typedef Vector<Real> V; // ROL Vector
119  typedef LinearOperator<Real> OP; // Linear Operator
120 
121 private:
122 
123  Teuchos::RCP<OP> bkop_;
124 
125 public:
126  BlockOperator2( Teuchos::RCP<OP> &a11, Teuchos::RCP<OP> &a21,
127  Teuchos::RCP<OP> &a12, Teuchos::RCP<OP> &a22 ) {
128 
129  using std::vector;
130  using Teuchos::RCP;
131  using Teuchos::rcp;
132 
133  RCP<vector<RCP<OP> > > ops = rcp( new vector<RCP<OP> > );
134  ops->push_back(a11);
135  ops->push_back(a21);
136  ops->push_back(a12);
137  ops->push_back(a22);
138 
139  bkop_ = rcp( new BlockOperator<Real>(ops) );
140  }
141 
142  void apply( V &Hv, const V &v, Real &tol ) const {
143  bkop_->apply(Hv,v,tol);
144  }
145 };
146 
147 } // namespace ROL
148 
149 #endif // ROL_BLOCKOPERATOR_H
Defines the linear algebra of vector space on a generic partitioned vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
LinearOperator< Real > OP
Teuchos::RCP< const Vector< Real > > get(size_type i) const
PartitionedVector< Real > PV
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
virtual void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
BlockOperator(const Teuchos::RCP< OpVec > &blocks)
Provides the interface to apply a linear operator.
Teuchos::RCP< OpVec > blocks_
Provides the interface to apply a block operator to a partitioned vector.
LinearOperator< Real > OP
OpVec::size_type uint
BlockOperator2(Teuchos::RCP< OP > &a11, Teuchos::RCP< OP > &a21, Teuchos::RCP< OP > &a12, Teuchos::RCP< OP > &a22)
Teuchos::RCP< OP > bkop_
std::vector< Teuchos::RCP< OP > > OpVec