MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_EvalNewPointTailoredApproachOrthogonal_Step.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 "MoochoPack_EvalNewPointTailoredApproachOrthogonal_Step.hpp"
43 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
44 #include "NLPInterfacePack_NLPDirect.hpp"
45 #include "AbstractLinAlgPack_MatrixComposite.hpp"
46 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
47 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
48 #include "AbstractLinAlgPack_VectorSpace.hpp"
49 #include "AbstractLinAlgPack_VectorStdOps.hpp"
50 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
51 #include "AbstractLinAlgPack_AssertOp.hpp"
52 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
53 #include "Teuchos_dyn_cast.hpp"
54 #include "Teuchos_Assert.hpp"
55 
56 namespace MoochoPack {
57 
58 EvalNewPointTailoredApproachOrthogonal_Step::EvalNewPointTailoredApproachOrthogonal_Step(
59  const deriv_tester_ptr_t &deriv_tester
60  ,const bounds_tester_ptr_t &bounds_tester
61  ,EFDDerivTesting fd_deriv_testing
62  )
63  :EvalNewPointTailoredApproach_Step(deriv_tester,bounds_tester,fd_deriv_testing)
64 {}
65 
66 // protected
67 
69  MatrixOp *Y
70  ,MatrixOp *Uy
71  )
72 {
73  using Teuchos::dyn_cast;
74 
75  MatrixIdentConcatStd
76  *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y) : NULL;
77  MatrixComposite
78  *Uy_cpst = Uy ? &dyn_cast<MatrixComposite>(*Uy) : NULL;
79 
80  if(Y_orth)
81  Y_orth->set_uninitialized();
82  TEUCHOS_TEST_FOR_EXCEPT( !( Uy_cpst == NULL ) ); // ToDo: Implement for undecomposed equalities
83 }
84 
86  const NLPDirect &nlp
87  ,const D_ptr_t &D
88  ,VectorMutable *py
89  ,MatrixOp *Y
90  ,MatrixOp *Uy
91  ,EJournalOutputLevel olevel
92  ,std::ostream &out
93  )
94 {
95  namespace rcp = MemMngPack;
96  using Teuchos::dyn_cast;
97  using LinAlgOpPack::syrk;
98 
99  const size_type
100  n = nlp.n(),
101  r = nlp.r();
102  const Range1D
103  var_dep(1,r),
104  var_indep(r+1,n),
105  con_decomp = nlp.con_decomp(),
106  con_undecomp = nlp.con_undecomp();
107 
108  //
109  // Get pointers to concreate matrices
110  //
111 
112  MatrixIdentConcatStd
113  *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y) : NULL;
114  MatrixComposite
115  *Uy_cpst = Uy ? &dyn_cast<MatrixComposite>(*Uy) : NULL;
116 
117  //
118  // Initialize the matrices
119  //
120 
121  // Y
122  if(Y_orth) {
123  D_ptr_t D_ptr = D;
124 // if(mat_rel == MATRICES_INDEP_IMPS) {
125 // D_ptr = D->clone();
126 // TEUCHOS_TEST_FOR_EXCEPTION(
127 // D_ptr.get() == NULL, std::logic_error
128 // ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
129 // "The matrix class used for the direct sensitivity matrix D = inv(C)*N of type \'"
130 // << typeName(*D) << "\' must return return.get() != NULL from the clone() method "
131 // "since mat_rel == MATRICES_INDEP_IMPS!" );
132 // }
133  Y_orth->initialize(
134  nlp.space_x() // space_cols
135  ,nlp.space_x()->sub_space(var_dep)->clone() // space_rows
136  ,MatrixIdentConcatStd::BOTTOM // top_or_bottom
137  ,-1.0 // alpha
138  ,D_ptr // D_ptr
139  ,BLAS_Cpp::trans // D_trans
140  );
141  }
142 
143  // S
144  if(S_ptr_.get() == NULL) {
145  S_ptr_ = nlp.factory_S()->create();
146  }
147  // S = I + (D)'*(D')'
148  dyn_cast<MatrixSymInitDiag>(*S_ptr_).init_identity(D->space_rows());
149  syrk(*D,BLAS_Cpp::trans,1.0,1.0,S_ptr_.get());
150 
151  TEUCHOS_TEST_FOR_EXCEPT( !( Uy_cpst == NULL ) ); // ToDo: Implement for undecomposed equalities
152 
153  recalc_py(*D,py,olevel,out);
154 
155 }
156 
158  const MatrixOp &D
159  ,VectorMutable *py
160  ,EJournalOutputLevel olevel
161  ,std::ostream &out
162  )
163 {
164  using BLAS_Cpp::no_trans;
165  using BLAS_Cpp::trans;
168  using LinAlgOpPack::V_MtV;
169 
170  const MatrixSymOpNonsing &S = *S_ptr_;
171 
172  VectorSpace::vec_mut_ptr_t // ToDo: make workspace!
173  tIa = D.space_rows().create_member(),
174  tIb = D.space_rows().create_member();
175  //
176  // py = -inv(R)*c
177  // py = -((I - D*inv(S)*D')*inv(C))*c
178  // = -(I - D*inv(S)*D')*(-py)
179  // = py - D*inv(S)*D'*py
180  //
181  // =>
182  //
183  // tIa = D'*py
184  // tIb = inv(S)*tIa
185  // py += -D*tIb
186  //
187  V_MtV( tIa.get(), D, trans, *py ); // tIa = D'*py
188  V_InvMtV( tIb.get(), S, no_trans, *tIa ); // tIb = inv(S)*tIa
189  Vp_StMtV( py, -1.0, D, no_trans, *tIb ); // py += -D*tIb
190 }
191 
193  std::ostream& out, const std::string& L
194  ) const
195 {
196  out
197  << L << "*** Orthogonal decomposition\n"
198  << L << "py = inv(I + D*D') * py <: space_range\n"
199  << L << "Y = [ I ; -D' ] <: space_x|space_range\n"
200  << L << "Uy = ???\n"
201  ;
202 }
203 
204 } // end namespace MoochoPack
T_To & dyn_cast(T_From &from)
void print_calc_py_Y_Uy(std::ostream &out, const std::string &leading_str) const
T * get() const
void calc_py_Y_Uy(const NLPDirect &nlp, const D_ptr_t &D, VectorMutable *py, MatrixOp *Y, MatrixOp *Uy, EJournalOutputLevel olevel, std::ostream &out)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t size_type
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
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)
void recalc_py(const MatrixOp &D, VectorMutable *py, EJournalOutputLevel olevel, std::ostream &out)
Base class for evaluating a new point for the "Tailored Approach".
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)