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_QPSchurInitKKTSystemHessianFixedFree.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include "ConstrainedOptPack_QPSchurInitKKTSystemHessianFixedFree.hpp"
45 #include "ConstrainedOptPack_initialize_Q_R_Q_X.hpp"
46 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
47 #include "AbstractLinAlgPack_sparse_bounds.hpp"
48 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
49 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
50 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
51 #include "Midynamic_cast_verbose.h"
52 #include "MiWorkspacePack.h"
53 #include "Miprofile_hack.h"
54 
55 namespace LinAlgOpPack {
57 }
58 
59 namespace ConstrainedOptPack {
60 
62  const DVectorSlice& g
63  ,const MatrixOp& G
64  ,value_type etaL
65  ,const SpVectorSlice& dL
66  ,const SpVectorSlice& dU
67  ,const MatrixOp* F
68  ,BLAS_Cpp::Transp trans_F
69  ,const DVectorSlice* f
70  ,const DVectorSlice& d
71  ,const SpVectorSlice& nu
72  ,size_type* n_R
73  ,i_x_free_t* i_x_free
74  ,i_x_fixed_t* i_x_fixed
75  ,bnd_fixed_t* bnd_fixed
76  ,j_f_decomp_t* j_f_decomp
77  ,DVector* b_X
78  ,Ko_ptr_t* Ko
79  ,DVector* fo
80  ) const
81 {
82  using Teuchos::dyn_cast;
83  using LinAlgOpPack::V_mV;
84  namespace rcp = MemMngPack;
85  using Teuchos::Workspace;
87 
88 #ifdef PROFILE_HACK_ENABLED
89  ProfileHackPack::ProfileTiming profile_timing( "QPSchurInitKKTSystemHessianFixedFree::initialize_kkt_system(...)" );
90 #endif
91 
92  // Validate type of and convert G
93 #ifdef _WINDOWS
94  const MatrixSymOp&
95  G_sym = dynamic_cast<const MatrixSymOp&>(G);
96 #else
97  const MatrixSymOp&
98  G_sym = dyn_cast<const MatrixSymOp>(G);
99 #endif
100 
101  const size_type nd = g.size();
102 
103  // Determine the number of initially fixed variables
104  Workspace<EBounds> x_frfx(wss,nd);
105  std::fill_n( &x_frfx[0], nd, FREE ); // make all free initially
106  size_type
107  num_init_fixed = 0;
108  {
109  const value_type inf_bnd = std::numeric_limits<value_type>::max();
111  dLU_itr(
112  dL.begin(), dL.end(), dL.offset(),
113  dU.begin(), dU.end(), dU.offset(), inf_bnd );
114  SpVectorSlice::const_iterator
115  nu_itr = nu.begin(),
116  nu_end = nu.end();
117  const SpVector::difference_type o = nu.offset();
118  while( !dLU_itr.at_end() || nu_itr != nu_end ) {
119  if( dLU_itr.at_end() ) { // Add the rest of the elements in nu
120  for( ; nu_itr != nu_end; ++num_init_fixed, ++nu_itr )
121  x_frfx[nu_itr->indice() + o - 1] = ( nu_itr->value() > 0.0 ? UPPER : LOWER );
122  }
123  else { // Be carefull to add fixed dL(i) == dU(i)
124  // Add elements in nu up to the current dLU_itr.indice()
125  for( ; nu_itr != nu_end && nu_itr->indice() + o < dLU_itr.indice(); ++num_init_fixed, ++nu_itr )
126  x_frfx[nu_itr->indice() + o - 1] = ( nu_itr->value() > 0.0 ? UPPER : LOWER );
127  if( dLU_itr.lbound() == dLU_itr.ubound() ) {
128  // This is a permanently fixed variable!
129  x_frfx[dLU_itr.indice() - 1] = EQUALITY;
130  ++num_init_fixed;
131  // Don't add a duplicate entry in nu
132  if( nu_itr != nu_end && nu_itr->indice() + o == dLU_itr.indice() )
133  ++nu_itr;
134  }
135  ++dLU_itr;
136  }
137  }
138  }
139  TEUCHOS_TEST_FOR_EXCEPT( !( nd >= num_init_fixed ) );
140 
141  // n_R
142  *n_R = nd - num_init_fixed;
143 
144  // Set up i_x_free[], i_x_fixed[], bnd_fixed[], and b_X
145  i_x_free->resize(*n_R);
146  i_x_fixed->resize(num_init_fixed+1);
147  bnd_fixed->resize(num_init_fixed+1);
148  b_X->resize(num_init_fixed+1);
149  {
150  const value_type inf_bnd = std::numeric_limits<value_type>::max();
152  dLU_itr(
153  dL.begin(), dL.end(), dL.offset(),
154  dU.begin(), dU.end(), dU.offset(), inf_bnd );
155  size_type i_R = 0, i_X = 0;
156  for( size_type i = 1; i <= nd; ++i ) {
157  const EBounds
158  bnd_i = x_frfx[i-1];
159  if( bnd_i == FREE ) {
160  (*i_x_free)[i_R] = i;
161  ++i_R;
162  }
163  else {
164  (*i_x_fixed)[i_X] = i;
165  (*bnd_fixed)[i_X] = bnd_i;
166  TEUCHOS_TEST_FOR_EXCEPT( !( !dLU_itr.at_end() ) ); // find entry in b_X
167  while( dLU_itr.indice() < i )
168  ++dLU_itr;
169  TEUCHOS_TEST_FOR_EXCEPT( !( dLU_itr.indice() == i ) );
170  value_type b_X_val = 0.0;
171  switch( bnd_i ) {
172  case EQUALITY:
173  case LOWER:
174  b_X_val = dLU_itr.lbound();
175  break;
176  case UPPER:
177  b_X_val = dLU_itr.ubound();
178  break;
179  default:
180  TEUCHOS_TEST_FOR_EXCEPT(true); // Local error only?
181  }
182  (*b_X)[i_X] = b_X_val;
183  ++i_X;
184  }
185  }
186  (*i_x_fixed)[i_X] = nd+1; // built-in relaxation variable
187  (*bnd_fixed)[i_X] = LOWER;
188  (*b_X)[i_X] = etaL;
189  ++i_X;
190  }
191 
192  // j_f_decomp[] = empty
193  j_f_decomp->resize(0);
194 
195  // Initialize temporary Q_R and Q_X (not including extra relaxation variable)
196  Workspace<size_type>
197  Q_R_row_i(wss,*n_R),
198  Q_R_col_j(wss,*n_R),
199  Q_X_row_i(wss,num_init_fixed),
200  Q_X_col_j(wss,num_init_fixed);
201  GenPermMatrixSlice
202  Q_R, Q_X;
203  initialize_Q_R_Q_X(
204  *n_R,num_init_fixed,&(*i_x_free)[0],&(*i_x_fixed)[0],false
205  ,&Q_R_row_i[0],&Q_R_col_j[0],&Q_R
206  ,&Q_X_row_i[0],&Q_X_col_j[0],&Q_X
207  );
208 
209  //
210  // Create and initialize object for Ko = G_RR = Q_R'*G*Q_R
211  //
212 
213  // Compute the dense matrix G_RR
214  DMatrix G_RR_dense(*n_R,*n_R);
215  DMatrixSliceSym sym_G_RR_dense(G_RR_dense(),BLAS_Cpp::lower);
217  &sym_G_RR_dense, 1.0, MatrixSymOp::DUMMY_ARG
218  ,G_sym, Q_R, BLAS_Cpp::no_trans, 0.0 );
219  // Initialize a factorization object for this matrix
220  typedef Teuchos::RCP<MatrixSymPosDefCholFactor> G_RR_ptr_t;
221  G_RR_ptr_t
222  G_RR_ptr = new MatrixSymPosDefCholFactor();
223  G_RR_ptr->initialize(sym_G_RR_dense);
224 
225  *Ko = Teuchos::rcp_implicit_cast<Ko_ptr_t::element_type>(G_RR_ptr); // Ko is initialized!
226 
227  // ToDo: (2001/07/05) We could be more carefull about how memory is initialized and reused
228  // in the future but this implementation is just easier.
229 
230  // fo = - Q_R'*g - Q_R'*G*(Q_X*b_X)
231  LinAlgOpPack::V_StMtV( fo, -1.0, Q_R, BLAS_Cpp::trans, g );
232  if( num_init_fixed ) {
233  SpVector b_XX;
234  AbstractLinAlgPack::V_MtV( &b_XX, Q_X, BLAS_Cpp::no_trans, (*b_X)(1,num_init_fixed) );
235  AbstractLinAlgPack::Vp_StPtMtV( &(*fo)(), -1.0, Q_R, BLAS_Cpp::trans, G, BLAS_Cpp::no_trans, b_XX() );
236  }
237 
238 }
239 
240 } // end namesapce ConstrainedOptPack
241 
242 #endif // 0
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
T_To & dyn_cast(T_From &from)
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
const MatrixSymOpNonsing element_type
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t size_type
void initialize_kkt_system(const DVectorSlice &g, const MatrixOp &G, value_type etaL, const SpVectorSlice &dL, const SpVectorSlice &dU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const DVectorSlice *f, const DVectorSlice &d, const SpVectorSlice &nu, size_type *n_R, i_x_free_t *i_x_free, i_x_fixed_t *i_x_fixed, bnd_fixed_t *bnd_fixed, j_f_decomp_t *j_f_decomp, DVector *b_X, Ko_ptr_t *Ko, DVector *fo) const
Initialize the KKT system where initially fixed variables are removed and no equality constraints are...
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
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()