ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_DecompositionSystemVarReductPermStd.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 "ConstrainedOptPack_DecompositionSystemVarReductPermStd.hpp"
00043 #include "ConstrainedOptPack_DecompositionSystemVarReductImp.hpp"
00044 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00045 #include "AbstractLinAlgPack_BasisSystemPerm.hpp"
00046 #include "AbstractLinAlgPack_PermutationOut.hpp"
00047 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00048 #include "Teuchos_Assert.hpp"
00049 
00050 namespace ConstrainedOptPack {
00051 
00052 // Constructors / initializers
00053 
00054 DecompositionSystemVarReductPermStd::DecompositionSystemVarReductPermStd(
00055   const decomp_sys_imp_ptr_t&        decomp_sys_imp
00056   ,const basis_sys_ptr_t&            basis_sys
00057   ,bool                              basis_selected
00058   ,EExplicitImplicit                 D_imp
00059   ,EExplicitImplicit                 Uz_imp
00060   )
00061 {
00062   this->initialize(decomp_sys_imp,basis_sys,basis_selected,D_imp,Uz_imp);
00063 }
00064 
00065 void DecompositionSystemVarReductPermStd::initialize(
00066   const decomp_sys_imp_ptr_t&        decomp_sys_imp
00067   ,const basis_sys_ptr_t&            basis_sys
00068   ,bool                              basis_selected
00069   ,EExplicitImplicit                 D_imp
00070   ,EExplicitImplicit                 Uz_imp
00071   )
00072 {
00073   decomp_sys_imp_ = decomp_sys_imp;
00074   basis_sys_      = basis_sys;
00075   basis_selected_ = basis_selected;
00076   this->D_imp(D_imp);
00077   this->Uz_imp(Uz_imp);
00078 }
00079 
00080 // Overridden from DecompositionSystem
00081 
00082 size_type DecompositionSystemVarReductPermStd::n() const
00083 {
00084   return decomp_sys_imp()->n();
00085 }
00086 
00087 size_type DecompositionSystemVarReductPermStd::m() const
00088 {
00089   return decomp_sys_imp()->m();
00090 }
00091 
00092 size_type DecompositionSystemVarReductPermStd::r() const
00093 {
00094   return decomp_sys_imp()->r();
00095 }
00096 
00097 Range1D DecompositionSystemVarReductPermStd::equ_decomp() const
00098 {
00099   return decomp_sys_imp()->equ_decomp();
00100 }
00101 
00102 Range1D DecompositionSystemVarReductPermStd::equ_undecomp() const
00103 {
00104   return decomp_sys_imp()->equ_undecomp();
00105 }
00106 
00107 const VectorSpace::space_ptr_t
00108 DecompositionSystemVarReductPermStd::space_range() const
00109 {
00110   return decomp_sys_imp()->space_range();
00111 }
00112 
00113 const VectorSpace::space_ptr_t
00114 DecompositionSystemVarReductPermStd::space_null() const
00115 {
00116   return decomp_sys_imp()->space_null();
00117 }
00118 
00119 const DecompositionSystem::mat_fcty_ptr_t
00120 DecompositionSystemVarReductPermStd::factory_Z() const
00121 {
00122   return decomp_sys_imp()->factory_Z();
00123 }
00124 
00125 const DecompositionSystem::mat_fcty_ptr_t
00126 DecompositionSystemVarReductPermStd::factory_Y() const
00127 {
00128   return decomp_sys_imp()->factory_Y();
00129 }
00130 
00131 const DecompositionSystem::mat_nonsing_fcty_ptr_t
00132 DecompositionSystemVarReductPermStd::factory_R() const
00133 {
00134   mat_nonsing_fcty_ptr_t factory_R = decomp_sys_imp()->factory_R();
00135   if( factory_R.get() != NULL )
00136     return factory_R;
00137   // Else assume that R will just be the basis matrix (coordinate decomposition!)
00138   return basis_sys_->factory_C();
00139 }
00140 
00141 const DecompositionSystem::mat_fcty_ptr_t
00142 DecompositionSystemVarReductPermStd::factory_Uz() const
00143 {
00144   return decomp_sys_imp()->factory_Uz();
00145 }
00146 
00147 const DecompositionSystem::mat_fcty_ptr_t
00148 DecompositionSystemVarReductPermStd::factory_Uy() const
00149 {
00150   return decomp_sys_imp()->factory_Uy();
00151 }
00152 
00153 void DecompositionSystemVarReductPermStd::update_decomp(
00154   std::ostream          *out
00155   ,EOutputLevel         olevel
00156   ,ERunTests            test_what
00157   ,const MatrixOp       &Gc
00158   ,MatrixOp             *Z
00159   ,MatrixOp             *Y
00160   ,MatrixOpNonsing      *R
00161   ,MatrixOp             *Uz
00162   ,MatrixOp             *Uy
00163   ,EMatRelations        mat_rel
00164   ) const
00165 {
00166   assert_basis_selected();
00167   decomp_sys_imp()->update_decomp(
00168     out,olevel,test_what,Gc,Z,Y
00169     ,R,Uz,Uy,mat_rel
00170     );
00171 }
00172 
00173 void DecompositionSystemVarReductPermStd::print_update_decomp(
00174   std::ostream& out, const std::string& L ) const
00175 {
00176   // ToDo: Print basis permutation stuff also?
00177   decomp_sys_imp()->print_update_decomp(out,L);
00178 }
00179 
00180 // Overridden from DecompositionSystemVarReduct
00181 
00182 Range1D DecompositionSystemVarReductPermStd::var_indep() const
00183 {
00184   return basis_sys_.get() ? basis_sys_->var_indep() : Range1D::Invalid;
00185 }
00186 
00187 Range1D DecompositionSystemVarReductPermStd::var_dep() const
00188 {
00189   return basis_sys_.get() ? basis_sys_->var_dep() : Range1D::Invalid;
00190 }
00191 
00192 // @name Overridden from DecompositionSystemVarReductPerm
00193 
00194 const DecompositionSystemVarReductPerm::perm_fcty_ptr_t
00195 DecompositionSystemVarReductPermStd::factory_P_var() const
00196 {
00197   return basis_sys()->factory_P_var();
00198 }
00199 
00200 const DecompositionSystemVarReductPerm::perm_fcty_ptr_t
00201 DecompositionSystemVarReductPermStd::factory_P_equ() const
00202 {
00203   return basis_sys()->factory_P_equ();
00204 }
00205 
00206 bool DecompositionSystemVarReductPermStd::has_basis() const
00207 {
00208   return basis_selected_;
00209 }
00210 
00211 void DecompositionSystemVarReductPermStd::set_decomp(
00212   std::ostream          *out
00213   ,EOutputLevel         olevel
00214   ,ERunTests            test_what
00215   ,const Permutation    &P_var
00216   ,const Range1D        &var_dep
00217   ,const Permutation    *P_equ
00218   ,const Range1D        *equ_decomp
00219   ,const MatrixOp       &Gc
00220   ,MatrixOp             *Z
00221   ,MatrixOp             *Y
00222   ,MatrixOpNonsing      *R
00223   ,MatrixOp             *Uz
00224   ,MatrixOp             *Uy
00225   ,EMatRelations        mat_rel
00226   )
00227 {
00228   // Forward these setting on to the implementation.
00229   decomp_sys_imp_->D_imp(  this->D_imp()  );
00230   decomp_sys_imp_->Uz_imp( this->Uz_imp() );
00231   // Get smart pointers to the basis matrix and the direct sensistivity matrices
00232   // and remove references to these matrix objects from the other decomposition
00233   // matrices by uninitializing them.
00234   Teuchos::RCP<MatrixOpNonsing>  C_ptr;
00235   Teuchos::RCP<MatrixOp>         D_ptr;
00236   decomp_sys_imp_->get_basis_matrices(
00237     out, olevel, test_what
00238     ,Z, Y, R, Uz, Uy
00239     ,&C_ptr
00240     ,&D_ptr // May return D_ptr.get() == NULL if not explicit chosen
00241     );
00242   // Tell the basis system object to set this basis
00243   try {
00244     basis_sys_->set_basis(
00245       P_var, var_dep
00246       ,P_equ, equ_decomp
00247       ,Gc
00248       ,C_ptr.get()
00249       ,D_ptr.get() // May be NULL
00250       ,this->Uz_imp() == MAT_IMP_EXPLICIT ? Uz : NULL
00251       ,(mat_rel == MATRICES_INDEP_IMPS
00252         ? BasisSystem::MATRICES_INDEP_IMPS : BasisSystem::MATRICES_ALLOW_DEP_IMPS )
00253       ,out
00254       );
00255   }
00256   catch( const BasisSystem::SingularBasis& except ) {
00257     if(out && olevel >= PRINT_BASIC_INFO)
00258       *out << "Passed in basis is singular, throwing SingularDecomposition: "
00259          << except.what() << std::endl;
00260     TEUCHOS_TEST_FOR_EXCEPTION(
00261       true, SingularDecomposition
00262       ,"DecompositionSystemVarReductPermStd::set_decomp(...): Passed in basis selection "
00263       "gave a singular basis matrix! : " << except.what() );
00264   }
00265   // If we get here the passed in basis selection is nonsingular and the basis matrices
00266   // are updated.  Now give them back to the decomp_sys_imp object and update the rest
00267   // of the decomposition matrices.
00268   const size_type
00269     m  = Gc.cols(),
00270     r  = C_ptr->rows();
00271   decomp_sys_imp_->set_basis_matrices(
00272     out, olevel, test_what
00273     ,C_ptr
00274     ,D_ptr // D_ptr.get() may be NULL
00275     ,r > m ? Uz : NULL
00276     ,basis_sys_ // Always reset
00277     );
00278   C_ptr = Teuchos::null;
00279   D_ptr = Teuchos::null;
00280   decomp_sys_imp()->update_decomp(
00281     out,olevel,test_what,Gc,Z,Y,R
00282     ,r > m ? Uz : NULL
00283     ,r > m ? Uy : NULL
00284     ,mat_rel
00285     );
00286   // We have a basis!
00287   basis_selected_ = true;
00288 }
00289 
00290 void DecompositionSystemVarReductPermStd::select_decomp(
00291   std::ostream          *out
00292   ,EOutputLevel         olevel
00293   ,ERunTests            test_what
00294   ,const Vector         *nu
00295   ,MatrixOp             *Gc
00296   ,Permutation          *P_var
00297   ,Range1D              *var_dep
00298   ,Permutation          *P_equ
00299   ,Range1D              *equ_decomp
00300   ,MatrixOp             *Z
00301   ,MatrixOp             *Y
00302   ,MatrixOpNonsing      *R
00303   ,MatrixOp             *Uz
00304   ,MatrixOp             *Uy
00305   ,EMatRelations        mat_rel
00306   )
00307 {
00308   // Forward these setting on to the implementation.
00309   decomp_sys_imp_->D_imp(  this->D_imp()  );
00310   decomp_sys_imp_->Uz_imp( this->Uz_imp() );
00311   // Get smart pointers to the basis matrix and the direct sensistivity matrices
00312   // and remove references to these matrix objects from the other decomposition
00313   // matrices by uninitializing them.
00314   Teuchos::RCP<MatrixOpNonsing>  C_ptr;
00315   Teuchos::RCP<MatrixOp>             D_ptr;
00316   //const bool unintialized_basis = decomp_sys_imp_->basis_sys()->var_dep().size() == 0;
00317   decomp_sys_imp_->get_basis_matrices(
00318     out, olevel, test_what
00319     ,Z, Y, R, Uz, Uy
00320     ,&C_ptr
00321     ,&D_ptr // May return D_ptr.get() == NULL if not explicit chosen
00322     );
00323   // Ask the basis system object to select a basis
00324   basis_sys_->select_basis(
00325     nu
00326     ,Gc
00327     ,P_var, var_dep
00328     ,P_equ, equ_decomp
00329     ,C_ptr.get()
00330     ,D_ptr.get() // May be NULL
00331     ,this->Uz_imp() == MAT_IMP_EXPLICIT ? Uz : NULL
00332     ,(mat_rel == MATRICES_INDEP_IMPS
00333       ? BasisSystem::MATRICES_INDEP_IMPS : BasisSystem::MATRICES_ALLOW_DEP_IMPS )
00334     ,out
00335     );
00336 
00337   if( out && (int)olevel >= (int)PRINT_BASIC_INFO ) {
00338     const Range1D var_indep = basis_sys_->var_indep(), equ_undecomp = basis_sys_->equ_undecomp();
00339     *out
00340       << "\nSelected a new basis\n"
00341       << "\nbs.var_dep()            = ["<<var_dep->lbound()<<","<<var_dep->ubound()<<"]"
00342       << "\nds.var_indep()          = ["<<var_indep.lbound()<<","<<var_indep.ubound()<<"]"
00343       << "\nds.equ_decomp()         = ["<<equ_decomp->lbound()<<","<<equ_decomp->ubound()<<"]"
00344       << "\nds.equ_undecomp()       = ["<<equ_undecomp.lbound()<<","<<equ_undecomp.ubound()<<"]"
00345       << std::endl;
00346   }
00347   if( out && (int)olevel >= (int)PRINT_VECTORS ) {
00348     *out
00349       << "\nP_var =\n" << *P_var
00350       << "\nP_equ =\n" << *P_equ
00351       ;
00352   }
00353   if( out && (int)olevel >= (int)PRINT_EVERY_THING ) {
00354     *out
00355       << "\nGc =\n" << *Gc;
00356   }
00357 
00358   // If we get here a nonsinguar basis selection has been made and the basis matrices
00359   // are updated.  Now give them back to the decomp_sys_imp object and update the rest
00360   // of the decomposition matrices.
00361   const size_type
00362     //n  = Gc->rows(),
00363     m  = Gc->cols(),
00364     r  = C_ptr->rows();
00365   decomp_sys_imp_->set_basis_matrices(
00366     out, olevel, test_what
00367     ,C_ptr
00368     ,D_ptr // D_ptr.get() may be NULL
00369     ,r > m ? Uz : NULL
00370     ,basis_sys_ // Always reset
00371     );
00372   C_ptr = Teuchos::null;
00373   D_ptr = Teuchos::null;
00374   decomp_sys_imp()->update_decomp(
00375     out,olevel,test_what,*Gc,Z,Y,R
00376     ,r > m ? Uz : NULL
00377     ,r > m ? Uy : NULL
00378     ,mat_rel
00379     );
00380   // We have a basis!
00381   basis_selected_ = true;
00382 }
00383 
00384 // private
00385 
00386 void DecompositionSystemVarReductPermStd::assert_basis_selected() const
00387 {
00388   TEUCHOS_TEST_FOR_EXCEPTION(
00389     !basis_selected_, std::logic_error
00390     ,"DecompositionSystemVarReductPermStd::assert_basis_selected(): Error, "
00391     "the methods set_decomp() or select_decomp() must be called first!" );
00392 }
00393 
00394 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends