ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_QPSchurInitKKTSystemHessianFixedFree.cpp
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include "ConstrainedOptPack_QPSchurInitKKTSystemHessianFixedFree.hpp"
00045 #include "ConstrainedOptPack_initialize_Q_R_Q_X.hpp"
00046 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
00047 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00048 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00049 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00050 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00051 #include "Midynamic_cast_verbose.h"
00052 #include "MiWorkspacePack.h"
00053 #include "Miprofile_hack.h"
00054 
00055 namespace LinAlgOpPack {
00056     using AbstractLinAlgPack::Vp_StMtV;
00057 }
00058 
00059 namespace ConstrainedOptPack {
00060 
00061 void QPSchurInitKKTSystemHessianFixedFree::initialize_kkt_system(
00062   const DVectorSlice&    g
00063   ,const MatrixOp&  G
00064   ,value_type           etaL
00065   ,const SpVectorSlice& dL
00066   ,const SpVectorSlice& dU
00067   ,const MatrixOp*  F
00068   ,BLAS_Cpp::Transp     trans_F
00069   ,const DVectorSlice*   f
00070   ,const DVectorSlice&   d
00071   ,const SpVectorSlice& nu
00072   ,size_type*           n_R
00073   ,i_x_free_t*          i_x_free
00074   ,i_x_fixed_t*         i_x_fixed
00075   ,bnd_fixed_t*         bnd_fixed
00076   ,j_f_decomp_t*        j_f_decomp
00077   ,DVector*              b_X
00078   ,Ko_ptr_t*            Ko
00079   ,DVector*              fo
00080   ) const
00081 {
00082   using Teuchos::dyn_cast;
00083   using LinAlgOpPack::V_mV;
00084   namespace rcp = MemMngPack;
00085   using Teuchos::Workspace;
00086   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00087 
00088 #ifdef PROFILE_HACK_ENABLED
00089   ProfileHackPack::ProfileTiming profile_timing( "QPSchurInitKKTSystemHessianFixedFree::initialize_kkt_system(...)" );
00090 #endif
00091 
00092   // Validate type of and convert G
00093 #ifdef _WINDOWS
00094   const MatrixSymOp&
00095     G_sym = dynamic_cast<const MatrixSymOp&>(G);
00096 #else
00097   const MatrixSymOp&
00098     G_sym = dyn_cast<const MatrixSymOp>(G);
00099 #endif
00100 
00101   const size_type nd = g.size();
00102 
00103   // Determine the number of initially fixed variables
00104   Workspace<EBounds> x_frfx(wss,nd);
00105   std::fill_n( &x_frfx[0], nd, FREE ); // make all free initially
00106   size_type
00107     num_init_fixed = 0;
00108   {
00109     const value_type inf_bnd = std::numeric_limits<value_type>::max();
00110     AbstractLinAlgPack::sparse_bounds_itr
00111       dLU_itr(
00112         dL.begin(), dL.end(), dL.offset(),
00113         dU.begin(), dU.end(), dU.offset(), inf_bnd );
00114     SpVectorSlice::const_iterator
00115       nu_itr = nu.begin(),
00116       nu_end = nu.end();
00117     const SpVector::difference_type o = nu.offset();
00118     while( !dLU_itr.at_end() || nu_itr != nu_end ) {
00119       if( dLU_itr.at_end() ) { // Add the rest of the elements in nu
00120         for( ; nu_itr != nu_end; ++num_init_fixed, ++nu_itr )
00121           x_frfx[nu_itr->indice() + o - 1] = ( nu_itr->value() > 0.0 ? UPPER : LOWER );
00122       }
00123       else { // Be carefull to add fixed dL(i) == dU(i)
00124         // Add elements in nu up to the current dLU_itr.indice()
00125         for( ; nu_itr != nu_end && nu_itr->indice() + o < dLU_itr.indice(); ++num_init_fixed, ++nu_itr )
00126           x_frfx[nu_itr->indice() + o - 1] = ( nu_itr->value() > 0.0 ? UPPER : LOWER );
00127         if( dLU_itr.lbound() == dLU_itr.ubound() ) {
00128           // This is a permanently fixed variable!
00129           x_frfx[dLU_itr.indice() - 1] = EQUALITY;
00130           ++num_init_fixed;
00131           // Don't add a duplicate entry in nu
00132           if( nu_itr != nu_end && nu_itr->indice() + o == dLU_itr.indice() )
00133             ++nu_itr;
00134         }
00135         ++dLU_itr;
00136       }
00137     }
00138   }
00139   TEUCHOS_TEST_FOR_EXCEPT( !(  nd >= num_init_fixed  ) );
00140 
00141   // n_R
00142   *n_R = nd - num_init_fixed;
00143   
00144   // Set up i_x_free[], i_x_fixed[], bnd_fixed[], and b_X
00145   i_x_free->resize(*n_R);
00146   i_x_fixed->resize(num_init_fixed+1);
00147   bnd_fixed->resize(num_init_fixed+1);
00148   b_X->resize(num_init_fixed+1);
00149   {
00150     const value_type inf_bnd = std::numeric_limits<value_type>::max();
00151     AbstractLinAlgPack::sparse_bounds_itr
00152       dLU_itr(
00153         dL.begin(), dL.end(), dL.offset(),
00154         dU.begin(), dU.end(), dU.offset(), inf_bnd );
00155     size_type i_R = 0, i_X = 0;
00156     for( size_type i = 1; i <= nd; ++i ) {
00157       const EBounds
00158         bnd_i = x_frfx[i-1];
00159       if( bnd_i == FREE ) {
00160         (*i_x_free)[i_R] = i;
00161         ++i_R;
00162       }
00163       else {
00164         (*i_x_fixed)[i_X] = i;
00165         (*bnd_fixed)[i_X] = bnd_i;
00166         TEUCHOS_TEST_FOR_EXCEPT( !(  !dLU_itr.at_end()  ) );    // find entry in b_X
00167         while( dLU_itr.indice() < i )
00168           ++dLU_itr;
00169         TEUCHOS_TEST_FOR_EXCEPT( !(  dLU_itr.indice() == i  ) );
00170         value_type b_X_val = 0.0;
00171         switch( bnd_i ) {
00172           case EQUALITY:
00173           case LOWER:
00174             b_X_val = dLU_itr.lbound();
00175             break;
00176           case UPPER:
00177             b_X_val = dLU_itr.ubound();
00178             break;
00179           default:
00180             TEUCHOS_TEST_FOR_EXCEPT(true); // Local error only?
00181         }
00182         (*b_X)[i_X] = b_X_val;
00183         ++i_X;
00184       }
00185     }
00186     (*i_x_fixed)[i_X] = nd+1;   // built-in relaxation variable
00187     (*bnd_fixed)[i_X] = LOWER;
00188     (*b_X)[i_X]       = etaL;
00189     ++i_X;
00190   }
00191   
00192   // j_f_decomp[] = empty
00193   j_f_decomp->resize(0);
00194 
00195   // Initialize temporary Q_R and Q_X (not including extra relaxation variable)
00196   Workspace<size_type>
00197     Q_R_row_i(wss,*n_R),
00198     Q_R_col_j(wss,*n_R),
00199     Q_X_row_i(wss,num_init_fixed),
00200     Q_X_col_j(wss,num_init_fixed);
00201   GenPermMatrixSlice
00202     Q_R, Q_X;
00203   initialize_Q_R_Q_X(
00204     *n_R,num_init_fixed,&(*i_x_free)[0],&(*i_x_fixed)[0],false
00205     ,&Q_R_row_i[0],&Q_R_col_j[0],&Q_R
00206     ,&Q_X_row_i[0],&Q_X_col_j[0],&Q_X
00207     );
00208 
00209   //
00210   // Create and initialize object for Ko = G_RR = Q_R'*G*Q_R
00211   //
00212 
00213   // Compute the dense matrix G_RR
00214   DMatrix G_RR_dense(*n_R,*n_R);
00215   DMatrixSliceSym sym_G_RR_dense(G_RR_dense(),BLAS_Cpp::lower);
00216   AbstractLinAlgPack::Mp_StPtMtP(
00217     &sym_G_RR_dense, 1.0, MatrixSymOp::DUMMY_ARG
00218     ,G_sym, Q_R, BLAS_Cpp::no_trans, 0.0 );
00219   // Initialize a factorization object for this matrix
00220   typedef Teuchos::RCP<MatrixSymPosDefCholFactor> G_RR_ptr_t;
00221   G_RR_ptr_t
00222     G_RR_ptr = new MatrixSymPosDefCholFactor();
00223   G_RR_ptr->initialize(sym_G_RR_dense);
00224   
00225   *Ko = Teuchos::rcp_implicit_cast<Ko_ptr_t::element_type>(G_RR_ptr); // Ko is initialized!
00226 
00227   // ToDo: (2001/07/05) We could be more carefull about how memory is initialized and reused
00228   // in the future but this implementation is just easier.
00229 
00230   // fo = - Q_R'*g - Q_R'*G*(Q_X*b_X)
00231   LinAlgOpPack::V_StMtV( fo, -1.0, Q_R, BLAS_Cpp::trans, g );
00232   if( num_init_fixed ) {
00233     SpVector b_XX;
00234     AbstractLinAlgPack::V_MtV( &b_XX, Q_X, BLAS_Cpp::no_trans, (*b_X)(1,num_init_fixed) );
00235     AbstractLinAlgPack::Vp_StPtMtV( &(*fo)(), -1.0, Q_R, BLAS_Cpp::trans, G, BLAS_Cpp::no_trans, b_XX() );
00236   }
00237 
00238 }
00239 
00240 } // end namesapce ConstrainedOptPack
00241 
00242 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends