ConstrainedOptPack_QPSchurInitKKTSystemHessianFixedFree.cpp

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

Generated on Wed May 12 21:52:29 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7