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_DecompositionSystemVarReductPermStd.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_DecompositionSystemVarReductPermStd.hpp"
43 #include "ConstrainedOptPack_DecompositionSystemVarReductImp.hpp"
44 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
45 #include "AbstractLinAlgPack_BasisSystemPerm.hpp"
46 #include "AbstractLinAlgPack_PermutationOut.hpp"
47 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
48 #include "Teuchos_Assert.hpp"
49 
50 namespace ConstrainedOptPack {
51 
52 // Constructors / initializers
53 
54 DecompositionSystemVarReductPermStd::DecompositionSystemVarReductPermStd(
55  const decomp_sys_imp_ptr_t& decomp_sys_imp
56  ,const basis_sys_ptr_t& basis_sys
57  ,bool basis_selected
58  ,EExplicitImplicit D_imp
59  ,EExplicitImplicit Uz_imp
60  )
61 {
62  this->initialize(decomp_sys_imp,basis_sys,basis_selected,D_imp,Uz_imp);
63 }
64 
66  const decomp_sys_imp_ptr_t& decomp_sys_imp
67  ,const basis_sys_ptr_t& basis_sys
68  ,bool basis_selected
69  ,EExplicitImplicit D_imp
70  ,EExplicitImplicit Uz_imp
71  )
72 {
73  decomp_sys_imp_ = decomp_sys_imp;
74  basis_sys_ = basis_sys;
75  basis_selected_ = basis_selected;
76  this->D_imp(D_imp);
77  this->Uz_imp(Uz_imp);
78 }
79 
80 // Overridden from DecompositionSystem
81 
83 {
84  return decomp_sys_imp()->n();
85 }
86 
88 {
89  return decomp_sys_imp()->m();
90 }
91 
93 {
94  return decomp_sys_imp()->r();
95 }
96 
98 {
99  return decomp_sys_imp()->equ_decomp();
100 }
101 
103 {
104  return decomp_sys_imp()->equ_undecomp();
105 }
106 
107 const VectorSpace::space_ptr_t
109 {
110  return decomp_sys_imp()->space_range();
111 }
112 
113 const VectorSpace::space_ptr_t
115 {
116  return decomp_sys_imp()->space_null();
117 }
118 
121 {
122  return decomp_sys_imp()->factory_Z();
123 }
124 
127 {
128  return decomp_sys_imp()->factory_Y();
129 }
130 
133 {
135  if( factory_R.get() != NULL )
136  return factory_R;
137  // Else assume that R will just be the basis matrix (coordinate decomposition!)
138  return basis_sys_->factory_C();
139 }
140 
143 {
144  return decomp_sys_imp()->factory_Uz();
145 }
146 
149 {
150  return decomp_sys_imp()->factory_Uy();
151 }
152 
154  std::ostream *out
155  ,EOutputLevel olevel
156  ,ERunTests test_what
157  ,const MatrixOp &Gc
158  ,MatrixOp *Z
159  ,MatrixOp *Y
160  ,MatrixOpNonsing *R
161  ,MatrixOp *Uz
162  ,MatrixOp *Uy
163  ,EMatRelations mat_rel
164  ) const
165 {
166  assert_basis_selected();
167  decomp_sys_imp()->update_decomp(
168  out,olevel,test_what,Gc,Z,Y
169  ,R,Uz,Uy,mat_rel
170  );
171 }
172 
174  std::ostream& out, const std::string& L ) const
175 {
176  // ToDo: Print basis permutation stuff also?
177  decomp_sys_imp()->print_update_decomp(out,L);
178 }
179 
180 // Overridden from DecompositionSystemVarReduct
181 
183 {
184  return basis_sys_.get() ? basis_sys_->var_indep() : Range1D::Invalid;
185 }
186 
188 {
189  return basis_sys_.get() ? basis_sys_->var_dep() : Range1D::Invalid;
190 }
191 
192 // @name Overridden from DecompositionSystemVarReductPerm
193 
196 {
197  return basis_sys()->factory_P_var();
198 }
199 
202 {
203  return basis_sys()->factory_P_equ();
204 }
205 
207 {
208  return basis_selected_;
209 }
210 
212  std::ostream *out
213  ,EOutputLevel olevel
214  ,ERunTests test_what
215  ,const Permutation &P_var
216  ,const Range1D &var_dep
217  ,const Permutation *P_equ
218  ,const Range1D *equ_decomp
219  ,const MatrixOp &Gc
220  ,MatrixOp *Z
221  ,MatrixOp *Y
222  ,MatrixOpNonsing *R
223  ,MatrixOp *Uz
224  ,MatrixOp *Uy
225  ,EMatRelations mat_rel
226  )
227 {
228  // Forward these setting on to the implementation.
229  decomp_sys_imp_->D_imp( this->D_imp() );
230  decomp_sys_imp_->Uz_imp( this->Uz_imp() );
231  // Get smart pointers to the basis matrix and the direct sensistivity matrices
232  // and remove references to these matrix objects from the other decomposition
233  // matrices by uninitializing them.
236  decomp_sys_imp_->get_basis_matrices(
237  out, olevel, test_what
238  ,Z, Y, R, Uz, Uy
239  ,&C_ptr
240  ,&D_ptr // May return D_ptr.get() == NULL if not explicit chosen
241  );
242  // Tell the basis system object to set this basis
243  try {
244  basis_sys_->set_basis(
245  P_var, var_dep
246  ,P_equ, equ_decomp
247  ,Gc
248  ,C_ptr.get()
249  ,D_ptr.get() // May be NULL
250  ,this->Uz_imp() == MAT_IMP_EXPLICIT ? Uz : NULL
251  ,(mat_rel == MATRICES_INDEP_IMPS
252  ? BasisSystem::MATRICES_INDEP_IMPS : BasisSystem::MATRICES_ALLOW_DEP_IMPS )
253  ,out
254  );
255  }
256  catch( const BasisSystem::SingularBasis& except ) {
257  if(out && olevel >= PRINT_BASIC_INFO)
258  *out << "Passed in basis is singular, throwing SingularDecomposition: "
259  << except.what() << std::endl;
262  ,"DecompositionSystemVarReductPermStd::set_decomp(...): Passed in basis selection "
263  "gave a singular basis matrix! : " << except.what() );
264  }
265  // If we get here the passed in basis selection is nonsingular and the basis matrices
266  // are updated. Now give them back to the decomp_sys_imp object and update the rest
267  // of the decomposition matrices.
268  const size_type
269  m = Gc.cols(),
270  r = C_ptr->rows();
271  decomp_sys_imp_->set_basis_matrices(
272  out, olevel, test_what
273  ,C_ptr
274  ,D_ptr // D_ptr.get() may be NULL
275  ,r > m ? Uz : NULL
276  ,basis_sys_ // Always reset
277  );
278  C_ptr = Teuchos::null;
279  D_ptr = Teuchos::null;
280  decomp_sys_imp()->update_decomp(
281  out,olevel,test_what,Gc,Z,Y,R
282  ,r > m ? Uz : NULL
283  ,r > m ? Uy : NULL
284  ,mat_rel
285  );
286  // We have a basis!
287  basis_selected_ = true;
288 }
289 
291  std::ostream *out
292  ,EOutputLevel olevel
293  ,ERunTests test_what
294  ,const Vector *nu
295  ,MatrixOp *Gc
296  ,Permutation *P_var
297  ,Range1D *var_dep
298  ,Permutation *P_equ
299  ,Range1D *equ_decomp
300  ,MatrixOp *Z
301  ,MatrixOp *Y
302  ,MatrixOpNonsing *R
303  ,MatrixOp *Uz
304  ,MatrixOp *Uy
305  ,EMatRelations mat_rel
306  )
307 {
308  // Forward these setting on to the implementation.
309  decomp_sys_imp_->D_imp( this->D_imp() );
310  decomp_sys_imp_->Uz_imp( this->Uz_imp() );
311  // Get smart pointers to the basis matrix and the direct sensistivity matrices
312  // and remove references to these matrix objects from the other decomposition
313  // matrices by uninitializing them.
316  //const bool unintialized_basis = decomp_sys_imp_->basis_sys()->var_dep().size() == 0;
317  decomp_sys_imp_->get_basis_matrices(
318  out, olevel, test_what
319  ,Z, Y, R, Uz, Uy
320  ,&C_ptr
321  ,&D_ptr // May return D_ptr.get() == NULL if not explicit chosen
322  );
323  // Ask the basis system object to select a basis
324  basis_sys_->select_basis(
325  nu
326  ,Gc
327  ,P_var, var_dep
328  ,P_equ, equ_decomp
329  ,C_ptr.get()
330  ,D_ptr.get() // May be NULL
331  ,this->Uz_imp() == MAT_IMP_EXPLICIT ? Uz : NULL
332  ,(mat_rel == MATRICES_INDEP_IMPS
333  ? BasisSystem::MATRICES_INDEP_IMPS : BasisSystem::MATRICES_ALLOW_DEP_IMPS )
334  ,out
335  );
336 
337  if( out && (int)olevel >= (int)PRINT_BASIC_INFO ) {
338  const Range1D var_indep = basis_sys_->var_indep(), equ_undecomp = basis_sys_->equ_undecomp();
339  *out
340  << "\nSelected a new basis\n"
341  << "\nbs.var_dep() = ["<<var_dep->lbound()<<","<<var_dep->ubound()<<"]"
342  << "\nds.var_indep() = ["<<var_indep.lbound()<<","<<var_indep.ubound()<<"]"
343  << "\nds.equ_decomp() = ["<<equ_decomp->lbound()<<","<<equ_decomp->ubound()<<"]"
344  << "\nds.equ_undecomp() = ["<<equ_undecomp.lbound()<<","<<equ_undecomp.ubound()<<"]"
345  << std::endl;
346  }
347  if( out && (int)olevel >= (int)PRINT_VECTORS ) {
348  *out
349  << "\nP_var =\n" << *P_var
350  << "\nP_equ =\n" << *P_equ
351  ;
352  }
353  if( out && (int)olevel >= (int)PRINT_EVERY_THING ) {
354  *out
355  << "\nGc =\n" << *Gc;
356  }
357 
358  // If we get here a nonsinguar basis selection has been made and the basis matrices
359  // are updated. Now give them back to the decomp_sys_imp object and update the rest
360  // of the decomposition matrices.
361  const size_type
362  //n = Gc->rows(),
363  m = Gc->cols(),
364  r = C_ptr->rows();
365  decomp_sys_imp_->set_basis_matrices(
366  out, olevel, test_what
367  ,C_ptr
368  ,D_ptr // D_ptr.get() may be NULL
369  ,r > m ? Uz : NULL
370  ,basis_sys_ // Always reset
371  );
372  C_ptr = Teuchos::null;
373  D_ptr = Teuchos::null;
374  decomp_sys_imp()->update_decomp(
375  out,olevel,test_what,*Gc,Z,Y,R
376  ,r > m ? Uz : NULL
377  ,r > m ? Uy : NULL
378  ,mat_rel
379  );
380  // We have a basis!
381  basis_selected_ = true;
382 }
383 
384 // private
385 
386 void DecompositionSystemVarReductPermStd::assert_basis_selected() const
387 {
389  !basis_selected_, std::logic_error
390  ,"DecompositionSystemVarReductPermStd::assert_basis_selected(): Error, "
391  "the methods set_decomp() or select_decomp() must be called first!" );
392 }
393 
394 } // end namespace ConstrainedOptPack
void initialize(const decomp_sys_imp_ptr_t &decomp_sys_imp, const basis_sys_ptr_t &basis_sys, bool basis_selected=false, EExplicitImplicit D_imp=MAT_IMP_AUTO, EExplicitImplicit Uz_imp=MAT_IMP_AUTO)
Initialize given decomposition system and basis system objects.
EOutputLevel
Enumeration for the amount of output to create from update_decomp().
ERunTests
Enumeration for if to run internal tests or not.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void set_decomp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const Permutation &P_var, const Range1D &var_dep, const Permutation *P_equ, const Range1D *equ_decomp, const MatrixOp &Gc, MatrixOp *Z, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uz, MatrixOp *Uy, EMatRelations mat_rel)
T * get() const
void update_decomp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const MatrixOp &Gc, MatrixOp *Z, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uz, MatrixOp *Uy, EMatRelations mat_rel) const
void print_update_decomp(std::ostream &out, const std::string &leading_str) const
size_t size_type
void select_decomp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const Vector *nu, MatrixOp *Gc, Permutation *P_var, Range1D *var_dep, Permutation *P_equ, Range1D *equ_decomp, MatrixOp *Z, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uz, MatrixOp *Uy, EMatRelations mat_rel)