ConstrainedOptPack_DecompositionSystemOrthogonal.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 
00033 #include "ConstrainedOptPack_DecompositionSystemOrthogonal.hpp"
00034 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
00035 #include "ConstrainedOptPack_MatrixDecompRangeOrthog.hpp"
00036 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00038 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00039 #include "AbstractLinAlgPack_MatrixOpSubView.hpp"
00040 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00041 #include "Teuchos_AbstractFactoryStd.hpp"
00042 #include "Teuchos_dyn_cast.hpp"
00043 #include "Teuchos_TestForException.hpp"
00044 
00045 namespace ConstrainedOptPack {
00046 
00047 DecompositionSystemOrthogonal::DecompositionSystemOrthogonal(
00048   const VectorSpace::space_ptr_t           &space_x
00049   ,const VectorSpace::space_ptr_t          &space_c
00050   ,const basis_sys_ptr_t                   &basis_sys
00051   ,const basis_sys_tester_ptr_t            &basis_sys_tester
00052   ,EExplicitImplicit                       D_imp
00053   ,EExplicitImplicit                       Uz_imp
00054   )
00055   :DecompositionSystemVarReductImp(
00056     space_x, space_c, basis_sys, basis_sys_tester
00057     ,D_imp, Uz_imp )
00058 {}
00059 
00060 // Overridden from DecompositionSystem
00061 
00062 const DecompositionSystem::mat_fcty_ptr_t
00063 DecompositionSystemOrthogonal::factory_Y() const
00064 {
00065   namespace rcp = MemMngPack;
00066   return Teuchos::rcp(
00067     new Teuchos::AbstractFactoryStd<MatrixOp,MatrixIdentConcatStd>()
00068     );
00069 }
00070 
00071 const DecompositionSystem::mat_nonsing_fcty_ptr_t
00072 DecompositionSystemOrthogonal::factory_R() const
00073 {
00074   namespace rcp = MemMngPack;
00075   return Teuchos::rcp(
00076     new Teuchos::AbstractFactoryStd<MatrixOpNonsing,MatrixDecompRangeOrthog>()
00077     );
00078 }
00079 
00080 const DecompositionSystem::mat_fcty_ptr_t
00081 DecompositionSystemOrthogonal::factory_Uy() const
00082 {
00083   return Teuchos::rcp(  new Teuchos::AbstractFactoryStd<MatrixOp,MatrixOpSubView>() );
00084 }
00085 
00086 // Overridden from DecompositionSystemVarReductImp
00087 
00088 void DecompositionSystemOrthogonal::update_D_imp_used(EExplicitImplicit *D_imp_used) const
00089 {
00090   *D_imp_used = MAT_IMP_EXPLICIT;
00091 }
00092 
00093 DecompositionSystem::mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t
00094 DecompositionSystemOrthogonal::uninitialize_matrices(
00095   std::ostream                                       *out
00096   ,EOutputLevel                                      olevel
00097   ,MatrixOp                                          *Y
00098   ,MatrixOpNonsing                                   *R
00099   ,MatrixOp                                          *Uy
00100   ) const
00101 {
00102   namespace rcp = MemMngPack;
00103   using Teuchos::dyn_cast;
00104   typedef DecompositionSystem::mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t
00105     C_ptr_t;
00106 
00107   //
00108   // Get pointers to concreate matrices
00109   //
00110   
00111   MatrixIdentConcatStd
00112     *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y)    : NULL;
00113   MatrixDecompRangeOrthog
00114     *R_orth = R ? &dyn_cast<MatrixDecompRangeOrthog>(*R) : NULL;
00115   MatrixOpSubView
00116     *Uy_cpst = Uy ? &dyn_cast<MatrixOpSubView>(*Uy)  : NULL;      
00117 
00118   //
00119   // Get the smart pointer to the basis matrix object C and the
00120   // matrix S = I + D'*D
00121   //
00122   
00123   C_ptr_t C_ptr = Teuchos::null;
00124   if(R_orth) {
00125     C_ptr  = Teuchos::rcp_const_cast<MatrixOpNonsing>(    R_orth->C_ptr() ); // This could be NULL!
00126     S_ptr_ = Teuchos::rcp_const_cast<MatrixSymOpNonsing>( R_orth->S_ptr() ); // ""
00127   }
00128   
00129   //
00130   // Uninitialize all of the matrices to remove references to C, D etc.
00131   //
00132 
00133   if(Y_orth)
00134     Y_orth->set_uninitialized();
00135   if(R_orth)
00136     R_orth->set_uninitialized();
00137   if(Uy_cpst)
00138     Uy_cpst->initialize(Teuchos::null);
00139 
00140   //
00141   // Return the owned? basis matrix object C
00142   //
00143 
00144   return C_ptr;
00145 
00146 }
00147 
00148 void DecompositionSystemOrthogonal::initialize_matrices(
00149   std::ostream                                           *out
00150   ,EOutputLevel                                          olevel
00151   ,const mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t &C
00152   ,const mat_fcty_ptr_t::element_type::obj_ptr_t         &D
00153   ,MatrixOp                                              *Y
00154   ,MatrixOpNonsing                                       *R
00155   ,MatrixOp                                              *Uy
00156   ,EMatRelations                                         mat_rel
00157   ) const
00158 {
00159   namespace rcp = MemMngPack;
00160   using Teuchos::dyn_cast;
00161   using LinAlgOpPack::syrk;
00162   typedef DecompositionSystem::mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t
00163     C_ptr_t;
00164   typedef DecompositionSystem::mat_fcty_ptr_t::element_type::obj_ptr_t
00165     D_ptr_t;
00166 
00167   const size_type
00168     //n = this->n(),
00169     r = this->r();
00170   const Range1D
00171     var_dep(1,r);
00172 
00173   //
00174   // Get pointers to concreate matrices
00175   //
00176   
00177   MatrixIdentConcatStd
00178     *Y_orth = Y ? &dyn_cast<MatrixIdentConcatStd>(*Y)    : NULL;
00179   MatrixDecompRangeOrthog
00180     *R_orth = R ? &dyn_cast<MatrixDecompRangeOrthog>(*R) : NULL;
00181 
00182   //
00183   // Initialize the matrices
00184   //
00185 
00186   if(Y_orth) {
00187     D_ptr_t  D_ptr = D;
00188     if(mat_rel == MATRICES_INDEP_IMPS) {
00189       D_ptr = D->clone();
00190       TEST_FOR_EXCEPTION(
00191         D_ptr.get() == NULL, std::logic_error
00192         ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
00193         "The matrix class used for the direct sensitivity matrix D = inv(C)*N of type \'"
00194         << typeName(*D) << "\' must return return.get() != NULL from the clone() method "
00195         "since mat_rel == MATRICES_INDEP_IMPS!" );
00196     }
00197     Y_orth->initialize(
00198       space_x()                                         // space_cols
00199       ,space_x()->sub_space(var_dep)->clone()           // space_rows
00200       ,MatrixIdentConcatStd::BOTTOM                     // top_or_bottom
00201       ,-1.0                                             // alpha
00202       ,D_ptr                                            // D_ptr
00203       ,BLAS_Cpp::trans                                  // D_trans
00204       );
00205   }
00206   if(R_orth) {
00207     C_ptr_t  C_ptr = C;
00208     if(mat_rel == MATRICES_INDEP_IMPS) {
00209       C_ptr = C->clone_mwons();
00210       TEST_FOR_EXCEPTION(
00211         C_ptr.get() == NULL, std::logic_error
00212         ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
00213         "The matrix class used for the basis matrix C of type \'"
00214         << typeName(*C) << "\' must return return.get() != NULL from the clone_mwons() method "
00215         "since mat_rel == MATRICES_INDEP_IMPS!" );
00216     }
00217     D_ptr_t  D_ptr = D;
00218     if(mat_rel == MATRICES_INDEP_IMPS) {
00219       D_ptr = D->clone();
00220       TEST_FOR_EXCEPTION(
00221         D_ptr.get() == NULL, std::logic_error
00222         ,"DecompositionSystemOrthogonal::update_decomp(...) : Error, "
00223         "The matrix class used for the direct sensitivity matrix D = inv(C)*N of type \'"
00224         << typeName(*D) << "\' must return return.get() != NULL from the clone() method "
00225         "since mat_rel == MATRICES_INDEP_IMPS!" );
00226     }
00227     if(S_ptr_.get() == NULL) {
00228       S_ptr_ = basis_sys()->factory_S()->create();
00229     }
00230     try {
00231       // S = I + (D)'*(D')'
00232       dyn_cast<MatrixSymInitDiag>(*S_ptr_).init_identity(D_ptr->space_rows());
00233       syrk(*D_ptr,BLAS_Cpp::trans,1.0,1.0,S_ptr_.get());
00234     }
00235     catch( const MatrixNonsing::SingularMatrix& except ) {
00236       TEST_FOR_EXCEPTION(
00237         true, SingularDecomposition
00238         ,"DecompositionSystemOrthogonal::initialize_matrices(...) : Error, update of S failed : "
00239         << except.what() );
00240     }
00241     R_orth->initialize(C_ptr,D_ptr,S_ptr_);
00242   }
00243   // ToDo: Implement for undecomposed equalities and general inequalities
00244 }
00245 
00246 void DecompositionSystemOrthogonal::print_update_matrices(
00247   std::ostream& out, const std::string& L ) const
00248 {
00249   out
00250     << L << "*** Orthogonal decompositon Y, R and Uy matrices (class DecompositionSystemOrthogonal)\n"
00251     << L << "Y  = [ I; -D' ] (using class MatrixIdentConcatStd)\n"
00252     << L << "R  = C*(I + D*D')\n"
00253     << L << "Uy = E - F*D'\n"
00254     ;
00255 }
00256     
00257 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:10:57 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3