ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_MatrixDecompRangeOrthog.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_MatrixDecompRangeOrthog.hpp"
00043 #include "AbstractLinAlgPack_VectorSpace.hpp"
00044 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00045 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
00046 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00047 #include "AbstractLinAlgPack_AssertOp.hpp"
00048 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00049 #include "Teuchos_Assert.hpp"
00050 #include "Teuchos_FancyOStream.hpp"
00051 
00052 namespace ConstrainedOptPack {
00053 
00054 // Constructors/initializers
00055 
00056 MatrixDecompRangeOrthog::MatrixDecompRangeOrthog()
00057 {}
00058 
00059 MatrixDecompRangeOrthog::MatrixDecompRangeOrthog(
00060   const C_ptr_t   &C_ptr
00061   ,const D_ptr_t  &D_ptr
00062   ,const S_ptr_t  &S_ptr
00063   )
00064 {
00065   this->initialize(C_ptr,D_ptr,S_ptr);
00066 }
00067 
00068 void MatrixDecompRangeOrthog::initialize(
00069   const C_ptr_t   &C_ptr
00070   ,const D_ptr_t  &D_ptr
00071   ,const S_ptr_t  &S_ptr
00072   )
00073 {
00074   const char func_name[] = "MatrixDecompRangeOrthog::initialize(...)";
00075   TEUCHOS_TEST_FOR_EXCEPTION(
00076     C_ptr.get() == NULL, std::invalid_argument
00077     ,func_name << " : Error!" );
00078   TEUCHOS_TEST_FOR_EXCEPTION(
00079     D_ptr.get() == NULL, std::invalid_argument
00080     ,func_name << " : Error!" );
00081   TEUCHOS_TEST_FOR_EXCEPTION(
00082     S_ptr.get() == NULL, std::invalid_argument
00083     ,func_name << " : Error!" );
00084 #ifdef ABSTRACT_LIN_ALG_PACK_CHECK_VEC_SPCS
00085   bool is_compatible = C_ptr->space_rows().is_compatible(D_ptr->space_cols());
00086   TEUCHOS_TEST_FOR_EXCEPTION(
00087     !is_compatible, VectorSpace::IncompatibleVectorSpaces
00088     ,func_name << " : Error, C_ptr->space_rows().is_compatible(D_ptr->space_cols()) == false!" );
00089   is_compatible = S_ptr->space_cols().is_compatible(D_ptr->space_rows());
00090   TEUCHOS_TEST_FOR_EXCEPTION(
00091     !is_compatible, VectorSpace::IncompatibleVectorSpaces
00092     ,func_name << " : Error, S_ptr->space_cols().is_compatible(D_ptr->space_rows()) == false!" );
00093 #endif  
00094   C_ptr_ = C_ptr;
00095   D_ptr_ = D_ptr;
00096   S_ptr_ = S_ptr;
00097 }
00098 
00099 void MatrixDecompRangeOrthog::set_uninitialized()
00100 {
00101   namespace rcp = MemMngPack;
00102   C_ptr_ = Teuchos::null;
00103   D_ptr_ = Teuchos::null;
00104   S_ptr_ = Teuchos::null;
00105 }
00106 
00107 // Overridden from MatrixOp
00108 
00109 size_type MatrixDecompRangeOrthog::rows() const
00110 {
00111   return C_ptr_.get() ? C_ptr_->rows() : 0;
00112 }
00113 
00114 size_type MatrixDecompRangeOrthog::cols() const
00115 {
00116   return C_ptr_.get() ? C_ptr_->cols() : 0;
00117 }
00118 
00119 const VectorSpace& MatrixDecompRangeOrthog::space_cols() const
00120 {
00121   return C_ptr_->space_cols();
00122 }
00123 
00124 const VectorSpace& MatrixDecompRangeOrthog::space_rows() const
00125 {
00126   return C_ptr_->space_rows();
00127 }
00128 
00129 std::ostream& MatrixDecompRangeOrthog::output(std::ostream& out_arg) const
00130 {
00131   Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00132   Teuchos::OSTab tab(out);
00133   *out
00134     << "This is a " << this->rows() << " x " << this->cols()
00135     << " nonsingular matrix (I + D'*D)*C with inverse inv(C')*(I + D*inv(S)*D') where C, D and S are:\n";
00136   tab.incrTab();
00137   *out << "C =\n" << *C_ptr();
00138   *out << "D =\n" << *D_ptr();
00139   *out << "S =\n" << *S_ptr();
00140   return out_arg;
00141 }
00142 
00143 void MatrixDecompRangeOrthog::Vp_StMtV(
00144   VectorMutable* y, value_type a, BLAS_Cpp::Transp R_trans
00145   , const Vector& x, value_type b
00146   ) const
00147 {
00148   using BLAS_Cpp::no_trans;
00149   using BLAS_Cpp::trans;
00150   using AbstractLinAlgPack::Vp_StV;
00151   using AbstractLinAlgPack::Vp_StMtV;
00152   using LinAlgOpPack::Vp_V;
00153   using LinAlgOpPack::V_MtV;
00154   using LinAlgOpPack::Vp_MtV;
00155   using LinAlgOpPack::V_StMtV;
00156 
00157   assert_initialized("MatrixDecompRangeOrthog::Vp_StMtV(...)");
00158   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,R_trans,x);
00159 
00160   const MatrixOpNonsing      &C = *C_ptr_;
00161   const MatrixOp             &D = *D_ptr_;
00162   const MatrixSymOpNonsing   &S = *S_ptr_;
00163   //
00164   // y = b*y + a*op(R)*x
00165   //
00166   VectorSpace::vec_mut_ptr_t               // ToDo: make workspace!
00167     tI = D.space_rows().create_member(),
00168     tD = D.space_cols().create_member();
00169   if(R_trans == no_trans) {
00170     //
00171     // y = b*y + a*C*(I + D*D')*x
00172     //
00173     // =>
00174     //
00175     // tI  = D'*x
00176     // tD  = x + D*tI
00177     // y   = b*y + a*C*tD
00178     //
00179     V_MtV( tI.get(), D, trans, x );          // tI  = D'*x
00180     *tD = x;                                 // tD = x + D*tI
00181     Vp_MtV( tD.get(), D, no_trans, *tI );    // ""
00182     Vp_StMtV( y, a, C, no_trans, *tD, b );   // y  = b*y + a*C*tD
00183   }
00184   else {
00185     //
00186     // y = b*y + a*(I + D*D')*C'*x
00187     //   = b*y + a*C'*x + D*D'*(a*C'*x)
00188     //
00189     // =>
00190     //
00191     // tD  = a*C'*x
00192     // tI  = D'*tD
00193     // y   = b*y + D*tI
00194     // y  += tD
00195     //
00196     V_StMtV( tD.get(), a, C, trans, x );      // tD  = a*C'*x
00197     V_MtV( tI.get(), D, trans, *tD );         // tI  = D'*tD
00198     Vp_MtV( y, D, no_trans, *tI, b );         // y   = b*y + D*tI
00199     Vp_V( y, *tD );                           // y  += tD
00200   }
00201 }
00202 
00203 // Overridden from MatrixOpNonsing
00204 
00205 void MatrixDecompRangeOrthog::V_InvMtV(
00206   VectorMutable* y, BLAS_Cpp::Transp R_trans
00207   , const Vector& x
00208   ) const
00209 {
00210   using BLAS_Cpp::no_trans;
00211   using BLAS_Cpp::trans;
00212   using AbstractLinAlgPack::Vp_StV;
00213   using AbstractLinAlgPack::Vp_StMtV;
00214   using AbstractLinAlgPack::V_InvMtV;
00215   using LinAlgOpPack::Vp_V;
00216   using LinAlgOpPack::V_MtV;
00217   using LinAlgOpPack::V_StMtV;
00218 
00219   assert_initialized("MatrixDecompRangeOrthog::V_InvMtV(...)");
00220   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,BLAS_Cpp::trans_not(R_trans),x);
00221 
00222   const MatrixOpNonsing      &C = *C_ptr_;
00223   const MatrixOp             &D = *D_ptr_;
00224   const MatrixSymOpNonsing   &S = *S_ptr_;
00225   //
00226   // y = inv(op(R))*x
00227   //
00228   VectorSpace::vec_mut_ptr_t               // ToDo: make workspace!
00229     tIa = D.space_rows().create_member(),
00230     tIb = D.space_rows().create_member(),
00231     tD  = D.space_cols().create_member();
00232   if(R_trans == no_trans) {
00233     //
00234     // y = (I - D*inv(S)*D')*inv(C)*x
00235     //   = inv(C)*x - D*inv(S)*D'*inv(C)*x
00236     //
00237     // =>
00238     //
00239     // y    = inv(C)*x
00240     // tIa  = D'*y
00241     // tIb  = inv(S)*tIa
00242     // y   += -D*tIb
00243     //
00244     V_InvMtV( y, C, no_trans, x );            // y    = inv(C)*x
00245     V_MtV( tIa.get(), D, trans, *y );         // tIa  = D'*y
00246     V_InvMtV( tIb.get(), S, no_trans, *tIa ); // tIb  = inv(S)*tIa
00247     Vp_StMtV( y, -1.0, D, no_trans, *tIb );   // y   += -D*tIb
00248   }
00249   else {
00250     //
00251     // y = inv(C')*(I - D*inv(S)*D')*x
00252     //
00253     // =>
00254     //
00255     // tIa  = D'*x
00256     // tIb  = inv(S)*tIa
00257     // tD   = x
00258     // tD  += -D*tIb
00259     // y    = Inv(C')*tD
00260     //
00261     V_MtV( tIa.get(), D, trans, x );                // tIa  = D'*x
00262     V_InvMtV( tIb.get(), S, no_trans, *tIa );       // tIb  = inv(S)*tIa
00263     *tD = x;                                        // tD   = x
00264     Vp_StMtV( tD.get(), -1.0, D, no_trans, *tIb );  // tD  += -D*tIb
00265     V_InvMtV( y, C, trans, *tD );                   // y    = Inv(C')*tD
00266   }
00267 }
00268 
00269 // private
00270 
00271 void MatrixDecompRangeOrthog::assert_initialized(const char func_name[]) const
00272 {
00273   TEUCHOS_TEST_FOR_EXCEPTION(
00274     C_ptr_.get() == NULL, std::logic_error
00275     ,func_name << " : Error, Must call initialize(...) first!" );
00276 }
00277 
00278 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends