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