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