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

Generated on Wed May 12 21:51:09 2010 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by  doxygen 1.4.7