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