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_DecompositionSystemOrthogonal.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 <assert.h>
43 
44 #include <typeinfo>
45 
46 #include "ConstrainedOptPack_DecompositionSystemOrthogonal.hpp"
47 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
48 #include "ConstrainedOptPack_MatrixDecompRangeOrthog.hpp"
49 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
50 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
51 #include "AbstractLinAlgPack_MatrixComposite.hpp"
52 #include "AbstractLinAlgPack_MatrixOpSubView.hpp"
53 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
54 #include "Teuchos_AbstractFactoryStd.hpp"
55 #include "Teuchos_dyn_cast.hpp"
56 #include "Teuchos_Assert.hpp"
57 
58 namespace ConstrainedOptPack {
59 
61  const VectorSpace::space_ptr_t &space_x
62  ,const VectorSpace::space_ptr_t &space_c
63  ,const basis_sys_ptr_t &basis_sys
64  ,const basis_sys_tester_ptr_t &basis_sys_tester
65  ,EExplicitImplicit D_imp
66  ,EExplicitImplicit Uz_imp
67  )
69  space_x, space_c, basis_sys, basis_sys_tester
70  ,D_imp, Uz_imp )
71 {}
72 
73 // Overridden from DecompositionSystem
74 
77 {
78  namespace rcp = MemMngPack;
79  return Teuchos::rcp(
81  );
82 }
83 
86 {
87  namespace rcp = MemMngPack;
88  return Teuchos::rcp(
90  );
91 }
92 
95 {
97 }
98 
99 // Overridden from DecompositionSystemVarReductImp
100 
102 {
103  *D_imp_used = MAT_IMP_EXPLICIT;
104 }
105 
106 DecompositionSystem::mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t
108  std::ostream *out
109  ,EOutputLevel olevel
110  ,MatrixOp *Y
111  ,MatrixOpNonsing *R
112  ,MatrixOp *Uy
113  ) const
114 {
115  namespace rcp = MemMngPack;
116  using Teuchos::dyn_cast;
117  typedef DecompositionSystem::mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t
118  C_ptr_t;
119 
120  //
121  // Get pointers to concreate matrices
122  //
123 
125  *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y) : NULL;
127  *R_orth = R ? &dyn_cast<MatrixDecompRangeOrthog>(*R) : NULL;
128  MatrixOpSubView
129  *Uy_cpst = Uy ? &dyn_cast<MatrixOpSubView>(*Uy) : NULL;
130 
131  //
132  // Get the smart pointer to the basis matrix object C and the
133  // matrix S = I + D'*D
134  //
135 
136  C_ptr_t C_ptr = Teuchos::null;
137  if(R_orth) {
138  C_ptr = Teuchos::rcp_const_cast<MatrixOpNonsing>( R_orth->C_ptr() ); // This could be NULL!
139  S_ptr_ = Teuchos::rcp_const_cast<MatrixSymOpNonsing>( R_orth->S_ptr() ); // ""
140  }
141 
142  //
143  // Uninitialize all of the matrices to remove references to C, D etc.
144  //
145 
146  if(Y_orth)
147  Y_orth->set_uninitialized();
148  if(R_orth)
149  R_orth->set_uninitialized();
150  if(Uy_cpst)
151  Uy_cpst->initialize(Teuchos::null);
152 
153  //
154  // Return the owned? basis matrix object C
155  //
156 
157  return C_ptr;
158 
159 }
160 
162  std::ostream *out
163  ,EOutputLevel olevel
164  ,const mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t &C
165  ,const mat_fcty_ptr_t::element_type::obj_ptr_t &D
166  ,MatrixOp *Y
167  ,MatrixOpNonsing *R
168  ,MatrixOp *Uy
169  ,EMatRelations mat_rel
170  ) const
171 {
172  namespace rcp = MemMngPack;
173  using Teuchos::dyn_cast;
174  using LinAlgOpPack::syrk;
175  typedef DecompositionSystem::mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t
176  C_ptr_t;
177  typedef DecompositionSystem::mat_fcty_ptr_t::element_type::obj_ptr_t
178  D_ptr_t;
179 
180  const size_type
181  //n = this->n(),
182  r = this->r();
183  const Range1D
184  var_dep(1,r);
185 
186  //
187  // Get pointers to concreate matrices
188  //
189 
191  *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y) : NULL;
193  *R_orth = R ? &dyn_cast<MatrixDecompRangeOrthog>(*R) : NULL;
194 
195  //
196  // Initialize the matrices
197  //
198 
199  if(Y_orth) {
200  D_ptr_t D_ptr = D;
201  if(mat_rel == MATRICES_INDEP_IMPS) {
202  D_ptr = D->clone();
204  D_ptr.get() == NULL, std::logic_error
205  ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
206  "The matrix class used for the direct sensitivity matrix D = inv(C)*N of type \'"
207  << typeName(*D) << "\' must return return.get() != NULL from the clone() method "
208  "since mat_rel == MATRICES_INDEP_IMPS!" );
209  }
210  Y_orth->initialize(
211  space_x() // space_cols
212  ,space_x()->sub_space(var_dep)->clone() // space_rows
213  ,MatrixIdentConcatStd::BOTTOM // top_or_bottom
214  ,-1.0 // alpha
215  ,D_ptr // D_ptr
216  ,BLAS_Cpp::trans // D_trans
217  );
218  }
219  if(R_orth) {
220  C_ptr_t C_ptr = C;
221  if(mat_rel == MATRICES_INDEP_IMPS) {
222  C_ptr = C->clone_mwons();
224  C_ptr.get() == NULL, std::logic_error
225  ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
226  "The matrix class used for the basis matrix C of type \'"
227  << typeName(*C) << "\' must return return.get() != NULL from the clone_mwons() method "
228  "since mat_rel == MATRICES_INDEP_IMPS!" );
229  }
230  D_ptr_t D_ptr = D;
231  if(mat_rel == MATRICES_INDEP_IMPS) {
232  D_ptr = D->clone();
234  D_ptr.get() == NULL, std::logic_error
235  ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
236  "The matrix class used for the direct sensitivity matrix D = inv(C)*N of type \'"
237  << typeName(*D) << "\' must return return.get() != NULL from the clone() method "
238  "since mat_rel == MATRICES_INDEP_IMPS!" );
239  }
240  if(S_ptr_.get() == NULL) {
241  S_ptr_ = basis_sys()->factory_S()->create();
242  }
243  try {
244  // S = I + (D)'*(D')'
245  dyn_cast<MatrixSymInitDiag>(*S_ptr_).init_identity(D_ptr->space_rows());
246  syrk(*D_ptr,BLAS_Cpp::trans,1.0,1.0,S_ptr_.get());
247  }
248  catch( const MatrixNonsing::SingularMatrix& except ) {
251  ,"DecompositionSystemOrthogonal::initialize_matrices(...) : Error, update of S failed : "
252  << except.what() );
253  }
254  R_orth->initialize(C_ptr,D_ptr,S_ptr_);
255  }
256  // ToDo: Implement for undecomposed equalities and general inequalities
257 }
258 
260  std::ostream& out, const std::string& L ) const
261 {
262  out
263  << L << "*** Orthogonal decompositon Y, R and Uy matrices (class DecompositionSystemOrthogonal)\n"
264  << L << "Y = [ I; -D' ] (using class MatrixIdentConcatStd)\n"
265  << L << "R = C*(I + D*D')\n"
266  << L << "Uy = E - F*D'\n"
267  ;
268 }
269 
270 } // end namespace ConstrainedOptPack
void print_update_matrices(std::ostream &out, const std::string &leading_str) const
EOutputLevel
Enumeration for the amount of output to create from update_decomp().
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T_To & dyn_cast(T_From &from)
const basis_sys_ptr_t & basis_sys() const
T * get() 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)
DecompositionSystemOrthogonal(const VectorSpace::space_ptr_t &space_x=Teuchos::null, const VectorSpace::space_ptr_t &space_c=Teuchos::null, const basis_sys_ptr_t &basis_sys=Teuchos::null, const basis_sys_tester_ptr_t &basis_sys_tester=Teuchos::null, EExplicitImplicit D_imp=MAT_IMP_EXPLICIT, EExplicitImplicit Uz_imp=MAT_IMP_EXPLICIT)
mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t uninitialize_matrices(std::ostream *out, EOutputLevel olevel, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uy) const
Concrete implementation class for a matrix vertically concatonated with an identity matrix...
size_t size_type
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
Matrix subclass for variable reduction orthogonal matrix R = Gc(:,con_decomp)'*Y. ...
virtual void initialize(const VectorSpace::space_ptr_t &space_cols, const VectorSpace::space_ptr_t &space_rows, ETopBottom top_or_bottom, value_type alpha, const D_ptr_t &D_ptr, BLAS_Cpp::Transp D_trans)
Setup with a matrix object.
size_type r() const
Returns this->basis_sys()->equ_decomp().size().
void initialize_matrices(std::ostream *out, EOutputLevel olevel, const mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t &C, const mat_fcty_ptr_t::element_type::obj_ptr_t &D, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uy, EMatRelations mat_rel) const
const VectorSpace::space_ptr_t & space_x() const
Specialization node implementation subclass of DecompositionSystem for variable reduction decompositi...
virtual void set_uninitialized()
Set the matrix to uninitialized.
std::string typeName(const T &t)