ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixDecompRangeOrthog.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 "ConstrainedOptPack_MatrixDecompRangeOrthog.hpp"
43 #include "AbstractLinAlgPack_VectorSpace.hpp"
44 #include "AbstractLinAlgPack_VectorStdOps.hpp"
45 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
46 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
47 #include "AbstractLinAlgPack_AssertOp.hpp"
48 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
49 #include "Teuchos_Assert.hpp"
50 #include "Teuchos_FancyOStream.hpp"
51 
52 namespace ConstrainedOptPack {
53 
54 // Constructors/initializers
55 
57 {}
58 
60  const C_ptr_t &C_ptr
61  ,const D_ptr_t &D_ptr
62  ,const S_ptr_t &S_ptr
63  )
64 {
65  this->initialize(C_ptr,D_ptr,S_ptr);
66 }
67 
69  const C_ptr_t &C_ptr
70  ,const D_ptr_t &D_ptr
71  ,const S_ptr_t &S_ptr
72  )
73 {
74  const char func_name[] = "MatrixDecompRangeOrthog::initialize(...)";
76  C_ptr.get() == NULL, std::invalid_argument
77  ,func_name << " : Error!" );
79  D_ptr.get() == NULL, std::invalid_argument
80  ,func_name << " : Error!" );
82  S_ptr.get() == NULL, std::invalid_argument
83  ,func_name << " : Error!" );
84 #ifdef ABSTRACT_LIN_ALG_PACK_CHECK_VEC_SPCS
85  bool is_compatible = C_ptr->space_rows().is_compatible(D_ptr->space_cols());
87  !is_compatible, VectorSpace::IncompatibleVectorSpaces
88  ,func_name << " : Error, C_ptr->space_rows().is_compatible(D_ptr->space_cols()) == false!" );
89  is_compatible = S_ptr->space_cols().is_compatible(D_ptr->space_rows());
91  !is_compatible, VectorSpace::IncompatibleVectorSpaces
92  ,func_name << " : Error, S_ptr->space_cols().is_compatible(D_ptr->space_rows()) == false!" );
93 #endif
94  C_ptr_ = C_ptr;
95  D_ptr_ = D_ptr;
96  S_ptr_ = S_ptr;
97 }
98 
100 {
101  namespace rcp = MemMngPack;
102  C_ptr_ = Teuchos::null;
103  D_ptr_ = Teuchos::null;
104  S_ptr_ = Teuchos::null;
105 }
106 
107 // Overridden from MatrixOp
108 
110 {
111  return C_ptr_.get() ? C_ptr_->rows() : 0;
112 }
113 
115 {
116  return C_ptr_.get() ? C_ptr_->cols() : 0;
117 }
118 
119 const VectorSpace& MatrixDecompRangeOrthog::space_cols() const
120 {
121  return C_ptr_->space_cols();
122 }
123 
124 const VectorSpace& MatrixDecompRangeOrthog::space_rows() const
125 {
126  return C_ptr_->space_rows();
127 }
128 
129 std::ostream& MatrixDecompRangeOrthog::output(std::ostream& out_arg) const
130 {
131  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
132  Teuchos::OSTab tab(out);
133  *out
134  << "This is a " << this->rows() << " x " << this->cols()
135  << " nonsingular matrix (I + D'*D)*C with inverse inv(C')*(I + D*inv(S)*D') where C, D and S are:\n";
136  tab.incrTab();
137  *out << "C =\n" << *C_ptr();
138  *out << "D =\n" << *D_ptr();
139  *out << "S =\n" << *S_ptr();
140  return out_arg;
141 }
142 
144  VectorMutable* y, value_type a, BLAS_Cpp::Transp R_trans
145  , const Vector& x, value_type b
146  ) const
147 {
148  using BLAS_Cpp::no_trans;
149  using BLAS_Cpp::trans;
152  using LinAlgOpPack::Vp_V;
153  using LinAlgOpPack::V_MtV;
154  using LinAlgOpPack::Vp_MtV;
155  using LinAlgOpPack::V_StMtV;
156 
157  assert_initialized("MatrixDecompRangeOrthog::Vp_StMtV(...)");
159 
160  const MatrixOpNonsing &C = *C_ptr_;
161  const MatrixOp &D = *D_ptr_;
162  const MatrixSymOpNonsing &S = *S_ptr_;
163  //
164  // y = b*y + a*op(R)*x
165  //
166  VectorSpace::vec_mut_ptr_t // ToDo: make workspace!
167  tI = D.space_rows().create_member(),
168  tD = D.space_cols().create_member();
169  if(R_trans == no_trans) {
170  //
171  // y = b*y + a*C*(I + D*D')*x
172  //
173  // =>
174  //
175  // tI = D'*x
176  // tD = x + D*tI
177  // y = b*y + a*C*tD
178  //
179  V_MtV( tI.get(), D, trans, x ); // tI = D'*x
180  *tD = x; // tD = x + D*tI
181  Vp_MtV( tD.get(), D, no_trans, *tI ); // ""
182  Vp_StMtV( y, a, C, no_trans, *tD, b ); // y = b*y + a*C*tD
183  }
184  else {
185  //
186  // y = b*y + a*(I + D*D')*C'*x
187  // = b*y + a*C'*x + D*D'*(a*C'*x)
188  //
189  // =>
190  //
191  // tD = a*C'*x
192  // tI = D'*tD
193  // y = b*y + D*tI
194  // y += tD
195  //
196  V_StMtV( tD.get(), a, C, trans, x ); // tD = a*C'*x
197  V_MtV( tI.get(), D, trans, *tD ); // tI = D'*tD
198  Vp_MtV( y, D, no_trans, *tI, b ); // y = b*y + D*tI
199  Vp_V( y, *tD ); // y += tD
200  }
201 }
202 
203 // Overridden from MatrixOpNonsing
204 
206  VectorMutable* y, BLAS_Cpp::Transp R_trans
207  , const Vector& x
208  ) const
209 {
210  using BLAS_Cpp::no_trans;
211  using BLAS_Cpp::trans;
215  using LinAlgOpPack::Vp_V;
216  using LinAlgOpPack::V_MtV;
217  using LinAlgOpPack::V_StMtV;
218 
219  assert_initialized("MatrixDecompRangeOrthog::V_InvMtV(...)");
221 
222  const MatrixOpNonsing &C = *C_ptr_;
223  const MatrixOp &D = *D_ptr_;
224  const MatrixSymOpNonsing &S = *S_ptr_;
225  //
226  // y = inv(op(R))*x
227  //
228  VectorSpace::vec_mut_ptr_t // ToDo: make workspace!
229  tIa = D.space_rows().create_member(),
230  tIb = D.space_rows().create_member(),
231  tD = D.space_cols().create_member();
232  if(R_trans == no_trans) {
233  //
234  // y = (I - D*inv(S)*D')*inv(C)*x
235  // = inv(C)*x - D*inv(S)*D'*inv(C)*x
236  //
237  // =>
238  //
239  // y = inv(C)*x
240  // tIa = D'*y
241  // tIb = inv(S)*tIa
242  // y += -D*tIb
243  //
244  V_InvMtV( y, C, no_trans, x ); // y = inv(C)*x
245  V_MtV( tIa.get(), D, trans, *y ); // tIa = D'*y
246  V_InvMtV( tIb.get(), S, no_trans, *tIa ); // tIb = inv(S)*tIa
247  Vp_StMtV( y, -1.0, D, no_trans, *tIb ); // y += -D*tIb
248  }
249  else {
250  //
251  // y = inv(C')*(I - D*inv(S)*D')*x
252  //
253  // =>
254  //
255  // tIa = D'*x
256  // tIb = inv(S)*tIa
257  // tD = x
258  // tD += -D*tIb
259  // y = Inv(C')*tD
260  //
261  V_MtV( tIa.get(), D, trans, x ); // tIa = D'*x
262  V_InvMtV( tIb.get(), S, no_trans, *tIa ); // tIb = inv(S)*tIa
263  *tD = x; // tD = x
264  Vp_StMtV( tD.get(), -1.0, D, no_trans, *tIb ); // tD += -D*tIb
265  V_InvMtV( y, C, trans, *tD ); // y = Inv(C')*tD
266  }
267 }
268 
269 // private
270 
271 void MatrixDecompRangeOrthog::assert_initialized(const char func_name[]) const
272 {
274  C_ptr_.get() == NULL, std::logic_error
275  ,func_name << " : Error, Must call initialize(...) first!" );
276 }
277 
278 } // end namespace ConstrainedOptPack
virtual const VectorSpace & space_rows() const =0
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
basic_OSTab< CharT, Traits > & incrTab(const int tabs=1)
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
void initialize(const C_ptr_t &C_ptr, const D_ptr_t &D_ptr, const S_ptr_t &S_ptr)
Initialize the matrix object.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
size_t size_type
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
Transp trans_not(Transp _trans)
void Vp_MtV_assert_compatibility(VectorMutable *v_lhs, const MatrixOp &m_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
virtual vec_mut_ptr_t create_member() const =0
Transp
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)