ConstrainedOptPack_QPSchurInitKKTSystemHessianSuperBasic.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_QPSchurInitKKTSystemHessianSuperBasic.hpp"
00030 #include "ConstrainedOptPack_MatrixHessianSuperBasic.hpp"
00031 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00032 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00033 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00034 #include "Midynamic_cast_verbose.h"
00035 
00036 namespace LinAlgOpPack {
00037   using AbstractLinAlgPack::Vp_StMtV;
00038 }
00039 
00040 namespace ConstrainedOptPack {
00041 
00042 void QPSchurInitKKTSystemHessianSuperBasic::initialize_kkt_system(
00043   const DVectorSlice&    g
00044   ,const MatrixOp&  G
00045   ,value_type           etaL
00046   ,const SpVectorSlice& dL
00047   ,const SpVectorSlice& dU
00048   ,const MatrixOp*  F
00049   ,BLAS_Cpp::Transp     trans_F
00050   ,const DVectorSlice*   f
00051   ,const DVectorSlice&   d
00052   ,const SpVectorSlice& nu
00053   ,size_type*           n_R
00054   ,i_x_free_t*          i_x_free
00055   ,i_x_fixed_t*         i_x_fixed
00056   ,bnd_fixed_t*         bnd_fixed
00057   ,j_f_decomp_t*        j_f_decomp
00058   ,DVector*              b_X
00059   ,Ko_ptr_t*            Ko
00060   ,DVector*              fo
00061   ) const
00062 {
00063   using BLAS_Cpp::trans;
00064   using Teuchos::dyn_cast;
00065   using LinAlgOpPack::V_mV;
00066   using LinAlgOpPack::V_StMtV;
00067   using AbstractLinAlgPack::Vp_StMtV;
00068 
00069   // Validate type of and convert G
00070   const MatrixHessianSuperBasic
00071     *G_super_ptr = dynamic_cast<const MatrixHessianSuperBasic*>(&G);
00072 
00073   if( G_super_ptr == NULL ) {
00074     init_kkt_full_.initialize_kkt_system(
00075       g,G,etaL,dL,dU,F,trans_F,f,d,nu,n_R,i_x_free,i_x_fixed,bnd_fixed
00076       ,j_f_decomp,b_X,Ko,fo);
00077     return;
00078   }
00079 
00080   const MatrixHessianSuperBasic
00081     &G_super = *G_super_ptr;
00082 
00083   // get some stuff
00084   const GenPermMatrixSlice
00085     &Q_R = G_super.Q_R(),
00086     &Q_X = G_super.Q_X();
00087   const size_type
00088     nd   = G_super.rows(),
00089     nd_R = Q_R.cols(),
00090     nd_X = Q_X.cols();
00091   assert( nd_R + nd_X == nd );
00092 
00093   // Setup output arguments
00094 
00095   // n_R = nd_R
00096   *n_R = nd_R;
00097   // i_x_free[(G.Q_R.begin()+l-1)->col_j()-1] = (G.Q_R.begin()+l-1)->row_i(), l = 1...nd_R
00098   i_x_free->resize( Q_R.is_identity() ? 0: nd_R );
00099   if( nd_R && !Q_R.is_identity() ) {
00100     GenPermMatrixSlice::const_iterator
00101       Q_itr = Q_R.begin();
00102     for( ; Q_itr != Q_R.end(); ++Q_itr ) {
00103       const size_type i = Q_itr->row_i();
00104       const size_type k = Q_itr->col_j();
00105       assert( 0 < i && i <= nd );
00106       assert( 0 < k && k <= nd_R );
00107       (*i_x_free)[k-1] = i;
00108     }
00109   }
00110   // i_x_fixed[]
00111   i_x_fixed->resize(nd_X+1);
00112   if(nd_X) {
00113     // i_x_fixed[(G.Q_X.begin()+l-1)->col_j()-1] = (G.Q_X.begin()+l-1)->row_i(), l = 1...nd_X
00114     GenPermMatrixSlice::const_iterator
00115       Q_itr = Q_X.begin();
00116     for( ; Q_itr != Q_X.end(); ++Q_itr ) {
00117       const size_type i = Q_itr->row_i();
00118       const size_type k = Q_itr->col_j();
00119       assert( 0 < i && i <= nd );
00120       assert( 0 < k && k <= nd_X );
00121       (*i_x_fixed)[k-1] = i;
00122     }
00123   }
00124   (*i_x_fixed)[nd_X] = nd+1; // relaxation is always initially active
00125   // bnd_fixed[]
00126   bnd_fixed->resize(nd_X+1);
00127   if(nd_X) {
00128     // bnd_fixed[l-1] = G.bnd_fixed[l-1], l = 1...nd_X
00129     typedef MatrixHessianSuperBasic MHSB;
00130     const MHSB::bnd_fixed_t
00131       &bnd_fixed_from = G_super.bnd_fixed();
00132     assert(bnd_fixed_from.size() == nd_X);
00133     MHSB::bnd_fixed_t::const_iterator
00134       bnd_from_itr = bnd_fixed_from.begin();
00135     bnd_fixed_t::iterator
00136       bnd_to_itr = bnd_fixed->begin();
00137     std::copy( bnd_from_itr, bnd_fixed_from.end(), bnd_to_itr );
00138   }
00139   (*bnd_fixed)[nd_X] = LOWER; // relaxation is always initially active
00140   // j_f_decomp[]
00141   j_f_decomp->resize(0);
00142   // b_X
00143   b_X->resize(nd_X+1);
00144   if(nd_X) {
00145     // b_X[l-1] = { dL(i) if bnd_fixed[l-1] == LOWER or EQUALITY
00146     //              dU(i) if bnd_fixed[l-1] == UPPER }
00147     //             , l = 1...nd_X
00148     //             (where i = i_x_fixed[l-1])
00149     bnd_fixed_t::const_iterator
00150       bnd_itr     = const_cast<const bnd_fixed_t&>(*bnd_fixed).begin(),
00151       bnd_itr_end = const_cast<const bnd_fixed_t&>(*bnd_fixed).begin() + nd_X;
00152     i_x_fixed_t::const_iterator
00153       i_x_itr     = const_cast<const i_x_fixed_t&>(*i_x_fixed).begin();
00154     DVector::iterator
00155       b_X_itr     = b_X->begin();
00156     const SpVectorSlice::element_type
00157       *ele = NULL;
00158     for( ; bnd_itr != bnd_itr_end; ++bnd_itr, ++i_x_itr, ++b_X_itr ) {
00159       const size_type i = *i_x_itr;
00160       switch(*bnd_itr) {
00161           case LOWER:
00162           case EQUALITY:
00163           *b_X_itr = (ele = dL.lookup_element(i))->value(); // Should not be null!
00164           break;
00165           case UPPER:
00166           *b_X_itr = (ele = dU.lookup_element(i))->value(); // Should not be null!
00167           break;
00168           default:
00169           TEST_FOR_EXCEPT(true);
00170       }
00171     }
00172   }
00173   (*b_X)[nd_X] = etaL; // relaxation is always initially active
00174   // Ko = G.B_RR
00175   *Ko = G_super.B_RR_ptr(); // now B_RR is a shared object
00176   // fo = - G.Q_R'*g - op(G.B_RX)*b_X(1:nd_X)
00177   V_StMtV( fo, -1.0, Q_R, trans, g );
00178   if( nd_X && G_super.B_RX_ptr().get() )
00179     Vp_StMtV( &(*fo)(), -1.0, *G_super.B_RX_ptr(), G_super.B_RX_trans(), (*b_X)(1,nd_X) );
00180 
00181 }
00182 
00183 } // end namesapce ConstrainedOptPack

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