ConstrainedOptPack_QPSchurInitKKTSystemHessianFixedFree.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "ConstrainedOptPack_QPSchurInitKKTSystemHessianFixedFree.hpp"
00030 #include "ConstrainedOptPack_initialize_Q_R_Q_X.hpp"
00031 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
00032 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00033 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00034 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00035 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00036 #include "Midynamic_cast_verbose.h"
00037 #include "MiWorkspacePack.h"
00038 #include "Miprofile_hack.h"
00039 
00040 namespace LinAlgOpPack {
00041     using AbstractLinAlgPack::Vp_StMtV;
00042 }
00043 
00044 namespace ConstrainedOptPack {
00045 
00046 void QPSchurInitKKTSystemHessianFixedFree::initialize_kkt_system(
00047   const DVectorSlice&    g
00048   ,const MatrixOp&  G
00049   ,value_type           etaL
00050   ,const SpVectorSlice& dL
00051   ,const SpVectorSlice& dU
00052   ,const MatrixOp*  F
00053   ,BLAS_Cpp::Transp     trans_F
00054   ,const DVectorSlice*   f
00055   ,const DVectorSlice&   d
00056   ,const SpVectorSlice& nu
00057   ,size_type*           n_R
00058   ,i_x_free_t*          i_x_free
00059   ,i_x_fixed_t*         i_x_fixed
00060   ,bnd_fixed_t*         bnd_fixed
00061   ,j_f_decomp_t*        j_f_decomp
00062   ,DVector*              b_X
00063   ,Ko_ptr_t*            Ko
00064   ,DVector*              fo
00065   ) const
00066 {
00067   using Teuchos::dyn_cast;
00068   using LinAlgOpPack::V_mV;
00069   namespace rcp = MemMngPack;
00070   using Teuchos::Workspace;
00071   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00072 
00073 #ifdef PROFILE_HACK_ENABLED
00074   ProfileHackPack::ProfileTiming profile_timing( "QPSchurInitKKTSystemHessianFixedFree::initialize_kkt_system(...)" );
00075 #endif
00076 
00077   // Validate type of and convert G
00078 #ifdef _WINDOWS
00079   const MatrixSymOp&
00080     G_sym = dynamic_cast<const MatrixSymOp&>(G);
00081 #else
00082   const MatrixSymOp&
00083     G_sym = dyn_cast<const MatrixSymOp>(G);
00084 #endif
00085 
00086   const size_type nd = g.size();
00087 
00088   // Determine the number of initially fixed variables
00089   Workspace<EBounds> x_frfx(wss,nd);
00090   std::fill_n( &x_frfx[0], nd, FREE ); // make all free initially
00091   size_type
00092     num_init_fixed = 0;
00093   {
00094     const value_type inf_bnd = std::numeric_limits<value_type>::max();
00095     AbstractLinAlgPack::sparse_bounds_itr
00096       dLU_itr(
00097         dL.begin(), dL.end(), dL.offset(),
00098         dU.begin(), dU.end(), dU.offset(), inf_bnd );
00099     SpVectorSlice::const_iterator
00100       nu_itr = nu.begin(),
00101       nu_end = nu.end();
00102     const SpVector::difference_type o = nu.offset();
00103     while( !dLU_itr.at_end() || nu_itr != nu_end ) {
00104       if( dLU_itr.at_end() ) { // Add the rest of the elements in nu
00105         for( ; nu_itr != nu_end; ++num_init_fixed, ++nu_itr )
00106           x_frfx[nu_itr->indice() + o - 1] = ( nu_itr->value() > 0.0 ? UPPER : LOWER );
00107       }
00108       else { // Be carefull to add fixed dL(i) == dU(i)
00109         // Add elements in nu up to the current dLU_itr.indice()
00110         for( ; nu_itr != nu_end && nu_itr->indice() + o < dLU_itr.indice(); ++num_init_fixed, ++nu_itr )
00111           x_frfx[nu_itr->indice() + o - 1] = ( nu_itr->value() > 0.0 ? UPPER : LOWER );
00112         if( dLU_itr.lbound() == dLU_itr.ubound() ) {
00113           // This is a permanently fixed variable!
00114           x_frfx[dLU_itr.indice() - 1] = EQUALITY;
00115           ++num_init_fixed;
00116           // Don't add a duplicate entry in nu
00117           if( nu_itr != nu_end && nu_itr->indice() + o == dLU_itr.indice() )
00118             ++nu_itr;
00119         }
00120         ++dLU_itr;
00121       }
00122     }
00123   }
00124   assert( nd >= num_init_fixed );
00125 
00126   // n_R
00127   *n_R = nd - num_init_fixed;
00128   
00129   // Set up i_x_free[], i_x_fixed[], bnd_fixed[], and b_X
00130   i_x_free->resize(*n_R);
00131   i_x_fixed->resize(num_init_fixed+1);
00132   bnd_fixed->resize(num_init_fixed+1);
00133   b_X->resize(num_init_fixed+1);
00134   {
00135     const value_type inf_bnd = std::numeric_limits<value_type>::max();
00136     AbstractLinAlgPack::sparse_bounds_itr
00137       dLU_itr(
00138         dL.begin(), dL.end(), dL.offset(),
00139         dU.begin(), dU.end(), dU.offset(), inf_bnd );
00140     size_type i_R = 0, i_X = 0;
00141     for( size_type i = 1; i <= nd; ++i ) {
00142       const EBounds
00143         bnd_i = x_frfx[i-1];
00144       if( bnd_i == FREE ) {
00145         (*i_x_free)[i_R] = i;
00146         ++i_R;
00147       }
00148       else {
00149         (*i_x_fixed)[i_X] = i;
00150         (*bnd_fixed)[i_X] = bnd_i;
00151         assert( !dLU_itr.at_end() );    // find entry in b_X
00152         while( dLU_itr.indice() < i )
00153           ++dLU_itr;
00154         assert( dLU_itr.indice() == i );
00155         value_type b_X_val = 0.0;
00156         switch( bnd_i ) {
00157           case EQUALITY:
00158           case LOWER:
00159             b_X_val = dLU_itr.lbound();
00160             break;
00161           case UPPER:
00162             b_X_val = dLU_itr.ubound();
00163             break;
00164           default:
00165             TEST_FOR_EXCEPT(true); // Local error only?
00166         }
00167         (*b_X)[i_X] = b_X_val;
00168         ++i_X;
00169       }
00170     }
00171     (*i_x_fixed)[i_X] = nd+1;   // built-in relaxation variable
00172     (*bnd_fixed)[i_X] = LOWER;
00173     (*b_X)[i_X]       = etaL;
00174     ++i_X;
00175   }
00176   
00177   // j_f_decomp[] = empty
00178   j_f_decomp->resize(0);
00179 
00180   // Initialize temporary Q_R and Q_X (not including extra relaxation variable)
00181   Workspace<size_type>
00182     Q_R_row_i(wss,*n_R),
00183     Q_R_col_j(wss,*n_R),
00184     Q_X_row_i(wss,num_init_fixed),
00185     Q_X_col_j(wss,num_init_fixed);
00186   GenPermMatrixSlice
00187     Q_R, Q_X;
00188   initialize_Q_R_Q_X(
00189     *n_R,num_init_fixed,&(*i_x_free)[0],&(*i_x_fixed)[0],false
00190     ,&Q_R_row_i[0],&Q_R_col_j[0],&Q_R
00191     ,&Q_X_row_i[0],&Q_X_col_j[0],&Q_X
00192     );
00193 
00194   //
00195   // Create and initialize object for Ko = G_RR = Q_R'*G*Q_R
00196   //
00197 
00198   // Compute the dense matrix G_RR
00199   DMatrix G_RR_dense(*n_R,*n_R);
00200   DMatrixSliceSym sym_G_RR_dense(G_RR_dense(),BLAS_Cpp::lower);
00201   AbstractLinAlgPack::Mp_StPtMtP(
00202     &sym_G_RR_dense, 1.0, MatrixSymOp::DUMMY_ARG
00203     ,G_sym, Q_R, BLAS_Cpp::no_trans, 0.0 );
00204   // Initialize a factorization object for this matrix
00205   typedef Teuchos::RefCountPtr<MatrixSymPosDefCholFactor> G_RR_ptr_t;
00206   G_RR_ptr_t
00207     G_RR_ptr = new MatrixSymPosDefCholFactor();
00208   G_RR_ptr->initialize(sym_G_RR_dense);
00209   
00210   *Ko = Teuchos::rcp_implicit_cast<Ko_ptr_t::element_type>(G_RR_ptr); // Ko is initialized!
00211 
00212   // ToDo: (2001/07/05) We could be more carefull about how memory is initialized and reused
00213   // in the future but this implementation is just easier.
00214 
00215   // fo = - Q_R'*g - Q_R'*G*(Q_X*b_X)
00216   LinAlgOpPack::V_StMtV( fo, -1.0, Q_R, BLAS_Cpp::trans, g );
00217   if( num_init_fixed ) {
00218     SpVector b_XX;
00219     AbstractLinAlgPack::V_MtV( &b_XX, Q_X, BLAS_Cpp::no_trans, (*b_X)(1,num_init_fixed) );
00220     AbstractLinAlgPack::Vp_StPtMtV( &(*fo)(), -1.0, Q_R, BLAS_Cpp::trans, G, BLAS_Cpp::no_trans, b_XX() );
00221   }
00222 
00223 }
00224 
00225 } // end namesapce ConstrainedOptPack

Generated on Thu Sep 18 12:35:15 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1