MOOCHO Version of the Day
NLPInterfacePack_ExampleNLPFirstOrder.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 <stdexcept>
00045 
00046 #include "NLPInterfacePack_ExampleNLPFirstOrder.hpp"
00047 #include "NLPInterfacePack_ExampleBasisSystem.hpp"
00048 #include "AbstractLinAlgPack_BasisSystemComposite.hpp"
00049 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00050 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00051 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00052 #include "ReleaseResource_ref_count_ptr.hpp"
00053 #include "Teuchos_dyn_cast.hpp"
00054 #include "Teuchos_Assert.hpp"
00055 
00056 namespace NLPInterfacePack {
00057 
00058 ExampleNLPFirstOrder::ExampleNLPFirstOrder(
00059   const VectorSpace::space_ptr_t&  vec_space
00060   ,value_type                      xo
00061   ,bool                            has_bounds
00062   ,bool                            dep_bounded
00063   )
00064   :ExampleNLPObjGrad(vec_space,xo,has_bounds,dep_bounded)
00065   ,initialized_(false)
00066 {
00067   basis_sys_ = Teuchos::rcp(
00068     new ExampleBasisSystem(
00069       this->space_x()
00070       ,this->var_dep()
00071       ,this->var_indep()
00072       )
00073     );
00074 }
00075 
00076 // Overridden public members from NLPFirstOrder
00077 
00078 void ExampleNLPFirstOrder::set_Gc(MatrixOp* Gc)
00079 {
00080   if(Gc) // Throw an exception if this matrix is the wrong type!
00081     Teuchos::dyn_cast<AbstractLinAlgPack::MatrixComposite>(*Gc);
00082   NLPFirstOrder::set_Gc(Gc);
00083 }
00084 
00085 const NLPFirstOrder::mat_fcty_ptr_t
00086 ExampleNLPFirstOrder::factory_Gc() const
00087 {
00088   return factory_Gc_;
00089 }
00090 
00091 const NLPFirstOrder::basis_sys_ptr_t
00092 ExampleNLPFirstOrder::basis_sys() const
00093 {
00094   return basis_sys_;
00095 }
00096 
00097 // Overridden public members from NLP
00098 
00099 void ExampleNLPFirstOrder::initialize(bool test_setup)
00100 {
00101   namespace rcp = MemMngPack;
00102 
00103   ExampleNLPObjGrad::initialize(test_setup);
00104   NLPFirstOrder::initialize(test_setup);
00105 
00106   factory_Gc_ = BasisSystemComposite::factory_Gc();
00107   
00108   initialized_ = true;
00109 }
00110 
00111 bool ExampleNLPFirstOrder::is_initialized() const
00112 {
00113   return initialized_;
00114 }
00115 
00116 // Overridden protected members from NLPFirstOrder
00117 
00118 void ExampleNLPFirstOrder::imp_calc_Gc(
00119   const Vector& x, bool newx, const FirstOrderInfo& first_order_info) const
00120 {
00121   namespace rcp = MemMngPack;
00122   using Teuchos::dyn_cast;
00123   using AbstractLinAlgPack::Vp_S; // Should not have to do this!
00124 
00125   const index_type
00126     n = this->n(),
00127     m = this->m();
00128   const Range1D
00129     var_dep   = this->var_dep(),
00130     var_indep = this->var_indep();
00131 
00132   // Get references to aggregate C and N matrices (if allocated)
00133   MatrixOpNonsing
00134     *C_aggr = NULL;
00135   MatrixOp
00136     *N_aggr = NULL;
00137   BasisSystemComposite::get_C_N(
00138     first_order_info.Gc, &C_aggr, &N_aggr ); // Will return NULLs if Gc is not initialized
00139 
00140   // Allocate C and N matrix objects if not done yet!
00141   Teuchos::RCP<MatrixOpNonsing>
00142     C_ptr = Teuchos::null;
00143   Teuchos::RCP<MatrixOp>
00144     N_ptr = Teuchos::null;
00145   if( C_aggr == NULL ) {
00146     const VectorSpace::space_ptr_t
00147       space_x  = this->space_x(),
00148       space_xD = space_x->sub_space(var_dep);
00149     C_ptr  = Teuchos::rcp(new MatrixSymDiagStd(space_xD->create_member()));
00150     N_ptr  = Teuchos::rcp(new MatrixSymDiagStd(space_xD->create_member()));
00151     C_aggr = C_ptr.get();
00152     N_aggr = N_ptr.get();
00153   }
00154 
00155   // Get references to concreate C and N matrices
00156   MatrixSymDiagStd
00157     &C = dyn_cast<MatrixSymDiagStd>(*C_aggr);
00158   MatrixSymDiagStd
00159     &N = dyn_cast<MatrixSymDiagStd>(*N_aggr);
00160   // Get x = [ x_D' x_I ]
00161   Vector::vec_ptr_t
00162     x_D = x.sub_view(var_dep),
00163     x_I = x.sub_view(var_indep);
00164   // Set the diagonals of C and N (this is the only computation going on here)
00165   C.diag() = *x_I;          // C.diag = x_I - 1.0
00166   Vp_S( &C.diag(),  -1.0 ); // ...
00167   N.diag() = *x_D;          // N.diag = x_D - 10.0
00168   Vp_S( &N.diag(), -10.0 ); // ...
00169   // Initialize the matrix object Gc if not done so yet
00170   if( C_ptr.get() != NULL ) {
00171     BasisSystemComposite::initialize_Gc(
00172       this->space_x(), var_dep, var_indep
00173       ,this->space_c()
00174       ,C_ptr, N_ptr
00175       ,first_order_info.Gc
00176       );
00177   }
00178 }
00179 
00180 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends