MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_EvalNewPointTailoredApproachOrthogonal_Step.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 "MoochoPack_EvalNewPointTailoredApproachOrthogonal_Step.hpp"
00043 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
00044 #include "NLPInterfacePack_NLPDirect.hpp"
00045 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00046 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00047 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00048 #include "AbstractLinAlgPack_VectorSpace.hpp"
00049 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00050 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00051 #include "AbstractLinAlgPack_AssertOp.hpp"
00052 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00053 #include "Teuchos_dyn_cast.hpp"
00054 #include "Teuchos_Assert.hpp"
00055 
00056 namespace MoochoPack {
00057 
00058 EvalNewPointTailoredApproachOrthogonal_Step::EvalNewPointTailoredApproachOrthogonal_Step(
00059   const deriv_tester_ptr_t                &deriv_tester
00060   ,const bounds_tester_ptr_t              &bounds_tester
00061   ,EFDDerivTesting                        fd_deriv_testing
00062   )
00063   :EvalNewPointTailoredApproach_Step(deriv_tester,bounds_tester,fd_deriv_testing)
00064 {}
00065 
00066 // protected
00067 
00068 void EvalNewPointTailoredApproachOrthogonal_Step::uninitialize_Y_Uy(
00069   MatrixOp         *Y
00070   ,MatrixOp        *Uy
00071   )
00072 {
00073   using Teuchos::dyn_cast;
00074 
00075   MatrixIdentConcatStd
00076     *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y)  : NULL;
00077   MatrixComposite
00078     *Uy_cpst = Uy ? &dyn_cast<MatrixComposite>(*Uy) : NULL;     
00079 
00080   if(Y_orth)
00081     Y_orth->set_uninitialized();
00082   TEUCHOS_TEST_FOR_EXCEPT( !( Uy_cpst == NULL ) ); // ToDo: Implement for undecomposed equalities
00083 }
00084 
00085 void EvalNewPointTailoredApproachOrthogonal_Step::calc_py_Y_Uy(
00086   const NLPDirect       &nlp
00087   ,const D_ptr_t        &D
00088   ,VectorMutable        *py
00089   ,MatrixOp             *Y
00090   ,MatrixOp             *Uy
00091   ,EJournalOutputLevel  olevel
00092   ,std::ostream         &out
00093   )
00094 {
00095   namespace rcp = MemMngPack;
00096   using Teuchos::dyn_cast;
00097   using LinAlgOpPack::syrk;
00098 
00099   const size_type
00100     n = nlp.n(),
00101     r = nlp.r();
00102   const Range1D
00103     var_dep(1,r),
00104     var_indep(r+1,n),
00105     con_decomp   = nlp.con_decomp(),
00106     con_undecomp = nlp.con_undecomp();
00107 
00108   //
00109   // Get pointers to concreate matrices
00110   //
00111   
00112   MatrixIdentConcatStd
00113     *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y)  : NULL;
00114   MatrixComposite
00115     *Uy_cpst = Uy ? &dyn_cast<MatrixComposite>(*Uy) : NULL;     
00116 
00117   //
00118   // Initialize the matrices
00119   //
00120 
00121   // Y
00122   if(Y_orth) {
00123     D_ptr_t  D_ptr = D;
00124 //    if(mat_rel == MATRICES_INDEP_IMPS) {
00125 //      D_ptr = D->clone();
00126 //      TEUCHOS_TEST_FOR_EXCEPTION(
00127 //        D_ptr.get() == NULL, std::logic_error
00128 //        ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
00129 //        "The matrix class used for the direct sensitivity matrix D = inv(C)*N of type \'"
00130 //        << typeName(*D) << "\' must return return.get() != NULL from the clone() method "
00131 //        "since mat_rel == MATRICES_INDEP_IMPS!" );
00132 //    }
00133     Y_orth->initialize(
00134       nlp.space_x()                                     // space_cols
00135       ,nlp.space_x()->sub_space(var_dep)->clone()       // space_rows
00136       ,MatrixIdentConcatStd::BOTTOM                     // top_or_bottom
00137       ,-1.0                                             // alpha
00138       ,D_ptr                                            // D_ptr
00139       ,BLAS_Cpp::trans                                  // D_trans
00140       );
00141   }
00142 
00143   // S
00144   if(S_ptr_.get() == NULL) {
00145     S_ptr_ = nlp.factory_S()->create();
00146   }
00147   // S = I + (D)'*(D')'
00148   dyn_cast<MatrixSymInitDiag>(*S_ptr_).init_identity(D->space_rows());
00149   syrk(*D,BLAS_Cpp::trans,1.0,1.0,S_ptr_.get());
00150 
00151   TEUCHOS_TEST_FOR_EXCEPT( !( Uy_cpst == NULL ) ); // ToDo: Implement for undecomposed equalities
00152 
00153   recalc_py(*D,py,olevel,out);
00154 
00155 }
00156 
00157 void EvalNewPointTailoredApproachOrthogonal_Step::recalc_py(
00158   const MatrixOp           &D
00159   ,VectorMutable           *py
00160   ,EJournalOutputLevel     olevel
00161   ,std::ostream            &out
00162   )
00163 {
00164   using BLAS_Cpp::no_trans;
00165   using BLAS_Cpp::trans;
00166   using AbstractLinAlgPack::Vp_StMtV;
00167   using AbstractLinAlgPack::V_InvMtV;
00168   using LinAlgOpPack::V_MtV;
00169 
00170   const MatrixSymOpNonsing   &S = *S_ptr_;
00171 
00172   VectorSpace::vec_mut_ptr_t               // ToDo: make workspace!
00173     tIa = D.space_rows().create_member(),
00174     tIb = D.space_rows().create_member();
00175   //
00176   // py = -inv(R)*c
00177   // py = -((I - D*inv(S)*D')*inv(C))*c
00178   //    = -(I - D*inv(S)*D')*(-py)
00179   //    = py - D*inv(S)*D'*py
00180   //
00181   // =>
00182   //
00183   // tIa  = D'*py
00184   // tIb  = inv(S)*tIa
00185   // py   += -D*tIb
00186   //
00187   V_MtV( tIa.get(), D, trans, *py );              // tIa  = D'*py
00188   V_InvMtV( tIb.get(), S, no_trans, *tIa );       // tIb  = inv(S)*tIa
00189   Vp_StMtV( py, -1.0, D, no_trans, *tIb );        // py   += -D*tIb
00190 }
00191 
00192 void EvalNewPointTailoredApproachOrthogonal_Step::print_calc_py_Y_Uy(
00193   std::ostream& out, const std::string& L
00194   ) const
00195 {
00196   out
00197     << L << "*** Orthogonal decomposition\n"
00198     << L << "py = inv(I + D*D') * py <: space_range\n"
00199     << L << "Y = [ I ; -D' ] <: space_x|space_range\n"
00200     << L << "Uy = ???\n"
00201     ;
00202 }
00203 
00204 } // end namespace MoochoPack 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends