ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixIdentConcat.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
43 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
44 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
45 #include "AbstractLinAlgPack_VectorSpace.hpp"
46 #include "AbstractLinAlgPack_VectorMutable.hpp"
47 #include "AbstractLinAlgPack_VectorStdOps.hpp"
48 #include "AbstractLinAlgPack_SpVectorClass.hpp"
49 #include "AbstractLinAlgPack_AssertOp.hpp"
50 #include "AbstractLinAlgPack_SpVectorOp.hpp" // RAB: 12/20/2002: Must be after DenseLinAlgPack_LinAlgOpPack.hpp
51  // for Intel C++ 5.0 (Windows 2000). Namespaces
52  // and lookup rules don't work properly with this compiler.
53 #include "Teuchos_FancyOStream.hpp"
54 
55 namespace {
56 
57 // Get a view of a vector (two versions)
58 
59 inline
61 get_view(
63  ,const RangePack::Range1D &rng
64  )
65 {
66  return v.sub_view(rng);
67 }
68 
69 inline
71 get_view(
72  const AbstractLinAlgPack::SpVectorSlice &v
73  ,const RangePack::Range1D &rng
74  )
75 {
76  return Teuchos::rcp( new AbstractLinAlgPack::SpVectorSlice( v(rng) ) );
77 }
78 
79 // Matrix-vector multiplication
80 
81 template<class V>
82 void mat_vec(
84  ,AbstractLinAlgPack::value_type a
85  ,AbstractLinAlgPack::value_type alpha
87  ,BLAS_Cpp::Transp D_trans
88  ,const DenseLinAlgPack::Range1D &D_rng
89  ,const DenseLinAlgPack::Range1D &I_rng
90  ,BLAS_Cpp::Transp M_trans
91  ,const V &x
92  ,AbstractLinAlgPack::value_type b
93  )
94 {
95  using BLAS_Cpp::no_trans;
96  using BLAS_Cpp::trans;
97  using BLAS_Cpp::trans_not;
98 
99  //
100  // M = [ alpha*op(D) ] or [ I ]
101  // [ I ] [ alpha*op(D) ]
102  //
103  if( M_trans == no_trans ) {
104  //
105  // y = b*y + a*M*x
106  // =>
107  // y(D_rng) = b*y(D_rng) + a*alpha*op(D)*x
108  // y(I_rng) = b*y(I_rng) + a*x
109  //
111  y_D = y->sub_view(D_rng),
112  y_I = y->sub_view(I_rng);
113  // y(D_rng) = b*y(D_rng) + a*alpha*op(D)*x
114  AbstractLinAlgPack::Vp_StMtV( y_D.get(), a*alpha, D, D_trans, x, b );
115  // y(I_rng) = b*y(I_rng) + a*x
116  LinAlgOpPack::Vt_S( y_I.get(), b );
117  LinAlgOpPack::Vp_StV( y_I.get(), a, x );
118  }
119  else {
120  //
121  // y = b*y + a*M'*x
122  // =>
123  // y = b*y + a*alpha*op(D')*x(D_rng) + a*x(I_rng)
124  //
125  LinAlgOpPack::Vp_StMtV( y, a*alpha, D, trans_not(D_trans), *get_view(x,D_rng), b );
126  LinAlgOpPack::Vp_StV( y, a, *get_view(x,I_rng) );
127  }
128 }
129 
130 } // end namespace
131 
132 namespace ConstrainedOptPack {
133 
134 // Overridden from MatrixBase
135 
137 {
138  return this->D_rng().size() + this->I_rng().size();
139 }
140 
142 {
143  return this->I_rng().size();
144 }
145 
147 {
148  const MatrixOp& D = this->D();
149  return D.nz() + BLAS_Cpp::cols(D.rows(),D.cols(),D_trans()); // D and I
150 }
151 
152 // Overridden from MatrixOp
153 
154 std::ostream& MatrixIdentConcat::output(std::ostream& out_arg) const
155 {
156  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
157  Teuchos::OSTab tab(out);
158  const Range1D D_rng = this->D_rng();
159  const BLAS_Cpp::Transp D_trans = this->D_trans();
160  *out << "Converted to dense =\n";
161  MatrixOp::output(*out);
162  *out << "This is a " << this->rows() << " x " << this->cols()
163  << " general matrix / identity matrix concatenated matrix object ";
164  if( D_rng.lbound() == 1 ) {
165  if( D_trans == BLAS_Cpp::no_trans )
166  *out << "[ alpha*D; I ]";
167  else
168  *out << "[ alpha*D'; I ]";
169  }
170  else {
171  if( D_trans == BLAS_Cpp::no_trans )
172  *out << "[ I; alpha*D ]";
173  else
174  *out << "[ I; alpha*D' ]";
175  }
176  *out << " where alpha and D are:";
177  tab.incrTab();
178  *out << "\nalpha = " << this->alpha();
179  *out << "\nD =\n" << D();
180  return out_arg;
181 }
182 
184  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
185  , const Vector& x, value_type b
186  ) const
187 {
189  mat_vec( y, a, alpha(), D(), D_trans(), D_rng(), I_rng(), M_trans, x, b );
190 }
191 
193  VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
194  , const SpVectorSlice& x, value_type b
195  ) const
196 {
197 
199  mat_vec( y, a, alpha(), D(), D_trans(), D_rng(), I_rng(), M_trans, x, b );
200 }
201 
202 } // end namespace ConstrainedOptPack
virtual const MatrixOp & D() const =0
void Vp_StMtV(VectorMutable *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &vs_rhs2, value_type beta) const
virtual vec_mut_ptr_t sub_view(const Range1D &rng)
vec_ptr_t sub_view(const index_type &l, const index_type &u) const
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
basic_OSTab< CharT, Traits > & incrTab(const int tabs=1)
T * get() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Range1D I_rng() const =0
size_t size_type
virtual value_type alpha() const =0
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
virtual Range1D D_rng() const =0
Transp trans_not(Transp _trans)
virtual BLAS_Cpp::Transp D_trans() const =0
void Vp_MtV_assert_compatibility(VectorMutable *v_lhs, const MatrixOp &m_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)