MOOCHO (Single Doxygen Collection) Version of the Day
ConstrainedOptPack_MatrixIdentConcat.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 // 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_MatrixIdentConcat.hpp"
00043 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00044 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00045 #include "AbstractLinAlgPack_VectorSpace.hpp"
00046 #include "AbstractLinAlgPack_VectorMutable.hpp"
00047 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00048 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00049 #include "AbstractLinAlgPack_AssertOp.hpp"
00050 #include "AbstractLinAlgPack_SpVectorOp.hpp" // RAB: 12/20/2002: Must be after DenseLinAlgPack_LinAlgOpPack.hpp
00051                                                                         // for Intel C++ 5.0 (Windows 2000).  Namespaces
00052                                                                         // and lookup rules don't work properly with this compiler.
00053 #include "Teuchos_FancyOStream.hpp"
00054 
00055 namespace {
00056 
00057 // Get a view of a vector (two versions)
00058 
00059 inline
00060 Teuchos::RCP<const AbstractLinAlgPack::Vector>
00061 get_view(
00062   const AbstractLinAlgPack::Vector  &v
00063   ,const RangePack::Range1D         &rng
00064   )
00065 {
00066    return v.sub_view(rng);
00067 }
00068 
00069 inline
00070 Teuchos::RCP<const AbstractLinAlgPack::SpVectorSlice>
00071 get_view(
00072   const AbstractLinAlgPack::SpVectorSlice &v
00073   ,const RangePack::Range1D               &rng
00074   )
00075 {
00076   return Teuchos::rcp( new AbstractLinAlgPack::SpVectorSlice( v(rng) ) );
00077 }
00078 
00079 // Matrix-vector multiplication
00080 
00081 template<class V>
00082 void mat_vec(
00083   AbstractLinAlgPack::VectorMutable        *y
00084   ,AbstractLinAlgPack::value_type          a
00085   ,AbstractLinAlgPack::value_type          alpha
00086   ,const AbstractLinAlgPack::MatrixOp      &D
00087   ,BLAS_Cpp::Transp                        D_trans
00088   ,const DenseLinAlgPack::Range1D          &D_rng
00089   ,const DenseLinAlgPack::Range1D          &I_rng
00090   ,BLAS_Cpp::Transp                        M_trans
00091   ,const V                                 &x
00092   ,AbstractLinAlgPack::value_type          b
00093   )
00094 {
00095   using BLAS_Cpp::no_trans;
00096   using BLAS_Cpp::trans;
00097   using BLAS_Cpp::trans_not;
00098 
00099   //
00100   // M = [ alpha*op(D) ] or [      I       ]
00101   //     [     I       ]    [ alpha*op(D)  ]
00102   //
00103   if( M_trans == no_trans ) {
00104     // 
00105     // y = b*y + a*M*x
00106     // =>
00107     // y(D_rng) = b*y(D_rng) + a*alpha*op(D)*x
00108     // y(I_rng) = b*y(I_rng) + a*x
00109     //
00110     AbstractLinAlgPack::VectorSpace::vec_mut_ptr_t
00111       y_D = y->sub_view(D_rng),
00112       y_I = y->sub_view(I_rng);
00113     // y(D_rng) = b*y(D_rng) + a*alpha*op(D)*x
00114     AbstractLinAlgPack::Vp_StMtV( y_D.get(), a*alpha, D, D_trans, x, b );
00115     // y(I_rng) = b*y(I_rng) + a*x
00116     LinAlgOpPack::Vt_S( y_I.get(), b );
00117     LinAlgOpPack::Vp_StV( y_I.get(), a, x );
00118   }
00119   else {
00120     //
00121     // y = b*y + a*M'*x
00122     // =>
00123     // y = b*y + a*alpha*op(D')*x(D_rng) + a*x(I_rng)
00124     //
00125     LinAlgOpPack::Vp_StMtV( y, a*alpha, D, trans_not(D_trans), *get_view(x,D_rng), b );
00126     LinAlgOpPack::Vp_StV( y, a, *get_view(x,I_rng) );
00127   }
00128 }
00129 
00130 } // end namespace
00131 
00132 namespace ConstrainedOptPack {
00133 
00134 // Overridden from MatrixBase
00135 
00136 size_type MatrixIdentConcat::rows() const
00137 {
00138   return this->D_rng().size() + this->I_rng().size();
00139 }
00140 
00141 size_type MatrixIdentConcat::cols() const
00142 {
00143   return this->I_rng().size();
00144 }
00145 
00146 size_type MatrixIdentConcat::nz() const
00147 {
00148   const MatrixOp& D = this->D();
00149   return D.nz() + BLAS_Cpp::cols(D.rows(),D.cols(),D_trans()); // D and I
00150 }
00151 
00152 // Overridden from MatrixOp
00153 
00154 std::ostream& MatrixIdentConcat::output(std::ostream& out_arg) const
00155 {
00156   Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00157   Teuchos::OSTab tab(out);
00158   const Range1D           D_rng   = this->D_rng();
00159   const BLAS_Cpp::Transp  D_trans = this->D_trans();
00160   *out << "Converted to dense =\n";
00161   MatrixOp::output(*out);
00162   *out << "This is a " << this->rows() << " x " << this->cols()
00163     << " general matrix / identity matrix concatenated matrix object ";
00164   if( D_rng.lbound() == 1 ) {
00165     if( D_trans == BLAS_Cpp::no_trans )
00166       *out << "[ alpha*D; I ]";
00167     else
00168       *out << "[ alpha*D'; I ]";
00169   }
00170   else {
00171     if( D_trans == BLAS_Cpp::no_trans )
00172       *out << "[ I; alpha*D ]";
00173     else
00174       *out << "[ I; alpha*D' ]";
00175   }
00176   *out << " where alpha and D are:";
00177   tab.incrTab();
00178   *out << "\nalpha = " << this->alpha();
00179   *out << "\nD =\n" << D();
00180   return out_arg;
00181 }
00182 
00183 void MatrixIdentConcat::Vp_StMtV(
00184   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00185   , const Vector& x, value_type b
00186   ) const
00187 {
00188   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,M_trans,x);
00189   mat_vec( y, a, alpha(), D(), D_trans(), D_rng(), I_rng(), M_trans, x, b );
00190 }
00191 
00192 void MatrixIdentConcat::Vp_StMtV(
00193   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00194   , const SpVectorSlice& x, value_type b
00195   ) const
00196 {
00197 
00198   AbstractLinAlgPack::Vp_MtV_assert_compatibility(y,*this,M_trans,x);
00199   mat_vec( y, a, alpha(), D(), D_trans(), D_rng(), I_rng(), M_trans, x, b );
00200 }
00201 
00202 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines