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

Generated on Tue Jul 13 09:30:51 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7