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

Generated on Tue Jul 13 09:30:51 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7