NLPInterfacePack: C++ Interfaces and Implementation for Non-Linear Programs Version of the Day
NLPInterfacePack_NLPSerialPreprocessExplJac.cpp
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <assert.h>
00043 
00044 #include <typeinfo>
00045 #include <algorithm>
00046 
00047 #include "NLPInterfacePack_NLPSerialPreprocessExplJac.hpp"
00048 #include "AbstractLinAlgPack_MatrixPermAggr.hpp"
00049 #include "AbstractLinAlgPack_BasisSystemFactory.hpp"
00050 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00051 #include "AbstractLinAlgPack_MatrixSparseCOORSerial.hpp"
00052 #include "AbstractLinAlgPack_PermutationSerial.hpp"
00053 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00054 #include "DenseLinAlgPack_DVectorOp.hpp"
00055 #include "DenseLinAlgPack_IVector.hpp"
00056 #include "DenseLinAlgPack_PermVecMat.hpp"
00057 #include "Teuchos_Assert.hpp"
00058 #include "Teuchos_dyn_cast.hpp"
00059 #include "Teuchos_AbstractFactoryStd.hpp"
00060 #include "OptionsFromStreamPack_OptionsFromStream.hpp"
00061 
00062 namespace NLPInterfacePack {
00063 
00064 // NLPSerialPreprocessExplJac
00065 
00066 // Constructors / initializers
00067 
00068 NLPSerialPreprocessExplJac::NLPSerialPreprocessExplJac(
00069   const basis_sys_fcty_ptr_t  &basis_sys_fcty
00070   ,const factory_mat_ptr_t    &factory_Gc_full
00071   )
00072   :initialized_(false),test_setup_(false)
00073 {
00074   this->set_basis_sys_fcty(basis_sys_fcty);
00075   this->set_factory_Gc_full(factory_Gc_full);
00076 }
00077 
00078 void NLPSerialPreprocessExplJac::set_factory_Gc_full(
00079   const factory_mat_ptr_t     &factory_Gc_full
00080   )
00081 {
00082   if(factory_Gc_full.get())
00083     factory_Gc_full_ = factory_Gc_full;
00084   else 
00085     factory_Gc_full_ = Teuchos::rcp(
00086       new Teuchos::AbstractFactoryStd<MatrixOp,MatrixSparseCOORSerial>() );
00087   factory_Gc_ = Teuchos::rcp( new Teuchos::AbstractFactoryStd<MatrixOp,MatrixPermAggr>() );
00088 }
00089 
00090 // Overridden public members from NLP
00091 
00092 void NLPSerialPreprocessExplJac::set_options( const options_ptr_t& options )
00093 {
00094   options_ = options;
00095 }
00096 
00097 const NLP::options_ptr_t&
00098 NLPSerialPreprocessExplJac::get_options() const
00099 {
00100   return options_;
00101 }
00102 
00103 void NLPSerialPreprocessExplJac::initialize(bool test_setup)
00104 {
00105   namespace mmp = MemMngPack;
00106 
00107   test_setup_ = test_setup;
00108 
00109   if( initialized_  && !imp_nlp_has_changed() ) {
00110     // The subclass NLP has not changed so we can just
00111     // slip this preprocessing.
00112     NLPFirstOrder::initialize(test_setup);
00113     NLPSerialPreprocess::initialize(test_setup);  // Some duplication but who cares!
00114     return;
00115   }
00116 
00117   // Initialize the base object first
00118   NLPFirstOrder::initialize(test_setup);
00119   NLPSerialPreprocess::initialize(test_setup);  // Some duplication but who cares!
00120 
00121   // Initialize the storage for the intermediate quanities
00122   Gc_nz_orig_ = imp_Gc_nz_orig();         // Get the estimated number of nonzeros in Gc
00123   Gc_val_orig_.resize(Gc_nz_orig_);
00124   Gc_ivect_orig_.resize(Gc_nz_orig_);
00125   Gc_jvect_orig_.resize(Gc_nz_orig_);
00126   Gh_nz_orig_ = imp_Gh_nz_orig();     // Get the estimated number of nonzeros in Gh
00127   Gh_val_orig_.resize(Gh_nz_orig_);
00128   Gh_ivect_orig_.resize(Gh_nz_orig_);
00129   Gh_jvect_orig_.resize(Gh_nz_orig_);
00130 
00131   Gc_perm_new_basis_updated_ = false;
00132 
00133   // If you get here then the initialization went Ok.
00134   initialized_ = true;
00135 }
00136 
00137 bool NLPSerialPreprocessExplJac::is_initialized() const {
00138   return initialized_;
00139 }
00140 
00141 // Overridden public members from NLPFirstOrder
00142 
00143 const NLPFirstOrder::mat_fcty_ptr_t
00144 NLPSerialPreprocessExplJac::factory_Gc() const
00145 {
00146   return factory_Gc_;
00147 }
00148 
00149 const NLPFirstOrder::basis_sys_ptr_t
00150 NLPSerialPreprocessExplJac::basis_sys() const
00151 {
00152   BasisSystemFactory &fcty = const_cast<NLPSerialPreprocessExplJac*>(this)->basis_sys_fcty();
00153   fcty.set_options(options_);
00154   return fcty.create();
00155 }
00156 
00157 void NLPSerialPreprocessExplJac::set_Gc(MatrixOp* Gc)
00158 {
00159   using Teuchos::dyn_cast;
00160   assert_initialized();
00161   if( Gc != NULL ) {
00162     dyn_cast<MatrixPermAggr>(*Gc); // With throw exception if not correct type!
00163   }
00164   NLPFirstOrder::set_Gc(Gc);
00165 }
00166 
00167 // Overridden public members from NLPVarReductPerm
00168 
00169 bool NLPSerialPreprocessExplJac::get_next_basis(
00170   Permutation*  P_var,   Range1D* var_dep
00171   ,Permutation* P_equ,   Range1D* equ_decomp
00172   )
00173 {
00174   const bool new_basis = NLPSerialPreprocess::get_next_basis(
00175     P_var,var_dep,P_equ,equ_decomp
00176     );
00177   if( new_basis ) {
00178     Gc_perm_new_basis_updated_ = false;
00179   }
00180   return new_basis;
00181 }
00182 
00183 void NLPSerialPreprocessExplJac::set_basis(
00184   const Permutation   &P_var,   const Range1D  &var_dep
00185   ,const Permutation  *P_equ,   const Range1D  *equ_decomp
00186   )
00187 {
00188   NLPSerialPreprocess::set_basis(
00189     P_var,var_dep,P_equ,equ_decomp
00190     );
00191   Gc_perm_new_basis_updated_ = false;
00192 }
00193 
00194 // Overridden protected members from NLPFirstOrder
00195 
00196 void NLPSerialPreprocessExplJac::imp_calc_Gc(
00197   const Vector& x, bool newx
00198   ,const FirstOrderInfo& first_order_info
00199   ) const
00200 {
00201   namespace mmp = MemMngPack;
00202   using Teuchos::dyn_cast;
00203 
00204   assert_initialized();
00205 
00206   const Range1D
00207     var_dep      = this->var_dep(),
00208     equ_decomp   = this->equ_decomp();
00209   // Get the dimensions of the original NLP
00210   const size_type
00211     n       = this->n(),
00212     n_orig  = this->imp_n_orig(),
00213     m_orig  = this->imp_m_orig(),
00214     mI_orig = this->imp_mI_orig();
00215   // Get the dimensions of the full matrices
00216   const size_type
00217     n_full  = n_orig + mI_orig,
00218     m_full  = m_orig + mI_orig;
00219   // Get the number of columns in the matrix being constructed here
00220   const size_type
00221     num_cols = m_full;
00222 
00223   //
00224   // Get references to the constituent objects
00225   //
00226 
00227   // Get the concrete type for the Jacobian matrix
00228   MatrixPermAggr
00229     &G_aggr = dyn_cast<MatrixPermAggr>( *first_order_info.Gc );
00230   // Get smart pointers to the constituent members
00231   Teuchos::RCP<MatrixOp>
00232     G_full = Teuchos::rcp_const_cast<MatrixOp>( G_aggr.mat_orig() );
00233   Teuchos::RCP<PermutationSerial>
00234     P_row = Teuchos::rcp_dynamic_cast<PermutationSerial>(
00235       Teuchos::rcp_const_cast<Permutation>( G_aggr.row_perm() ) );  // variable permutation
00236   Teuchos::RCP<PermutationSerial>
00237     P_col = Teuchos::rcp_dynamic_cast<PermutationSerial>(
00238       Teuchos::rcp_const_cast<Permutation>( G_aggr.col_perm() ) );  // constraint permutation
00239   Teuchos::RCP<const MatrixOp>
00240     G_perm = G_aggr.mat_perm();
00241   // Remove references to G_full, G_perm, P_row and P_col.
00242   G_aggr.set_uninitialized();
00243   // Allocate the original matrix object if not done so yet
00244   if( G_full.get() == NULL || G_full.count() > 1 )
00245     G_full = factory_Gc_full_->create();
00246   // Get reference to the MatrixLoadSparseElements interface
00247   MatrixLoadSparseElements
00248     &G_lse = dyn_cast<MatrixLoadSparseElements>(*G_full);
00249 
00250   //
00251   // Calcuate the full explicit Jacobian
00252   //
00253 
00254   set_x_full( VectorDenseEncap(x)(), newx, &x_full() );
00255   if( m_orig )
00256     imp_calc_Gc_orig( x_full(), newx, first_order_expl_info() );
00257   if( mI_orig )
00258     imp_calc_Gh_orig( x_full(), newx, first_order_expl_info() );
00259 
00260   // Now get the actual number of nonzeros
00261   const size_type nz_full
00262     = Gc_nz_orig_ + Gh_nz_orig_ + mI_orig;  // Gc_orig, Gh_orig, -I
00263 
00264   // Determine if we need to set the structure and the nonzeros or just the nonzero values
00265   const bool load_struct = (G_lse.nz() == 0);
00266 
00267   size_type G_nz_previous;
00268   if( load_struct ) {
00269     G_lse.reinitialize(n,num_cols,nz_full); // The actual number of nonzeros will be minus the fixed variables
00270   }
00271   else {
00272     G_nz_previous = G_lse.nz();
00273     G_lse.reset_to_load_values();           // Use row and column indexes already set (better be same insert order!)
00274   }
00275 
00276   //
00277   // Load the matrix entries where we remove variables fixed by bounds
00278   //
00279 
00280   // Get pointers to buffers to fill with nonzero entries
00281   value_type      *val    = NULL;
00282    index_type     *ivect  = NULL,
00283                     *jvect  = NULL;
00284   G_lse.get_load_nonzeros_buffers(
00285     nz_full // We may actually load less
00286     ,&val
00287     ,load_struct ? &ivect : NULL
00288     ,load_struct ? &jvect : NULL
00289     );
00290   // Pointers to the full COOR matrix just updated
00291   const value_type      *val_orig     = NULL;
00292   const value_type      *val_orig_end = NULL;
00293    const index_type      *ivect_orig   = NULL;
00294   const index_type      *jvect_orig   = NULL;
00295 
00296   index_type nz = 0;
00297   if( m_orig ) {
00298     // Load entries for Gc_orig
00299     val_orig    = &Gc_val_orig_[0];
00300     val_orig_end  = val_orig + Gc_nz_orig_;
00301     ivect_orig    = &Gc_ivect_orig_[0];
00302     jvect_orig    = &Gc_jvect_orig_[0];
00303     imp_fill_jacobian_entries(
00304       n, n_full, load_struct
00305       ,0 // column offset
00306       ,val_orig, val_orig_end, ivect_orig, jvect_orig
00307       ,&nz // This will be incremented
00308       ,val, ivect, jvect
00309       );
00310   }
00311   if( mI_orig > 0 ) {
00312     // Load entires for Gc_orig and -I
00313     val_orig    = &Gh_val_orig_[0];
00314     val_orig_end  = val_orig + Gh_nz_orig_;
00315     ivect_orig    = &Gh_ivect_orig_[0];
00316     jvect_orig    = &Gh_jvect_orig_[0];
00317     imp_fill_jacobian_entries(
00318       n, n_full, load_struct
00319       ,m_orig // column offset (i.e. [ Gc_orig, Gh_orig ] )
00320       ,val_orig, val_orig_end, ivect_orig, jvect_orig
00321       ,&nz // This will be incremented
00322       ,val + nz, ivect + nz, jvect + nz
00323       );
00324     // -I
00325     value_type         *val_itr    = val   + nz;
00326     index_type         *ivect_itr  = ivect + nz;
00327     index_type         *jvect_itr  = jvect + nz;
00328     const IVector& var_full_to_remove_fixed = this->var_full_to_remove_fixed();
00329     if( load_struct ) {
00330       // Fill values and i and j
00331       for( index_type k = 1; k <= mI_orig; ++k ) {
00332         size_type var_idx = var_full_to_remove_fixed(n_orig+k); // Knows about slacks
00333 #ifdef TEUCHOS_DEBUG
00334         TEUCHOS_TEST_FOR_EXCEPT( !(  0 < var_idx && var_idx <= n_full  ) );
00335 #endif
00336         if(var_idx <= n) {
00337           // This is not a fixed variable
00338           *val_itr++ = -1.0;
00339           *ivect_itr++ = var_idx;
00340           *jvect_itr++ = m_orig + k; // (i.e. [ 0,  -I ] )
00341           ++nz;
00342         }
00343       }
00344     }
00345     else {
00346       // Just fill values
00347       for( index_type k = 1; k <= mI_orig; ++k ) {
00348         size_type var_idx = var_full_to_remove_fixed(n_orig+k); // Knows about slacks
00349 #ifdef TEUCHOS_DEBUG
00350         TEUCHOS_TEST_FOR_EXCEPT( !(  0 < var_idx && var_idx <= n_full  ) );
00351 #endif
00352         if(var_idx <= n) {
00353           // This is not a fixed variable
00354           *val_itr++ = -1.0;
00355           ++nz;
00356         }
00357       }
00358     }
00359   }
00360 
00361   if( !load_struct ) {
00362     // Check that the number of nonzeros added matches the number of nonzeros in G
00363     TEUCHOS_TEST_FOR_EXCEPTION(
00364       G_nz_previous != nz, std::runtime_error
00365       ,"NLPSerialPreprocessExplJac::imp_calc_Gc(...): Error, "
00366       "The number of added nonzeros does not match the number of nonzeros "
00367       "in the previous matrix load!." );
00368   }
00369   
00370   // Commit the nonzeros
00371   G_lse.commit_load_nonzeros_buffers(
00372     nz  // The actual number of nonzeros to set
00373     ,&val
00374     ,load_struct ? &ivect : NULL
00375     ,load_struct ? &jvect : NULL
00376     );
00377   G_lse.finish_construction(test_setup_);
00378 
00379   //
00380   // Setup permuted view
00381   //
00382 
00383   // Setup row (variable) permutation
00384   if( P_row.get() == NULL || P_col.count() > 1 )
00385       P_row = Teuchos::rcp(new PermutationSerial());
00386   Teuchos::RCP<IVector>        var_perm;
00387   if( P_row->perm().get() == NULL )  var_perm = Teuchos::rcp(new IVector(n_full));
00388   else                               var_perm = Teuchos::rcp_const_cast<IVector>(P_row->perm());
00389   *var_perm = this->var_perm();
00390   P_row->initialize(var_perm,Teuchos::null);
00391   // Setup column (constraint) permutation
00392   if( P_col.get() == NULL || P_col.count() > 1 )
00393       P_col = Teuchos::rcp(new PermutationSerial());
00394   Teuchos::RCP<IVector>        con_perm;
00395   if( P_col->perm().get() == NULL )  con_perm = Teuchos::rcp(new IVector(m_full));
00396   else                               con_perm = Teuchos::rcp_const_cast<IVector>(P_col->perm());
00397   *con_perm = this->equ_perm();
00398   P_col->initialize(con_perm,Teuchos::null);
00399   // Setup G_perm
00400   int num_row_part, num_col_part;
00401   index_type row_part[3], col_part[3];
00402   if(var_dep.size()) {
00403     num_row_part = 2;
00404     row_part[0] = 1;
00405     row_part[1] = (var_dep.lbound() == 1 ? var_dep.ubound()+1 : var_dep.lbound());
00406     row_part[2] = n+1;
00407   }
00408   else {
00409     num_row_part = 1;
00410     row_part[0] = 1;
00411     row_part[1] = n+1;
00412   }
00413   if(equ_decomp.size()) {
00414     num_col_part = 2;
00415     col_part[0] = 1;
00416     col_part[1] = (equ_decomp.lbound() == 1 ? equ_decomp.ubound()+1 : equ_decomp.lbound());
00417     col_part[2] = m_full+1;
00418   }
00419   else {
00420     num_col_part = 1;
00421     col_part[0] = 1;
00422     col_part[1] = m_full+1;
00423   }
00424   if( G_perm.get() == NULL || !Gc_perm_new_basis_updated_ ) {
00425     G_perm = G_full->perm_view(
00426       P_row.get(),row_part,num_row_part
00427       ,P_col.get(),col_part,num_col_part
00428       );
00429   }
00430   else {
00431     G_perm = G_full->perm_view_update(
00432       P_row.get(),row_part,num_row_part
00433       ,P_col.get(),col_part,num_col_part
00434       ,G_perm
00435       );
00436   }
00437   Gc_perm_new_basis_updated_ = true;
00438 
00439   //
00440   // Reinitialize the aggregate matrix object.
00441   //
00442 
00443   G_aggr.initialize(G_full,P_row,P_col,G_perm);
00444 }
00445 
00446 // protected members
00447 
00448 void NLPSerialPreprocessExplJac::assert_initialized() const
00449 {
00450   TEUCHOS_TEST_FOR_EXCEPTION(
00451     !initialized_, UnInitialized
00452     ,"NLPSerialPreprocessExplJac : The nlp has not been initialized yet" );
00453 }
00454 
00455 // private
00456 
00457 void NLPSerialPreprocessExplJac::imp_fill_jacobian_entries(
00458   size_type           n
00459   ,size_type          n_full
00460   ,bool               load_struct
00461   ,const index_type   col_offset
00462   ,const value_type   *val_orig
00463   ,const value_type   *val_orig_end
00464   ,const index_type   *ivect_orig
00465   ,const index_type   *jvect_orig
00466   ,index_type         *nz
00467   ,value_type         *val_itr
00468   ,index_type         *ivect_itr
00469   ,index_type         *jvect_itr
00470   ) const
00471 {
00472   const IVector& var_full_to_remove_fixed = this->var_full_to_remove_fixed();
00473   if( load_struct ) {
00474     // Fill values and i and j
00475     for( ; val_orig != val_orig_end ; ++val_orig, ++ivect_orig, ++jvect_orig) {
00476 #ifdef TEUCHOS_DEBUG
00477       TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= *ivect_orig && *ivect_orig <= n_full  ) );
00478 #endif
00479       size_type var_idx = var_full_to_remove_fixed(*ivect_orig);
00480 #ifdef TEUCHOS_DEBUG
00481       TEUCHOS_TEST_FOR_EXCEPT( !(  0 < var_idx && var_idx <= n_full  ) );
00482 #endif
00483       if(var_idx <= n) {
00484         // This is not a fixed variable
00485         *val_itr++ = *val_orig;
00486         // Also fill the row and column indices
00487         *ivect_itr++ = var_idx;
00488         *jvect_itr++ = col_offset + (*jvect_orig);
00489         ++(*nz);
00490       }
00491     }
00492   }
00493   else {
00494     // Just fill values
00495     for( ; val_orig != val_orig_end ; ++val_orig, ++ivect_orig) {
00496 #ifdef TEUCHOS_DEBUG
00497       TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= *ivect_orig && *ivect_orig <= n_full  ) );
00498 #endif
00499       size_type var_idx = var_full_to_remove_fixed(*ivect_orig);
00500 #ifdef TEUCHOS_DEBUG
00501       TEUCHOS_TEST_FOR_EXCEPT( !(  0 < var_idx && var_idx <= n_full  ) );
00502 #endif
00503       if(var_idx <= n) {
00504         // This is not a fixed variable
00505         *val_itr++ = *val_orig;
00506         ++(*nz);
00507       }
00508     }
00509   }
00510 }
00511 
00512 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends