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

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