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_QPSchurInitKKTSystemHessianSuperBasic.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_QPSchurInitKKTSystemHessianSuperBasic.hpp"
45 #include "ConstrainedOptPack_MatrixHessianSuperBasic.hpp"
46 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
47 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
48 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
49 #include "Midynamic_cast_verbose.h"
50 
51 namespace LinAlgOpPack {
53 }
54 
55 namespace ConstrainedOptPack {
56 
58  const DVectorSlice& g
59  ,const MatrixOp& G
60  ,value_type etaL
61  ,const SpVectorSlice& dL
62  ,const SpVectorSlice& dU
63  ,const MatrixOp* F
64  ,BLAS_Cpp::Transp trans_F
65  ,const DVectorSlice* f
66  ,const DVectorSlice& d
67  ,const SpVectorSlice& nu
68  ,size_type* n_R
69  ,i_x_free_t* i_x_free
70  ,i_x_fixed_t* i_x_fixed
71  ,bnd_fixed_t* bnd_fixed
72  ,j_f_decomp_t* j_f_decomp
73  ,DVector* b_X
74  ,Ko_ptr_t* Ko
75  ,DVector* fo
76  ) const
77 {
78  using BLAS_Cpp::trans;
79  using Teuchos::dyn_cast;
80  using LinAlgOpPack::V_mV;
81  using LinAlgOpPack::V_StMtV;
83 
84  // Validate type of and convert G
85  const MatrixHessianSuperBasic
86  *G_super_ptr = dynamic_cast<const MatrixHessianSuperBasic*>(&G);
87 
88  if( G_super_ptr == NULL ) {
89  init_kkt_full_.initialize_kkt_system(
90  g,G,etaL,dL,dU,F,trans_F,f,d,nu,n_R,i_x_free,i_x_fixed,bnd_fixed
91  ,j_f_decomp,b_X,Ko,fo);
92  return;
93  }
94 
95  const MatrixHessianSuperBasic
96  &G_super = *G_super_ptr;
97 
98  // get some stuff
99  const GenPermMatrixSlice
100  &Q_R = G_super.Q_R(),
101  &Q_X = G_super.Q_X();
102  const size_type
103  nd = G_super.rows(),
104  nd_R = Q_R.cols(),
105  nd_X = Q_X.cols();
106  TEUCHOS_TEST_FOR_EXCEPT( !( nd_R + nd_X == nd ) );
107 
108  // Setup output arguments
109 
110  // n_R = nd_R
111  *n_R = nd_R;
112  // i_x_free[(G.Q_R.begin()+l-1)->col_j()-1] = (G.Q_R.begin()+l-1)->row_i(), l = 1...nd_R
113  i_x_free->resize( Q_R.is_identity() ? 0: nd_R );
114  if( nd_R && !Q_R.is_identity() ) {
115  GenPermMatrixSlice::const_iterator
116  Q_itr = Q_R.begin();
117  for( ; Q_itr != Q_R.end(); ++Q_itr ) {
118  const size_type i = Q_itr->row_i();
119  const size_type k = Q_itr->col_j();
120  TEUCHOS_TEST_FOR_EXCEPT( !( 0 < i && i <= nd ) );
121  TEUCHOS_TEST_FOR_EXCEPT( !( 0 < k && k <= nd_R ) );
122  (*i_x_free)[k-1] = i;
123  }
124  }
125  // i_x_fixed[]
126  i_x_fixed->resize(nd_X+1);
127  if(nd_X) {
128  // i_x_fixed[(G.Q_X.begin()+l-1)->col_j()-1] = (G.Q_X.begin()+l-1)->row_i(), l = 1...nd_X
129  GenPermMatrixSlice::const_iterator
130  Q_itr = Q_X.begin();
131  for( ; Q_itr != Q_X.end(); ++Q_itr ) {
132  const size_type i = Q_itr->row_i();
133  const size_type k = Q_itr->col_j();
134  TEUCHOS_TEST_FOR_EXCEPT( !( 0 < i && i <= nd ) );
135  TEUCHOS_TEST_FOR_EXCEPT( !( 0 < k && k <= nd_X ) );
136  (*i_x_fixed)[k-1] = i;
137  }
138  }
139  (*i_x_fixed)[nd_X] = nd+1; // relaxation is always initially active
140  // bnd_fixed[]
141  bnd_fixed->resize(nd_X+1);
142  if(nd_X) {
143  // bnd_fixed[l-1] = G.bnd_fixed[l-1], l = 1...nd_X
144  typedef MatrixHessianSuperBasic MHSB;
145  const MHSB::bnd_fixed_t
146  &bnd_fixed_from = G_super.bnd_fixed();
147  TEUCHOS_TEST_FOR_EXCEPT( !( bnd_fixed_from.size() == nd_X ) );
148  MHSB::bnd_fixed_t::const_iterator
149  bnd_from_itr = bnd_fixed_from.begin();
150  bnd_fixed_t::iterator
151  bnd_to_itr = bnd_fixed->begin();
152  std::copy( bnd_from_itr, bnd_fixed_from.end(), bnd_to_itr );
153  }
154  (*bnd_fixed)[nd_X] = LOWER; // relaxation is always initially active
155  // j_f_decomp[]
156  j_f_decomp->resize(0);
157  // b_X
158  b_X->resize(nd_X+1);
159  if(nd_X) {
160  // b_X[l-1] = { dL(i) if bnd_fixed[l-1] == LOWER or EQUALITY
161  // dU(i) if bnd_fixed[l-1] == UPPER }
162  // , l = 1...nd_X
163  // (where i = i_x_fixed[l-1])
164  bnd_fixed_t::const_iterator
165  bnd_itr = const_cast<const bnd_fixed_t&>(*bnd_fixed).begin(),
166  bnd_itr_end = const_cast<const bnd_fixed_t&>(*bnd_fixed).begin() + nd_X;
167  i_x_fixed_t::const_iterator
168  i_x_itr = const_cast<const i_x_fixed_t&>(*i_x_fixed).begin();
169  DVector::iterator
170  b_X_itr = b_X->begin();
171  const SpVectorSlice::element_type
172  *ele = NULL;
173  for( ; bnd_itr != bnd_itr_end; ++bnd_itr, ++i_x_itr, ++b_X_itr ) {
174  const size_type i = *i_x_itr;
175  switch(*bnd_itr) {
176  case LOWER:
177  case EQUALITY:
178  *b_X_itr = (ele = dL.lookup_element(i))->value(); // Should not be null!
179  break;
180  case UPPER:
181  *b_X_itr = (ele = dU.lookup_element(i))->value(); // Should not be null!
182  break;
183  default:
185  }
186  }
187  }
188  (*b_X)[nd_X] = etaL; // relaxation is always initially active
189  // Ko = G.B_RR
190  *Ko = G_super.B_RR_ptr(); // now B_RR is a shared object
191  // fo = - G.Q_R'*g - op(G.B_RX)*b_X(1:nd_X)
192  V_StMtV( fo, -1.0, Q_R, trans, g );
193  if( nd_X && G_super.B_RX_ptr().get() )
194  Vp_StMtV( &(*fo)(), -1.0, *G_super.B_RX_ptr(), G_super.B_RX_trans(), (*b_X)(1,nd_X) );
195 
196 }
197 
198 } // end namesapce ConstrainedOptPack
199 
200 #endif // 0
T_To & dyn_cast(T_From &from)
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
size_t size_type
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 initialize_kkt_system(const Vector &g, const MatrixOp &G, value_type etaL, const Vector *dL, const Vector *dU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, const Vector *d, const Vector *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 all variables (except the relaxation variable) are initially free and...
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 the variables are initiallly fixed and free and no constraints are in...
Transp
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)